Anzeige

C: Einträge in Protein Data Bank (PDB) - Files selektiv auslesen


JuK

Grünschnabel
#1
Guten Tag liebe Gemeinde,

bin ganz neu hier im Forum und auch mit C erst seit Anfang April am werkeln. Das Ganze findet im Rahmen eines Wahlpflichtmoduls meines Biochemiestudiums statt, entsprechend dünn ist das Fundament an Kenntnissen und Erfahrungen.

Aufgabe:

Es gilt ein Programm zu schreiben, das die Torsionswinkel zwischen den Atomen in Nukleinsäuren berechnet, wobei dynamische Speicheralloziierung stattfinden soll. Die Rohdaten dafür werden auch aus PDB-files bezogen, wobei es sich im Grunde nur um Textdateinen handelt (siehe Anhang).

Umstand:

Erst ab Zeile 315 beginnt die relevante Auflistung von den, an der Struktur beteiligten, Atomen. Das sieht dann folgendermaßen aus:

ATOM 1 O5' G A 1 40.972 8.748 5.279 1.00 29.42 O
ATOM 2 C5' G A 1 42.038 7.963 5.819 1.00 25.67 C
ATOM 3 C4' G A 1 41.482 6.699 6.422 1.00 24.13 C
ATOM 4 O4' G A 1 40.978 6.977 7.756 1.00 23.86 O
ATOM 5 C3' G A 1 40.286 6.109 5.685 1.00 21.93 C
ATOM 6 O3' G A 1 40.700 5.292 4.598 1.00 19.64 O
ATOM 7 C2' G A 1 39.606 5.296 6.780 1.00 22.18 C
ATOM 8 O2' G A 1 40.213 4.045 7.016 1.00 22.19 O
ATOM 9 C1' G A 1 39.814 6.197 7.997 1.00 21.83 C
ATOM 10 N9 G A 1 38.695 7.104 8.209 1.00 18.54 N
ATOM 11 C8 G A 1 38.707 8.473 8.103 1.00 17.96 C
ATOM 12 N7 G A 1 37.551 9.017 8.364 1.00 18.23 N
ATOM 13 C5 G A 1 36.724 7.941 8.655 1.00 17.17 C
ATOM 14 C6 G A 1 35.362 7.912 9.019 1.00 15.25 C
ATOM 15 O6 G A 1 34.602 8.855 9.166 1.00 15.04 O
ATOM 16 N1 G A 1 34.909 6.609 9.225 1.00 15.08 N
ATOM 17 C2 G A 1 35.681 5.478 9.105 1.00 16.78 C
ATOM 18 N2 G A 1 35.070 4.308 9.370 1.00 14.87 N
ATOM 19 N3 G A 1 36.961 5.496 8.760 1.00 16.30 N
ATOM 20 C4 G A 1 37.415 6.753 8.556 1.00 17.87 C
ATOM 21 P U A 2 39.955 5.413 3.180 1.00 19.59 P
ATOM 22 OP1 U A 2 40.555 4.430 2.256 1.00 20.04 O
ATOM 23 OP2 U A 2 39.887 6.839 2.784 1.00 20.94 O
ATOM 24 O5' U A 2 38.472 4.919 3.478 1.00 19.95 O
.
.
.
(Die fett markierten Einträge sind für das Lösen der Aufgabe erforderlich)

Ziel:
Erstes großes Ziel wäre die selektive Formatierung der Daten in eine übergeordnete Struktur, indem quasi folgende sprachlich formulierte Aussage C-kryptisch umgesetzt wird:

Falls das erste Wort einer Zeile "ATOM" lautet, sollen die Einträge 2-9 dieser Zeile in die spezifischen Einträge der Struktur "Molecule" geschrieben werden...

Ansatz:
... dazu dachte ich wäre fscanf erstmal sehr nützlich. Das Ganze habe ich mit zwei unendlichen While-Schleifen verknüpft, wobei die erste nur dazu dient "numcount" zu bestimmen. Das wäre dann die Anzahl an ATOM-Zeilen, die dann für die erste For-Schleife wichitg wäre. Diese befindet sich in der zweiten While-Schleife.

Den Kot den man bis jetzt produziert hat sieht so aus:

C++:
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>

#define NATOM 500

/*================================================================*/
typedef struct 
{
   int line;
   int num;
   char ATOM[NATOM];
   int serialnum[NATOM];
   char recname[NATOM];
   char residue[NATOM];
   char chain[NATOM];
   int residueserialnum[NATOM];
   double xyz[NATOM][3]; 
}Molecule;

