ZDOCK : advancements

src/INTBuilder.c

@@ -325,6 +325,21 @@ int main(int argc, char** argv){
                                 dock_res->trans_Y[target_conf],dock_res->trans_Z[target_conf],
                                 dock_res->alpha[target_conf],dock_res->beta[target_conf],
                                 dock_res->gamma[target_conf]);
+        }else if(ZDOCK){
+            printf("R %f %f %f\n",pdbR->centerX,pdbR->centerY,pdbR->centerZ);
+            printf("L %f %f %f\n",pdbL->centerX,pdbL->centerY,pdbL->centerZ);
+            init_ZDOCK(&pdbR,&pdbL);
+            printf("R %f %f %f\n",pdbR->centerX,pdbR->centerY,pdbR->centerZ);
+            printf("L %f %f %f\n",pdbL->centerX,pdbL->centerY,pdbL->centerZ);
+        if(constructPDB){
+            writePDB(pdbL, ligand,target_conf);
+            writePDB(pdbR, receptor,target_conf);
+        }
+
+            newPDB = rotate_ZDOCK(pdbL,newPDB,dock_res->trans_X[target_conf],
+                                dock_res->trans_Y[target_conf],dock_res->trans_Z[target_conf],
+                                dock_res->alpha[target_conf],dock_res->beta[target_conf],
+                                dock_res->gamma[target_conf]);
         }else{
             newPDB = rotate_global(pdbL,pdbR,newPDB,dock_res->distCenters[target_conf],
                                 dock_res->theta[target_conf],dock_res->phi[target_conf],
@@ -371,7 +386,10 @@ int main(int argc, char** argv){
                                     dock_res->alpha[i],dock_res->beta[i],
                                     dock_res->gamma[i]);
             }else if(ZDOCK){
+                printf("%f %f %f\n",pdbR->centerX,pdbR->centerY,pdbR->centerZ);
                 init_ZDOCK(&pdbR,&pdbL);
+                printf("%f %f %f\n",pdbR->centerX,pdbR->centerY,pdbR->centerZ);
+                exit(1);
 
                 newPDB = rotate_ZDOCK(pdbL,newPDB,dock_res->trans_X[i],
                                     dock_res->trans_Y[i],dock_res->trans_Z[i],

src/allocate.c

@@ -85,7 +85,7 @@ struct docking_results* allocate_dockingResults(int nbConf){
     }
 
     /* Only HEX uses translation values */
-    if(HEX){
+    if(HEX || ZDOCK){
         dock_res->trans_X = malloc((2*nbConf)*sizeof(int));
         if(dock_res->trans_X == NULL){
             perror("malloc");
@@ -128,7 +128,7 @@ struct docking_results* allocate_dockingResults(int nbConf){
         dock_res->beta[i]        = 0;
         dock_res->gamma[i]       = 0;
         dock_res->t_conf[i]      = 0;
-        if(HEX){
+        if(HEX || ZDOCK){
             dock_res->trans_X[i]    = 0;
             dock_res->trans_Y[i]    = 0;
             dock_res->trans_Z[i]    = 0;

src/fileIO.c

@@ -7,6 +7,7 @@
 
 struct docking_results* getDataForComplex(){
     char* buf = NULL;
+    char buf2[100] = "";
     FILE* condFile_stream = NULL;
     int numberOfConf = 0, rc = 0;
     int idConf = 0, buf4 = 0, i = 0;
@@ -18,16 +19,16 @@ struct docking_results* getDataForComplex(){
 
     if(HEX){
         for(i=0;i<NB_HEADER_LINE_HEX;i++){
-            getline(&buf,&len,condFile_stream);
+            rc = getline(&buf,&len,condFile_stream);
         }
     }
 
     if(ZDOCK){
-        getline(&buf,&len,condFile_stream);
-        fscanf(condFile_stream,"%f %f %f",&ZDOCK_REC_A_ROT,&ZDOCK_REC_B_ROT,&ZDOCK_REC_G_ROT);
-        fscanf(condFile_stream,"%f %f %f",&ZDOCK_LIG_A_ROT,&ZDOCK_LIG_B_ROT,&ZDOCK_LIG_G_ROT);
-        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_REC_X_TRANS,&ZDOCK_REC_Y_TRANS,&ZDOCK_REC_Z_TRANS);
-        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_LIG_X_TRANS,&ZDOCK_LIG_Y_TRANS,&ZDOCK_LIG_Z_TRANS);
+        rc = fscanf(condFile_stream,"%d %f %d",&ZDOCK_GRID_SIZE,&ZDOCK_GRID_UNIT,&ZDOCK_PROT_FIXED);
+        rc = fscanf(condFile_stream,"%f %f %f",&ZDOCK_REC_A_ROT,&ZDOCK_REC_B_ROT,&ZDOCK_REC_G_ROT);
+        rc = fscanf(condFile_stream,"%f %f %f",&ZDOCK_LIG_A_ROT,&ZDOCK_LIG_B_ROT,&ZDOCK_LIG_G_ROT);
+        rc = fscanf(condFile_stream,"%s %f %f %f",buf2,&ZDOCK_REC_X_TRANS,&ZDOCK_REC_Y_TRANS,&ZDOCK_REC_Z_TRANS);
+        rc = fscanf(condFile_stream,"%s %f %f %f",buf2,&ZDOCK_LIG_X_TRANS,&ZDOCK_LIG_Y_TRANS,&ZDOCK_LIG_Z_TRANS);
     }
 
     /* Get the number of conformations */
@@ -36,8 +37,11 @@ struct docking_results* getDataForComplex(){
     }else{
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
     }
-    if(buf3 > numberOfConf)
+    if(ZDOCK){
+        numberOfConf++;
+    }else if(buf3 > numberOfConf){
         numberOfConf = (int) buf3;
+    }
 
     while(!feof(condFile_stream)){
         if(HEX){
@@ -46,18 +50,27 @@ struct docking_results* getDataForComplex(){
             rc = fscanf(condFile_stream,"%f",&buf3);
         }
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
-        if(buf3 > numberOfConf)
+        if(ZDOCK){
+            if(feof(condFile_stream)){
+                break;
+            }
+            numberOfConf++;
+        }else if(buf3 > numberOfConf)
             numberOfConf = (int) buf3;
     }
-    //numberOfConf++;
+    if(ZDOCK)
+        numberOfConf--;
     rewind(condFile_stream);
 
     /* This is with the first MAXDo format */
     if(MAXDo){
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the first line */
+    }else if(ZDOCK){
+        for(i=0;i<NB_HEADER_LINE_ZDOCK;i++)
+            rc = getline(&buf,&len,condFile_stream);
     }else if(HEX){
         for(i=0;i<NB_HEADER_LINE_HEX;i++){
-            getline(&buf,&len,condFile_stream);
+            rc = getline(&buf,&len,condFile_stream);
         }
     }
 
@@ -73,7 +86,9 @@ struct docking_results* getDataForComplex(){
             rc = fscanf(condFile_stream,"%d %d",&idConf,&buf4);
         }else if(HEX){
             rc = fscanf(condFile_stream,"%d %d %f %f",&buf4,&idConf,&buf3,&buf3);
-        }else{
+        }else if(ZDOCK){
+            idConf++;
+        }else {
             /* Else, that means that we have a format where
              * values for distance between center follows directly
              * the ID of the conformation.
@@ -89,14 +104,16 @@ struct docking_results* getDataForComplex(){
             break;
         }
 
-
         if(HEX){
             rc = fscanf(condFile_stream,"%f %f %f %f %f %f",&dock_res->trans_X[idConf],&dock_res->trans_Y[idConf],&dock_res->trans_Z[idConf],&dock_res->alpha[idConf],&dock_res->beta[idConf],&dock_res->gamma[idConf]);
+        }else if(ZDOCK){
+            rc = fscanf(condFile_stream,"%f %f %f %f %f %f",&dock_res->alpha[idConf],&dock_res->beta[idConf],&dock_res->gamma[idConf],&dock_res->trans_X[idConf],&dock_res->trans_Y[idConf],&dock_res->trans_Z[idConf]);
         }else{
             rc = fscanf(condFile_stream,"%f %f %f %f %f %f",&dock_res->distCenters[idConf],&dock_res->theta[idConf],&dock_res->phi[idConf],&dock_res->alpha[idConf],&dock_res->beta[idConf],&dock_res->gamma[idConf]);
         }
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
 
+
         /* theta and phi angles are written in Radian already
          * However that is not the case for alpha, beta and gamma angles
          */
@@ -112,6 +129,11 @@ struct docking_results* getDataForComplex(){
             dock_res->beta[idConf]  = dock_res->beta[idConf]  * COEF_DEGREE_TO_RADIAN;
             dock_res->gamma[idConf] = dock_res->gamma[idConf] * COEF_DEGREE_TO_RADIAN;
         }
+
+        /* This happens when we arrive on the last line */
+        if(feof(condFile_stream)){
+            break;
+        }
     }
 
     dock_res->nbConf = numberOfConf;

src/param.c

@@ -8,10 +8,10 @@
 
 
 void print_usage() {
-    printf("Usage: \n"
+    printf("Usage:\n"
     "Mandatory\n"
     "=========\n\n"
-    "<-dockingFile>     MAXDo file that needs to be analysed\n"
+    "<-dockingFile>     docking file that needs to be analysed\n"
     "OR\n"
     "<-complex>         This option means that we look at two PDBs instead of looking\n"
     "                   at the dockingFile. Both PBDs for '-rec' and '-lig' must be in\n"
@@ -41,6 +41,8 @@ void print_usage() {
     "<-MAXDo>...........Means that the docking file is in the original MAXDo format\n"
     "OR\n"
     "<-HEX>.............Means that the docking file is in the HEX format\n"
+    "OR\n"
+    "<-ZDOCK>...........Means that the docking file is in the ZDOCK format\n"
     );
 }
 
@@ -70,12 +72,14 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
     target_conf      = -1;
     DIST_FOR_CONTACT = -1;
     ZDOCK_GRID_UNIT  = 0;
-    ZDOCK_LIG_X_ROT  = 0;
-    ZDOCK_LIG_Y_ROT  = 0;
-    ZDOCK_LIG_Z_ROT  = 0;
-    ZDOCK_REC_X_ROT  = 0;
-    ZDOCK_REC_Y_ROT  = 0;
-    ZDOCK_REC_Z_ROT  = 0;
+    ZDOCK_GRID_SIZE  = 0;
+    ZDOCK_PROT_FIXED = 0;
+    ZDOCK_LIG_A_ROT  = 0;
+    ZDOCK_LIG_B_ROT  = 0;
+    ZDOCK_LIG_G_ROT  = 0;
+    ZDOCK_REC_A_ROT  = 0;
+    ZDOCK_REC_B_ROT  = 0;
+    ZDOCK_REC_G_ROT  = 0;
     ZDOCK_LIG_X_TRANS= 0;
     ZDOCK_LIG_Y_TRANS= 0;
     ZDOCK_LIG_Z_TRANS= 0;

src/rotation.c

@@ -10,13 +10,13 @@ void init_ZDOCK(struct pdb_values** pdbR_addr, struct pdb_values** pdbL_addr){
     struct pdb_values* pdbL = *pdbL_addr;
     struct pdb_values* newPDB = NULL;
     
-    newPDB = rotate_ZDOCK(pdbR,newPDB,ZDOCK_REC_X_TRANS,ZDOCK_REC_Y_TRANS,ZDOCK_REC_Z_TRANS,ZDOCK_REC_A_ROT,ZDOCK_REC_B_ROT,ZDOCK_REC_G_ROT);
-    pdbR = newPDB;
+    newPDB = rotate_ZDOCK_init(pdbR,newPDB,ZDOCK_REC_X_TRANS,ZDOCK_REC_Y_TRANS,ZDOCK_REC_Z_TRANS,ZDOCK_REC_A_ROT,ZDOCK_REC_B_ROT,ZDOCK_REC_G_ROT);
+    (*pdbR_addr) = newPDB;
     newPDB = NULL;
     // TODO : free pdbR
 
-    newPDB = rotate_ZDOCK(pdbL,newPDB,ZDOCK_LIG_X_TRANS,ZDOCK_LIG_Y_TRANS,ZDOCK_LIG_Z_TRANS,ZDOCK_LIG_A_ROT,ZDOCK_LIG_B_ROT,ZDOCK_LIG_G_ROT);
-    pdbL = newPDB;
+    newPDB = rotate_ZDOCK_init(pdbL,newPDB,ZDOCK_LIG_X_TRANS,ZDOCK_LIG_Y_TRANS,ZDOCK_LIG_Z_TRANS,ZDOCK_LIG_A_ROT,ZDOCK_LIG_B_ROT,ZDOCK_LIG_G_ROT);
+    (*pdbL_addr) = newPDB;
     newPDB = NULL;
     // TODO : free pdbL
 }
@@ -46,6 +46,8 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
      * https://upload.wikimedia.org/wikipedia/commons/thumb/a/a1/Eulerangles.svg/300px-Eulerangles.svg.png
      */
 
+    printf("rotate : T X/Y/Z : %f _%f _%f\nR A/B/G : %f %f %f\n",trans_X,trans_Y,trans_Z,alpha,beta,gamma);
+
     float xi = 0, yi = 0, zi = 0;
     float x1i = 0, y1i = 0, z1i = 0;
     float x2i = 0, y2i = 0, z2i = 0;
@@ -61,14 +63,10 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
     float yN = s_a;
     float zN = 0;
 
-    /* The coordinates of the unit vector for the axis Z can be computed as such :
-     *  xZ : c_b
-     *  yZ : c_a * s_b
-     *  zZ : s_a * s_b
-     */
-    float xZ = c_b;
-    float yZ = c_a * s_b;
-    float zZ = s_a * s_b;
+    /* The coordinates of the unit vector for the axis Z can be computed as such : */
+    float xZ = yN*s_b;
+    float yZ = -xN*s_b;
+    float zZ = c_b;
 
     if(newPDB == NULL){
         newPDB = allocate_pdb(pdb->nbRes);
@@ -78,9 +76,9 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
 
 
     /* New center of the protein */
-    float newX = pdb->centerX + trans_X;
-    float newY = pdb->centerY + trans_Y;
-    float newZ = pdb->centerZ + trans_Z;
+    float newX = pdb->centerX + (trans_X-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
+    float newY = pdb->centerY + (trans_Y-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
+    float newZ = pdb->centerZ + (trans_Z-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
 
 
     /* This is now the part where we actually compute the rotations, given
@@ -114,13 +112,114 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
             /* 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;
+            newPDB->residues[i].x[j] = x3i; // + newX;
+            newPDB->residues[i].y[j] = y3i; // + newY;
+            newPDB->residues[i].z[j] = z3i; // + newZ;
+
+        }
+    }
+
+    newPDB->centerX = newX;
+    newPDB->centerY = newY;
+    newPDB->centerZ = newZ;
+    printf("NEW : %f %f %f\n",newPDB->centerX,newPDB->centerY,newPDB->centerZ);
+    printf("trans : %f %f %f\n",trans_X,trans_Y,trans_Z);
+
+    return newPDB;
+}
+
+struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha, float beta, float gamma){
+    /* This function aims at rotating the pdb parameter and gets its new coordinates
+     * in the newPDB parameter.
+     * Let the N axis the rotation of the x axis about the z axis by an alpha angle.
+     * Let the Z axis the rotation of the z axis about the N axis by a beta angle.
+     *
+     * The xyz system is centered on the center of mass of the pdb parameter.
+     * We then rotate the pdb by an alpha angle about the axis z, then by a beta
+     * angle about the axis N and finally by a gamma angle about the axis Z.
+     * The following image available on wikipedia can help to understand the scheme :
+     * https://upload.wikimedia.org/wikipedia/commons/thumb/a/a1/Eulerangles.svg/300px-Eulerangles.svg.png
+     */
+
+    float xi = 0, yi = 0, zi = 0;
+    float x1i = 0, y1i = 0, z1i = 0;
+    float x2i = 0, y2i = 0, z2i = 0;
+    float x3i = 0, y3i = 0, z3i = 0;
 
+    float c_a=cos(alpha),s_a=sin(alpha);
+    float c_b=cos(beta),s_b=sin(beta);
+    float c_g=cos(gamma),s_g=sin(gamma);
+
+
+    /* Coordinates of the N axis */
+    float xN = c_a;
+    float yN = s_a;
+    float zN = 0;
+    printf("N AXIS : %f %f %f\n",xN,yN,zN);
+
+    /* The coordinates of the unit vector for the axis Z can be computed as such :
+     */
+    float xZ = yN*s_b;
+    float yZ = -xN*s_b;
+    float zZ = c_b;
+    printf("Z AXIS : %f %f %f\n",xZ,yZ,zZ);
+
+    if(newPDB == NULL){
+        newPDB = allocate_pdb(pdb->nbRes);
+        newPDB->nbAtom = pdb->nbAtom;
+    }
+    memcpy(newPDB->residues,pdb->residues,pdb->nbRes*sizeof(struct residue));
+
+
+    /* New center of the protein */
+    float newX = 0;
+    float newY = 0;
+    float newZ = 0;
+
+
+    /* This is now the part where we actually compute the rotations, given
+     * the coordinates for the rotation axis N and Z
+     */
+    int i = 0, j = 0;
+    for(i=0;i<pdb->nbRes;i++){
+        for(j=0;j<pdb->residues[i].nbAtom;j++){
+
+            xi = pdb->residues[i].x[j];
+            yi = pdb->residues[i].y[j];
+            zi = pdb->residues[i].z[j];
+
+
+            /* alpha rotation of the residue */
+            x1i=c_a*xi-s_a*yi;
+            y1i=c_a*yi+s_a*xi;
+            z1i=zi;
+
+            /* beta rotation of the residue */
+            x2i = (xN*xN+(1.0-xN*xN)*c_b)*x1i + xN*yN*(1.0-c_b)*y1i         + yN*s_b*z1i;
+            y2i = xN*yN*(1.0-c_b)*x1i         + (yN*yN+(1.0-yN*yN)*c_b)*y1i - xN*s_b*z1i;
+            z2i = -yN*s_b*x1i                 + xN*s_b*y1i                  + c_b*z1i;
+
+            /* gamma rotation of the residue */
+            x3i = (xZ*xZ+(1.0-xZ*xZ)*c_g) *x2i + (xZ*yZ*(1.0-c_g)-zZ*s_g)*y2i + (xZ*zZ*(1.0-c_g)+yZ*s_g)*z2i;
+            y3i = (xZ*yZ*(1.0-c_g)+zZ*s_g)*x2i + (yZ*yZ+(1.0-yZ*yZ)*c_g) *y2i + (yZ*zZ*(1.0-c_g)-xZ*s_g)*z2i;
+            z3i = (xZ*zZ*(1.0-c_g)-yZ*s_g)*x2i + (yZ*zZ*(1.0-c_g)+xZ*s_g)*y2i + (zZ*zZ+(1.0-zZ*zZ)*c_g) *z2i;
+
+
+            /* Now the residue has been rotated, and we need to translate it
+             * to its final destination
+             */
+            newPDB->residues[i].x[j] = x3i - trans_X;
+            newPDB->residues[i].y[j] = y3i - trans_Y;
+            newPDB->residues[i].z[j] = z3i - trans_Z;
+            newX += newPDB->residues[i].x[j];
+            newY += newPDB->residues[i].y[j];
+            newZ += newPDB->residues[i].z[j];
         }
     }
 
+    newX /= newPDB->nbAtom;
+    newY /= newPDB->nbAtom;
+    newZ /= newPDB->nbAtom;
     newPDB->centerX = newX;
     newPDB->centerY = newY;
     newPDB->centerZ = newZ;
@@ -215,13 +314,10 @@ struct pdb_values* rotate_global(struct pdb_values* pdb, struct pdb_values* pdbR
     float zN = 0;
 
     /* The coordinates of the unit vector for the axis Z can be computed as such :
-     *  xZ : c_b
-     *  yZ : c_a * s_b
-     *  zZ : s_a * s_b
      */
-    float xZ = c_b;
-    float yZ = c_a * s_b;
-    float zZ = s_a * s_b;
+    float xZ = yN*s_b;
+    float yZ = -xN*s_b;
+    float zZ = c_b;
 
     if(newPDB == NULL){
         newPDB = allocate_pdb(pdb->nbRes);

src/rotation.h

@@ -13,6 +13,8 @@ struct pdb_values* rotate_global(struct pdb_values* pdb, struct pdb_values* pdbR
 
 struct pdb_values* rotate_HEX(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha,float beta,float gamma);
 
+struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha, float beta, float gamma);
+
 struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha, float beta, float gamma);
 
 void init_ZDOCK(struct pdb_values** pdbR_addr, struct pdb_values** pdbL_addr);

src/struct.h

@@ -16,7 +16,7 @@
 #define DIST_FOR_CLASH 2           /* Distance in Angstrom under which we consider
                                     * the contact is a clash */
 #define NB_HEADER_LINE_HEX 12        /* Number of header lines at the top of a hex docking file */
-
+#define NB_HEADER_LINE_ZDOCK 5        /* Number of header lines at the top of a ZDOCK docking file */
 /* This is the format of a line in a PDB file */
 #define FORMAT_LINE_PDB "ATOM  %5d  %3s %3s %c%5s   %8.3f%8.3f%8.3f\n"
 
@@ -46,6 +46,8 @@ float   ZDOCK_REC_X_TRANS;
 float   ZDOCK_REC_Y_TRANS;
 float   ZDOCK_REC_Z_TRANS;
 float   ZDOCK_GRID_UNIT;
+int     ZDOCK_GRID_SIZE;
+int     ZDOCK_PROT_FIXED;
 FILE*   verbose_file;       /* If verbose mode, the output file */
 char*   pdbDir;             /* Directory containing the PDBs */
 char*   receptor;           /* Name of the receptor */