Commenting each function in progress. Correcting a bug for the HCMD2 data read

src/INTBuilder.c

@@ -15,10 +15,13 @@ void getCandidatesForP1(struct pdb_values* pdb2, struct residue** t_candid1, str
     /* This function will update the t_candid1 array including all candidate residues
      * of P1. At the end of the function, the first *nbCand1 residues in t_candid1 will
      * correspond to P1 candidates
+     *
+     * This corresponds to the algorithm as described in the getInterface function
      */
 
     float x = 0, y = 0, z = 0;
     float centerX = 0, centerY = 0, centerZ = 0;
+    float farthestX = 0, farthestY = 0, farthestZ = 0;
     float dist = 0, minDist = -1, maxDist = -1, threshold_dist = 0;
     int i = 0, j = 0, idxCand = -1, idxAtom = -1;
     int nbCandInit1 = *nbCand1;
@@ -37,7 +40,6 @@ void getCandidatesForP1(struct pdb_values* pdb2, struct residue** t_candid1, str
             z = t_candid1[i]->z[j];
 
             dist = sqrt((x-centerX)*(x-centerX) + (y-centerY)*(y-centerY) + (z-centerZ)*(z-centerZ));
-                //printf("dist %f\n",dist);
             if(dist > maxDist){
                 maxDist = dist;
                 idxAtom = j;
@@ -48,25 +50,25 @@ void getCandidatesForP1(struct pdb_values* pdb2, struct residue** t_candid1, str
 
     /* Then we get the closest distance from that residue to a residue of the receptor */
 
-    /* Not center of P1, but farthest residue from P2 */
-    centerX = t_candid1[idxCand]->x[idxAtom];
-    centerY = t_candid1[idxCand]->y[idxAtom];
-    centerZ = t_candid1[idxCand]->z[idxAtom];
+    /* farthestX,Y,Z represent the farthest residue of P1 from P2 */
+    farthestX = t_candid1[idxCand]->x[idxAtom];
+    farthestY = t_candid1[idxCand]->y[idxAtom];
+    farthestZ = t_candid1[idxCand]->z[idxAtom];
     for(i=0;i<nbCandInit2;i++){
         for(j=0;j<t_candid2[i]->nbAtom;j++){
             x = t_candid2[i]->x[j];
             y = t_candid2[i]->y[j];
             z = t_candid2[i]->z[j];
 
-            dist = sqrt((x-centerX)*(x-centerX) + (y-centerY)*(y-centerY) + (z-centerZ)*(z-centerZ));
+            dist = sqrt((x-farthestX)*(x-farthestX) + (y-farthestY)*(y-farthestY) + (z-farthestZ)*(z-farthestZ));
             if(dist < minDist || minDist < 0){
                 minDist = dist;
             }
         }
     }
 
-    /* We remove every residue on pdb1 that is closer than treshold_dist from the center
-     * of mass
+    /* We remove every residue on pdb1 that is closer than treshold_dist from the farthest
+     * residue computed just before
      */
     threshold_dist = minDist - DIST_FOR_CONTACT;
     for(i=0;i<nbCandInit1;i++){
@@ -83,7 +85,7 @@ void getCandidatesForP1(struct pdb_values* pdb2, struct residue** t_candid1, str
             z = t_candid1[i]->z[j];
 
 
-            dist = sqrt((x-centerX)*(x-centerX) + (y-centerY)*(y-centerY) + (z-centerZ)*(z-centerZ));
+            dist = sqrt((x-farthestX)*(x-farthestX) + (y-farthestY)*(y-farthestY) + (z-farthestZ)*(z-farthestZ));
             if(dist > threshold_dist){
                 t_candid1[*nbCand1] = t_candid1[i];
                 (*nbCand1)++;
@@ -94,6 +96,22 @@ void getCandidatesForP1(struct pdb_values* pdb2, struct residue** t_candid1, str
 }
 
 void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
+    /* This function will output the docking interface between pdbR and pdbL 
+     * To compute the docking interface, we considered that any residue of P1 at or closer than 10A
+     * of a residue of P2 will represent a contact
+     *
+     * To avoid unecessary computations, we remove residues that can't be in the interface given the
+     * following algorithm :
+     * Take a point of P1
+     * let point2 be the farthest residue of P2 from the selected point in P1
+     * let d1 be the distance between point2 and its closest point of P1. 
+     * let d2 be the substraction to d1 the contact distance threshold (10A in this case)
+     * Remove all residues of P2 that are strictly closer than d2 to the selected point point2
+     *
+     * Iterate that algorithm multiple time over P1 and P2. This allows to greatly narrow down the 
+     * final quadratic complexity
+     */
+
     int i = 0, j = 0, k = 0, l = 0;
     int oldCandL = -1, oldCandR = -1;
     float x1 = 0, y1 = 0, z1 = 0;
@@ -116,19 +134,19 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
         exit(EXIT_FAILURE);
     }
 
-    //memset(t_candidateL,0,(1+pdbL->nbRes)*sizeof(struct residue));
-    //memset(t_candidateR,0,(1+pdbR->nbRes)*sizeof(struct residue));
-
+    /* Initialize all the potential candidates */
     for(i=0;i<pdbL->nbRes;i++)
         t_candidateL[i] = &pdbL->residues[i];
     for(i=0;i<pdbR->nbRes;i++)
         t_candidateR[i] = &pdbR->residues[i];
 
+    /* Start removing candidates */
     while(oldCandL != nbCandL && oldCandR != nbCandR){
-        k++;
         oldCandL = nbCandL;
         oldCandR = nbCandR;
         getCandidatesForP1(pdbR,t_candidateL,t_candidateR,&nbCandL,&nbCandR);
+
+        /* If there is no candidate remaining */
         if(nbCandL == 0 || nbCandR == 0){
             nbCandL = 0;
             nbCandR = 0;
@@ -136,17 +154,24 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
         }
 
         getCandidatesForP1(pdbL,t_candidateR,t_candidateL,&nbCandR,&nbCandL);
+
+        /* If there is no candidate remaining */
         if(nbCandL == 0 || nbCandR == 0){
             nbCandL = 0;
             nbCandR = 0;
             break;
         }
     }
-    //printf("NB_ITER %d\n",k);
 
+    /* The first nbCandR values of t_candidateR correspond to the possible candidates 
+     * for the receptor
+     */
     for(i=0;i<nbCandR;i++){
         t_candidateR[i]->isCandidate = 0;
     }
+    /* The first nbCandL values of t_candidateL correspond to the possible candidates 
+     * for the ligand
+     */
     for(i=0;i<nbCandL;i++){
         t_candidateL[i]->isCandidate = 0;
     }
@@ -156,12 +181,17 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
      * for each receptor resdiue candidate, we look at each of its atoms
      * for each ligand residue candidate, we llok at each of its atoms
      */
+
+    /* For each candidate residue of the receptor */
     for(i=0;i<nbCandR;i++){
         for(j=0;j<t_candidateR[i]->nbAtom;j++){
             x1 = t_candidateR[i]->x[j];
             y1 = t_candidateR[i]->y[j];
             z1 = t_candidateR[i]->z[j];
+
+            /* For each candidate residue of the ligand */
             for(k=0;k<nbCandL;k++){
+                /* If both the residues we are looking at are already selected */
                 if(t_candidateL[k]->isCandidate == 1 && t_candidateR[i]->isCandidate == 1){
                     continue;
                 }
@@ -181,17 +211,19 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
         }
     }
 
+    /* Output the ligand's interface */
     for(i=0;i<nbCandL;i++){
         if(t_candidateL[i]->isCandidate == 1){
             fprintf(outputFile_lig,"'%s' '%c' ",strtok(t_candidateL[i]->idRes," "),t_candidateL[i]->chain);
             fflush(outputFile_lig);
-        }else{
-        
         }
     }
+    /* Output the receptor's interface */
     for(i=0;i<nbCandR;i++){
-        if(t_candidateR[i]->isCandidate == 1)
+        if(t_candidateR[i]->isCandidate == 1){
             fprintf(outputFile_rec,"'%s' '%c' ",strtok(t_candidateR[i]->idRes," "),t_candidateR[i]->chain);
+            fflush(outputFile_lig);
+        }
     }
 
 
@@ -202,29 +234,15 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL){
 
 struct docking_results* getDataForComplex_HCMD2(){
     char* buf = NULL;
-    char buf2[5];
     FILE* condFile_stream = NULL;
     int numberOfConf = 0, rc;
+    float buf3 = 0;
     size_t len = 2000;
     struct docking_results* dock_res = NULL;
 
-
-    strncpy(buf2,receptor,4);
-    buf2[4] = '\0';
-    //sprintf(totalFilePath,"%s/%s/Cond.%s.%s.UB.global.dat",dockingFile,buf2,receptor,ligand);
-
-    if(access(dockingFile,F_OK) == -1){
-        /* File does not exists, 4 file in this case
-         * Return NULL, which should be treated accordingly
-         */
-        fprintf(stderr,"No such file : %s\n",dockingFile);
-        exit(EXIT_FAILURE);
-    }
-
     condFile_stream = fopen(dockingFile,"r");
 
     /* Get the number of conformations */
-    float buf3 = 0;
     rc = fscanf(condFile_stream,"%f",&buf3);
     rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
     if(buf3 > numberOfConf)
@@ -238,7 +256,6 @@ struct docking_results* getDataForComplex_HCMD2(){
     }
     numberOfConf++;
     rewind(condFile_stream);
-    rc = getline(&buf,&len,condFile_stream); /* Get rid of the first line */
 
     dock_res = allocate_dockingResults(numberOfConf);
 

src/allocate.c

@@ -10,6 +10,7 @@ void freePDB(struct pdb_values* pdb){
 }
 
 void freeDock(struct docking_results* dock_res){
+    /* Free the docking interface */
     free(dock_res->listEner);
     free(dock_res->distCenters);
     free(dock_res->theta);
@@ -22,6 +23,9 @@ void freeDock(struct docking_results* dock_res){
 }
 
 struct docking_results* allocate_dockingResults(int nbConf){
+    /* This function will return a pointer on structure docking_results
+     * of size nbconf for each of the subsequent pointers
+     */
     int i = 0;
     struct docking_results* dock_res = NULL;
 
@@ -92,6 +96,7 @@ struct docking_results* allocate_dockingResults(int nbConf){
 }
 
 void reset_dockingResults(struct docking_results* dock_res){
+
     int i = 0;
 
     for(i=0;i<dock_res->nbConf;i++){
@@ -102,12 +107,14 @@ void reset_dockingResults(struct docking_results* dock_res){
         dock_res->alpha[i]       = 0;
         dock_res->beta[i]        = 0;
         dock_res->gamma[i]       = 0;
-        
-    }   
+
+    }
     dock_res->nbConf = 0;
 }
 
 struct pdb_values* allocate_pdb(int nbRes){
+    /* Allocate a PDB structure containing nbRes residues */
+
     struct pdb_values* pdb = NULL;
 
     pdb = malloc(sizeof(struct pdb_values));
@@ -130,6 +137,9 @@ struct pdb_values* allocate_pdb(int nbRes){
 }
 
 void free_argLine(){
+    /* Free the global variables that were allocated at the
+     * start of the program
+     */
     free(outputDir);
     free(dockingFile);
     free(pdbDir);
@@ -146,4 +156,3 @@ void free_argLine(){
         fclose(verbose_file);
     }
 }
-

src/param.c

 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <unistd.h>
 #include "struct.h"
 
 
@@ -25,6 +26,12 @@ void print_usage() {
 }
 
 struct argLine* parseLineOfArgument(int argc, char** argv){
+    /* This function aims at parsing the argument line given
+     * to the program.
+     * All the values that are set here are global variable, and
+     * therefore will not change throughout the execution of the program
+     * The declaration of these variables can be found in struct.h
+     */
     int i = 0, print_help = 0;
     char buf[500];
     struct argLine* paramVals = NULL;
@@ -189,5 +196,14 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
         outputFile_lig = fopen(buf,"w");
     }
 
+
+    if(access(dockingFile,F_OK) == -1){
+        /* File does not exists, 4 file in this case
+         * Return NULL, which should be treated accordingly
+         */
+        fprintf(stderr,"No such file : %s\n",dockingFile);
+        exit(EXIT_FAILURE);
+    }
+
     return paramVals;
 }

src/rotation.c

@@ -7,12 +7,27 @@
 
 
 void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, float R, float theta, float phi){
+    /* Sets the values of x, y and z to the carthesian equivalent of the spherical
+     * coordinate given by R, theta, phi with a center x0, y0 and z0
+     */
     *x = x0 + R * sin(theta) * cos(phi);
     *y = y0 + R * sin(theta) * sin(phi);
     *z = z0 + R * cos(theta);
 }
 
 struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR, struct pdb_values* newPDB, float R, float theta, float phi, float alpha,float beta,float gamma, struct docking_results* dock_res, int nConf){
+    /* This function rotates the PDB pdb given in argument and sets the new coordinates in newPDB. 
+     * The final center of the repair is the geometric center of the PDB pdbR. 
+     * The angles alpha, beta, gamma in arguments correspond to the rotation angles that have to be applied to
+     * each residue. 
+     * In the struct dock_res, we can find the angles alpha0, beta0 and gamma0. alpha0 corresponds to the angles alpha
+     * and beta for the first CA (from the axis x and y) + the rotation angles dalpha and dbeta respectively
+     * gamma0 corresponds to the angle gamma from the fifth CA of the protein after being rotated of alpha (NOT dalpha)
+     * then by beta (NOT dbeta).
+     * c_a0 can be seen as the unit vector of the axis X rotated by alpha0
+     *
+     * The rotation matrices are available at : https://en.wikipedia.org/wiki/Rotation_matrix in "Nested Dimensions"
+     */
     
     float xi = 0, yi = 0, zi = 0;
     float x1i = 0, y1i = 0, z1i = 0;
@@ -27,6 +42,8 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
     float c_a0=cos(dock_res->alpha[nConf]),s_a0=sin(dock_res->alpha[nConf]);
     float c_b0=cos(dock_res->beta[nConf]), s_b0=sin(dock_res->beta[nConf]);
     
+    /* These values correspond to the final rotation axis obtained after rotating
+     * the Y axis by an angle of alpha0 then beta0*/
     float rx = 0, ry = 0, rz = 0;
 
     if(newPDB == NULL){
@@ -42,52 +59,57 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
     int i = 0, j = 0;
     for(i=0;i<pdb->nbRes;i++){
         for(j=0;j<pdb->residues[i].nbAtom;j++){
-           // printf("a/b/g %f %f %f\n",alpha,beta,gamma);
-           // printf("MAXa/b/g %f %f %f\n",dock_res->alpha[nConf],dock_res->beta[nConf],dock_res->gamma[nConf]);
+
             xi = pdb->residues[i].x[j] - pdb->centerX;
             yi = pdb->residues[i].y[j] - pdb->centerY;
             zi = pdb->residues[i].z[j] - pdb->centerZ;
 
-           // printf("xi/yi/zi %f %f %f\n",xi,yi,zi);
 
-            // alpha rotation
+            // alpha rotation of the residue
             x1i=c_a*xi-s_a*yi;
             y1i=c_a*yi+s_a*xi;
             z1i=zi;
-           // printf("x1i/y1i/z1i %f %f %f\n",x1i,y1i,z1i);
 
-            // beta rotation
+            /* beta rotation of the residue */
             x2i=(c_a0*c_a0+(1.0-c_a0*c_a0)*c_b)*x1i+c_a0*s_a0*(1.0-c_b)*y1i+s_a0*s_b*z1i;
             y2i=c_a0*s_a0*(1.0-c_b)*x1i+(s_a0*s_a0+(1.0-s_a0*s_a0)*c_b)*y1i-c_a0*s_b*z1i;
             z2i=-s_a0*s_b*x1i+c_a0*s_b*y1i+c_b*z1i;
-           // printf("x2i/y2i/z2i %f %f %f\n",x2i,y2i,z2i);
             
+            /* We compute the value for the last rotation axis
+             * from the angles alpha0, beta0 and gamma0 
+             */
             rx = -s_a0*c_b0;
             ry = c_a0*c_b0;
             rz = s_b0;
 
-            // gamma rotation
+            // gamma rotation of the residue
             x3i=(rx*rx+(1.0-rx*rx)*c_g)*x2i+(rx*ry*(1.0-c_g)-rz*s_g)*y2i+(rx*rz*(1.0-c_g)+ry*s_g)*z2i;
             y3i=(rx*ry*(1.0-c_g)+rz*s_g)*x2i+(ry*ry+(1.0-ry*ry)*c_g)*y2i+(ry*rz*(1.0-c_g)-rx*s_g)*z2i;
             z3i=(rx*rz*(1.0-c_g)-ry*s_g)*x2i+(ry*rz*(1.0-c_g)+rx*s_g)*y2i+(rz*rz+(1.0-rz*rz)*c_g)*z2i;
 
-           // printf("x3i/y3i/z3i %f %f %f\n",x3i,y3i,z3i);
+
+            /* Now the residue has been rotated, and we need to translate it 
+             * to its final destination
+             */
             newPDB->residues[i].x[j] = x3i + newX;
             newPDB->residues[i].y[j] = y3i + newY;
             newPDB->residues[i].z[j] = z3i + newZ;
-           // printf("x/y/z %f %f %f\n",pdb->residues[i].x[j],pdb->residues[i].y[j],pdb->residues[i].z[j]);
-           newPDB->xCA1 = (newPDB->xCA1 - pdb->centerX) + newX;
-           newPDB->yCA1 = (newPDB->yCA1 - pdb->centerY) + newY;
-           newPDB->zCA1 = (newPDB->zCA1 - pdb->centerZ) + newZ;
-
-           newPDB->xCA5 = (newPDB->xCA5 - pdb->centerX) + newX;
-           newPDB->yCA5 = (newPDB->yCA5 - pdb->centerY) + newY;
-           newPDB->zCA5 = (newPDB->zCA5 - pdb->centerZ) + newZ;
-
-           newPDB->centerX = newX;
-           newPDB->centerY = newY;
-           newPDB->centerZ = newZ;
-           // printf("x/y/z %f %f %f\n",newPDB->residues[i].x[j],newPDB->residues[i].y[j],newPDB->residues[i].z[j]);
+
+            /* We also translate the coordinates of the first and fifth CA, even
+             * though they are not used afterwards. This is for the sake of 
+             * having the right values at the right places
+             */
+            newPDB->xCA1 = (newPDB->xCA1 - pdb->centerX) + newX;
+            newPDB->yCA1 = (newPDB->yCA1 - pdb->centerY) + newY;
+            newPDB->zCA1 = (newPDB->zCA1 - pdb->centerZ) + newZ;
+
+            newPDB->xCA5 = (newPDB->xCA5 - pdb->centerX) + newX;
+            newPDB->yCA5 = (newPDB->yCA5 - pdb->centerY) + newY;
+            newPDB->zCA5 = (newPDB->zCA5 - pdb->centerZ) + newZ;
+
+            newPDB->centerX = newX;
+            newPDB->centerY = newY;
+            newPDB->centerZ = newZ;
         }
     }
     return newPDB;
@@ -96,6 +118,11 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
 
 
 void getAnglesFromCA1andCA5_sophie(struct pdb_values* pdb, float* alpha, float* beta, float* gamma){
+    /* The goal of this function is to find the angles alpha, beta and gamma that correspond 
+     * to the angles alpha and beta of the first CA of the protein
+     * We then have to rotate the fifth CA by the angle alpha, then beta in order to compute the gamma angle
+     */
+
     *alpha = 0;
     *beta = 0;
     *gamma = 0;
@@ -108,23 +135,28 @@ void getAnglesFromCA1andCA5_sophie(struct pdb_values* pdb, float* alpha, float*
     float c_b = 0, s_b = 0;
     float c_g = 0, s_g = 0;
 
-    // CA 1
+    // CA 1 coordinates
     x1 = pdb->xCA1 - pdb->centerX;
     y1 = pdb->yCA1 - pdb->centerY;
     z1 = pdb->zCA1 - pdb->centerZ;
 
-    // CA 5
+    // CA 5 coordinates
     x2 = pdb->xCA5 - pdb->centerX;
     y2 = pdb->yCA5 - pdb->centerY;
     z2 = pdb->zCA5 - pdb->centerZ;
     
     // Lenght of the projection on the plane X,Y
     p1 = sqrt(x1*x1+y1*y1);
+    
+    // Distance between the CA 1 and the geometric center of the ligand
     r1 = sqrt(x1*x1+y1*y1+z1*z1);
 
+    // cos and sin of the alpha angle
+    // The sin has to be computed negatively (for some reason)
     c_a = y1/p1;
     s_a = -x1/p1;
 
+    // cos and sin of the beta angle
     c_b = p1/r1;
     s_b = z1/r1;
 
@@ -137,15 +169,17 @@ void getAnglesFromCA1andCA5_sophie(struct pdb_values* pdb, float* alpha, float*
         *beta = -(*beta);
 
     
-    /* Clockwise rotation matrix for alpha */
+    /* Rotation matrix for alpha */
     x1i = x2*c_a - y2*s_a;
     y1i = y2*c_a + x2*s_a;
     z1i = z2;
 
+    /* Rotation matrix for beta */
     x2i = x1i;
     y2i = y1i*c_b + z1i*s_b;
     z2i = z1i*c_b - y1i*s_b;
 
+    /* From this point we are able to compute the gamma angle */
     p2 = sqrt(x2i*x2i+z2i*z2i);
     c_g = z2i/p2;
     s_g = x2i/p2;