/*================================================================*/
void check_args(int argc)
{
   if(argc != 2)
   {
      printf("Vorgegebene Eingabestruktur:\n");
      printf("[Programmname] [Inputdatei]\n\n");
   }
}

/*================================================================*/
void read_molecule(Molecule *mol, char *fname) {

   FILE *fpin;
   char buff[5];
   int i, numcount;

    if ( (fpin=fopen(fname, "r")) == NULL) {
           printf("Datei %s konnte nicht aufgerufen werden. ABBRUCH.\n\n", fname);     
           exit(EXIT_FAILURE);
       };

   while(1){
    fscanf(fpin, "%s", buff);
    if (buff == "ATOM"){
    numcount++;
    }
    if(buff == "END"){
        break;
    }
   printf("%d", numcount);

   }

   while (1){
    fscanf(fpin, "%s", (*mol).ATOM);
    if ((*mol).ATOM == "ATOM"){
        for(i=0; i<=(numcount+1); i++){ /*numcount+1, da die Auflistung der Atome durch eine TER-Zeile (line 482) getrennt wird*/
        
        fscanf(fpin, "%s","%d","%s","%s","%s","%d","%d","%d","%d", (*mol).ATOM[i], &(*mol).serialnum[i], (*mol).recname[i], (*mol).residue[i], (*mol).chain[i], &(*mol).residueserialnum[i], &(*mol).xyz[i][0], &(*mol).xyz[i][1], &(*mol).xyz[i][2]);      

            if((*mol).ATOM == "END"){
                break;
            }
        }
    }
   break;  
   }

   fclose(fpin);
}

/*================================================================*/
   int main(int argc, char *argv[])
   {
      Molecule mol;

    check_args (argc);
       read_molecule (&mol, argv[1]);

       return 0;
   }
Problem:
Auf Compilerisch:

warning: writing into constant object (argument 3) [-Wformat=]
fscanf(fpin, "%s","%d","%s","%s","%s","%d","%d","%d","%d", (*mol).ATOM, &(*mol).serialnum, (*mol).
^
warning: too many arguments for format [-Wformat-extra-args]
fscanf(fpin, "%s","%d","%s","%s","%s","%d","%d","%d","%d", (*mol).ATOM, &(*mol).serialnum, (*mol).

--> Wobei es bereits schön wäre wenn das Programm in seiner Ausführung soweit kommen würde, denn es hängt in der ersten While - Schleife fest, was man mit dem printf (line 54) prüfte.

Erkenntnisse:

Wie es so ist; während man das Problem erklärt fällt einem Einiges auf:
--> line 48-53:
C++:
    if (buff == "ATOM"){
    numcount++;
    }
    if(buff == "END"){
        break;
    }
wird nicht ausgeführt, sonst wäre numcount nicht 0 && der break müsste stattfinden.
Jetzt lese ich gerade, dass es an dem Unterschied zwischen String ("ATOM") und Char (buff) liegen könnte. Des Weiteren funktioniert der break in der while Schleife nicht, "return 0" hingegen schon, was mich dann halt nur aus der ganzen Funktion wirft, was natürlich zu vermeiden wäre.

--> Da Kette A und B in der PDB-Datei durch eine TER-Zeile (ln 482) getrennt sind könnte das bestimmt Probleme mit fscanf in der While - Schleife geben, da diese Schleife den ATOM-Block ja wie eine kontinuierliche Liste behandeln sollte.

An Euch:
Für jede Art von Hilfe, Inspiration, Tipps oder Erfahrungsaustausch wäre ich sehr dankbar, da ich mit meinem Latein adhuc am Ende bin.

Herzliche Grüße,
JuK
 

Anhänge

cwriter

Erfahrenes Mitglied
#2
Hi

warning: writing into constant object (argument 3) [-Wformat=]
fscanf(fpin, "%s","%d","%s","%s","%s","%d","%d","%d","%d", (*mol).ATOM, &(*mol).serialnum, (*mol).
Dies liegt an einem offensichtlichen Missverständnis darüber, wie fscanf funktioniert.
Das 2. Argument ist der Formatstring, bei dir nun "%s". Und das heisst: Lies einen String bis zum ersten Whitespace. Wohin? Nun, in das nächste Argument. Das nächste Argument ist aber "%d", und das ist ein String Literal, also ein "const char* const" (konstante Zeichen in einer nicht veränderbaren Adresse).
Daher die Warnung: "Writing into constant object".
Fix: All deine "%s", "%d" usw. müssen in einem String Literal sein, mit Leerzeichen abgetrennt.
Etwa so:
C:
fscanf(fpin, "%s %d %s %s %s %f %f %f %f", (*mol).ATOM[I], &(*mol).serialnum[I], (*mol).[/I][/I]
Ein kleiner Fix dazu: Um Gleitkommazahlen zu lesen, nimmt man statt "%d" (für Ganzzahlen) ein "%f".

Problem 2:
Dein buff ist definiert als
C:
char buf[5];
Aber was bedeutet das denn genau? Das bedeutet, dass du nach FILE* fpin (was Speicher braucht) weitere Variablen hinzufügst, nämlich genau 5 char-Elemente.
Also was sollte dann "==" vergleichen?
C ist eine dumme Sprache. Sie kennt eigentlich keine Strings. Für sie ist jede Art von Daten einfach nur eine Speicheradresse (mit einigen wenigen Ausnahmen).
Deshalb kannst du auch von fast jeder Variable herausfinden, wo sie ist:
C:
int i = 0;
int* wo = &i;
printf("i ist an Adresse: %p\n", wo);
Nun ist es mit Arrays ähnlich. Du hast einen Array. Wie kann man den am besten repräsentieren?
C wählte einfach Startadresse und Länge. Das reicht, den Array exakt zu kennen.
Nun sind das aber 2 Werte, und das ist schlecht zu verarbeiten. Also schmeisst man in vielen Operationen einfach die Länge weg. Das ist einfach, aber auch einer der vielen Gründe, warum C als "unsicher" angesehen wird: Man kann viele Fehler einbauen.
Ok. Du vergleichst nun
C:
buff == "ATOM"
buff wird zu einem char* umgewandelt. Dieser Pointer zeigt nun auf den Start des Arrays und hat keine Längeninformation mehr.
Das kann gar nie true sein. buff ist auf dem Stack, dein String-Literal wird normalerweise in ein Read-Only-Segment gelinkt. Der einfachste Weg, sich das zu merken: Es müssen ja andere Werte da reinpassen. Also kann nicht ein String Literal an derselben Adresse stehen (sonst könnte man den ja überschreiben, was ähnlich wäre wie 1 := 2)
Ok, viel Theorie. Die Lösung: In C vergleicht man Strings mit strcmp() oder strncmp().
In deinem Fall:
C:
if(strcmp(buff, "ATOM") == 0)
{

}
else if(strcmp(buff, "END") == 0)
{

}
strcmp gibt 0 zurück, wenn die Strings identisch sind, -1, wenn der erste vor dem 2. kommt (lexikographisch geordnet) und sonst 1.

Deine Art, die betroffenen Zeilen zu finden, ist etwas seltsam. Nutze doch stattdessen fgets(), um solange Zeilen zu lesen, bis strstr() sagt, dass die Zeile mit "ATOM" beginnt. fscanf() mit "%s" ist dafür eher nicht geeignet, da es ja nicht unbedingt nur einen String pro Zeile liest.

Andere Tips:
"(*x)." ist identisch zu "x->", letzteres aber leichter zu lesen.

Dein fscanf()-format ist ein bisschen seltsam. Du kannst durchaus auch sowas wie "%d.%d" für deine Fixed-Point Numbers nehmen. Damit wäre es wohl ein bisschen einfacher zu lesen und es erlaubt es, den Rückgabewert besser zur Fehlerprüfung zu nutzen. Das Problem bei Doubles ist die Präzision, es könnte zu Rundungsfehlern kommen. Wenn das gemäss Aufgabe ok ist, musst du das alles aber nicht selber basteln und kannst ganz einfach %f nehmen.

Übrigens: Eine sehr schön strukturierte Frage mit allem, was zum Helfen erforderlich ist. Sehr nett! :)

Gruss
cwriter
 

JuK

Grünschnabel
#3
Vielen Dank cwriter!

Nach ein paar Umstrukturierungen sieht das wesentliche Codefragment (ln 48 - 89) nun folgendermaßen aus:

C++:
/*================================================================*
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>

#define LEN 1000
#define NATOM 1000

/*================================================================*/
typedef struct 
{
   int line;
   int num[NATOM];
   char ATOM[NATOM][5];
   int serialnum[NATOM];
   char recname[NATOM][5];
   char residue[NATOM][5];
   char chain[NATOM][2];
   int residueserialnum[NATOM];
   float xyz[NATOM][3]; 
}Molecule;

/*================================================================*/
void check_args(int argc)
{
   if(argc != 2)
   {
      printf("Vorgegebene Eingabestruktur:\n");
      printf("[Programmname] [Inputdatei]\n\n");
   }
}

/*================================================================*/
int read_molecule(Molecule *mol, char *fname) {

   FILE *fpin;
   char line[LEN][10];
   int i=0, t, buff, numcount, linecount, endline;

    if ( (fpin=fopen(fname, "r")) == NULL) {
           printf("Datei %s konnte nicht aufgerufen werden. ABBRUCH.\n\n", fname);     
           exit(EXIT_FAILURE);
       };


   while(1) 
   {
      buff = fgetc(fpin);
      if (buff == '\n')  
      {
         linecount++;
      }

      else if (buff == EOF)
      {
         break;   
      }   
   }
   printf("Anzahl der Zeilen:                         %d\n", linecount);
   
   fseek(fpin, 0L, SEEK_SET);

   for (t=0; t<linecount; t++)
   {

      fgets(line[t], LEN, fpin);
      fscanf(fpin, "%s", (*mol).ATOM[t]);
   
      if (strcmp((*mol).ATOM[t], "ATOM") == 0)
      {
         numcount++;
         fscanf(fpin, "%d %s %s %s %d %f %f %f", 
             &(*mol).serialnum[i],
            (*mol).recname[i], (*mol).residue[i], (*mol).chain[i], 
            &(*mol).residueserialnum[i], &(*mol).xyz[i][0], &(*mol).xyz[i][1], 
            &(*mol).xyz[i][2]);      
         i++;
      }

   }
   printf("Anzahl ATOM als erstes Wort einer Zeile:   %d\n", numcount);
   
   for (i=0; i<numcount; i++)
   {
      printf(" %d %s %s %s %d %f %f %f \n", (*mol).serialnum[i],
         (*mol).recname[i], (*mol).residue[i], (*mol).chain[i], (*mol).residueserialnum[i], (*mol).xyz[i][0], (*mol).xyz[i][1], (*mol).xyz[i][2]);
   }

   fclose(fpin);
}

/*================================================================*/
   int main(int argc, char *argv[])
   {
      Molecule mol;

        check_args (argc);
       read_molecule (&mol, argv[1]);

       return 0;
   }
"(*x)." ist identisch zu "x->", letzteres aber leichter zu lesen.
Stimmt. Wurde aber nicht beherzigt, weil der Prof da etwas eigen ist :/

Die While-Schleife ermittelt erfolgreich die Zeilenanzahl, um die nächste For-Schleife in Ihrer Ausführung zu begrenzen, wobei ich auf die strstr() - Funktion verzichtet habe, da es zu einer seltsamen mir unerklärlichen Änderung des Outputs kommt.

Zudem wurde ein zweites Inkrement eingeführt (Zeile 79), damit das Auswertungsarray (Molecule) mit 1 und nicht mit 315 anfängt.

Auf Alternativen von fscanf habe ich ehrlich gesagt nicht getestet, geschweige denn danach gesucht.
Dabei ist mir aber eine interessante Entdeckung unter gekommen:
Das zweite fscanf (Zeile 74) funktioniert nur tadelfrei, falls man auf den Scan des ersten Wortes (ATOM) verzichtet:

fscanf(fpin, "%d %s %s %s %d %f %f %f",
&(*mol).serialnum, etc.

anstatt

fscanf(fpin, "%s %d %s %s %s %d %f %f %f",
(*mol).ATOM , &(*mol).serialnum, etc.


und ich habe nicht die geringste Idee warum.

Fragen:
"Kommunizieren" diese scan-Befehle untereinander, bzw. arbeiten sie kooperativ? So nach dem Motto: "Das wurde bereits gescannt. Mach was anderes!" . Oder wie kann man sich das vorstellen ?

Ist es möglich das "#define NATOM 1000" durch eine Variable zu ersetzen, die erst in einer Funktion bestimmt wird? Ziel wäre es hier beim Verwenden größerer PDB-Files nicht immer wieder nachzuschauen, ob der Wert für NATOM ausreicht um die gesamten Daten zu fassen.

Wie im ersten Beitrag erwähnt soll das Ganze am Ende dynamischer Speicheralloziierung folgen, wobei ich (noch) gar keine Ahnung habe wie sich das gestalten soll. Falls Euch jetzt schon was auf- oder einfällt, das Ihr bereit seit zu teilen wäre das superkalifragilistiexpialigetisch.

Einstweilen viele liebe Grüße,
JuK

P.S. und weiß jemand zufällig warum bei mir irgendwann alles kursiv wird ?
 

cwriter

Erfahrenes Mitglied
#4
Die While-Schleife ermittelt erfolgreich die Zeilenanzahl, um die nächste For-Schleife in Ihrer Ausführung zu begrenzen, wobei ich auf die strstr() - Funktion verzichtet habe, da es zu einer seltsamen mir unerklärlichen Änderung des Outputs kommt.
Das ist eigentlich nicht nötig und vor allem langsam (fgetc() liest ein Zeichen pro Aufruf...).
fgets() kann auch direkt Zeilen lesen und zählen. (fgets() liest ja gerade immer eine Zeile, sofern der Buffer gross genug ist. Wenn dieser jenes nicht ist, dann kannst du das ja auch behandeln (prüfe auf letztes Zeichen. Ist es '\n', so wurde das Ende der Zeile erreicht). Allerdings sollte ein 1024 Buffer für die allermeisten Textdateien reichen.

Die While-Schleife ermittelt erfolgreich die Zeilenanzahl, um die nächste For-Schleife in Ihrer Ausführung zu begrenzen, wobei ich auf die strstr() - Funktion verzichtet habe, da es zu einer seltsamen mir unerklärlichen Änderung des Outputs kommt.
strstr() verändert die Strings nicht (wie auch? Die Strings selbst sind const und die Pointer können innerhalb der Funktion nicht geändert werden).

Auf Alternativen von fscanf habe ich ehrlich gesagt nicht getestet, geschweige denn danach gesucht.
Dabei ist mir aber eine interessante Entdeckung unter gekommen:
Das zweite fscanf (Zeile 74) funktioniert nur tadelfrei, falls man auf den Scan des ersten Wortes (ATOM) verzichtet:
Naja, mit dem fscanf Zeile 69 hast du ja schon einen String konsumiert und der Cursor der Datei wurde entsprechend bewegt. Normalerweise will man dasselbe ja nicht 2x lesen :)
Die fscanfs kommunizieren nicht direkt aufeinander, aber sie arbeiten ja auf derselben Datei. Der State der Datei (nämlich die Cursorposition) ändert sich nach solchen Aufrufen, daher änderst du ja quasi indirekt einen Parameter (nämlich einen Wert, der im FILE-struct liegt (bzw. dahinter).

Ist es möglich das "#define NATOM 1000" durch eine Variable zu ersetzen, die erst in einer Funktion bestimmt wird? Ziel wäre es hier beim Verwenden größerer PDB-Files nicht immer wieder nachzuschauen, ob der Wert für NATOM ausreicht um die gesamten Daten zu fassen.
Möglich ist alles, nur einfach wird es nicht. Naja, geht schon: Du bräuchtest einen Konstruktor für jedes Struct-Object. Etwa so:
C:
typedef struct 
{
   int line;
   int* num;
   char** ATOM;
   int* serialnum;
   char** recname;
   char** residue;
   char** chain;
   int* residueserialnum;
   float** xyz; 
}Molecule;


bool initMolecule(Molecule* m, size_t NATOM)
{
   // #TODO: Was, wenn calloc() fehlschlägt?
   m->ATOM = calloc(sizeof(char*), NATOM);
   for(size_t i = 0; i < NATOM; ++i) m->ATOM[i] = calloc(sizeof(char), 5);

   //Ähnlich für die anderen char** / float**

   m->serialnum = calloc(sizeof(int), NATOM);

   //Ähnlich für num, residueserialnum
}

void freeMolecule(Molecule* m, size_t NATOM)
{
   for(size_t i = 0; i < NATOM; ++i) free(m->ATOM[i]);
   free(m->ATOM);

   //usw
}
Wie du siehst, gibt es hier einige Probleme, z.B. muss NATOM auch in der free-Operation bekannt sein. Eine Lösung ist, NATOM einfach auch im Struct zu speichern. Übrigens könnte man die **-Variablen durchaus vermeiden, dann würde sich aber die Syntax im restlichen Code ändern. (Man könnte es mit einem eindimensionalen Pointer machen, dann wäre aber z.B.
C:
ATOM[i][j];
//wird zu
ATOM[i * NATOM + j];
)

P.S. und weiß jemand zufällig warum bei mir irgendwann alles kursiv wird ?
Das liegt am BB-Code für kursiv [i]bla[/i].
In C im allgemeinen und in deinem Code im Speziellen kommt aber oft array[i] vor => der Index verschwindet und alles danach wird kursiv. Einfachster Fix: Einfach jeden Code in Codetags packen. Anderer Fix: [noparse][/ noparse] (ohne Abstand von [ zu /) verwenden, dann wird es nicht mehr interpretiert

Gruss
cwriter
 
Anzeige
Anzeige