src/contrib/do_multiprot.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "string2.h"
  46 #include "strdb.h"
  47 #include "macros.h"
  48 #include "smalloc.h"
  49 #include "mshift.h"
  50 #include "statutil.h"
  51 #include "copyrite.h"
  52 #include "pdbio.h"
  53 #include "gmx_fatal.h"
  54 #include "xvgr.h"
  55 #include "matio.h"
  56 #include "index.h"
  57 #include "gstat.h"
  58 #include "tpxio.h"
  59 #include "viewit.h"
  60 #include "gbutil.h"
  61 #include "vec.h"
  62 #include "confio.h"
  63 #include "gmxfio.h"
  64
  65 typedef struct {
  66     int resnr;
  67     int count;
  68 } t_countres;
  69
  70 static void process_multiprot_output(const char *fn, real *rmsd, int *nres, rvec rotangles,
  71                                      rvec translation, bool bCountres, t_countres *countres)
  72 {
  73     FILE       *mpoutput;
  74     char       line[256];
  75     char       *string;
  76     int        i=0,j=0,res;
  77
  78     (*rmsd)=-1;
  79     (*nres)=0;
  80     mpoutput=ffopen (fn,"r");
  81
  82     if (bCountres) {
  83         do {
  84             fgets(line, 256, mpoutput);
  85         } while (strstr(line,"Match List") == NULL);
  86         fgets(line, 256, mpoutput);
  87         do {
  88             string = strtok (line,".");
  89             while (i<2 && string != NULL) {
  90                 string = strtok (NULL,". ");
  91                 i++;
  92             }
  93             i=0;
  94             res=atoi(string);
  95             if (res > 0) {
  96                 while (countres[j].resnr!=res)
  97                     j++;
  98                 countres[j].count++;
  99             }
 100             fgets(line, 256, mpoutput);
 101         } while (strstr(line,"End of Match List") == NULL);
 102         rewind(mpoutput);
 103     }
 104     do {
 105         fgets(line, 256, mpoutput);
 106     } while (strstr(line,"Trans : ") == NULL);
 107
 108     string = strtok (line," :");
 109     string = strtok (NULL," ");
 110     while (i<3 && string != NULL) {
 111         string = strtok (NULL," ");
 112         rotangles[i]=atof(string);
 113         i++;
 114     }
 115     i=0;
 116     while (i<3 && string != NULL) {
 117         string = strtok (NULL," ");
 118         translation[i]=atof(string)/10;
 119         i++;
 120     }
 121     if (i!=3) {
 122         gmx_warning("Not enough values for rotation and translation vectors in the output of multiprot");
 123     }
 124
 125     rotangles[YY]=rotangles[YY]*(-1);
 126
 127     while ((*rmsd) <0) {
 128         fgets(line, 256, mpoutput);
 129         if (strstr(line,"RMSD : ") != NULL) {
 130             string = strtok (line,":");
 131             string = strtok (NULL,":");
 132             (*rmsd)=atof(string)/10;
 133         }
 134     }
 135     while (!(*nres)) {
 136         fgets(line,256, mpoutput);
 137         if (strstr(line,"Match List") != NULL) {
 138             string = strtok (line,":");
 139             string = strtok (NULL,":");
 140             (*nres) = atoi(string);
 141         }
 142     }
 143     ffclose(mpoutput);
 144 }
 145
 146 int main(int argc,char *argv[])
 147 {
 148     const char *desc[] = {
 149         "[TT]do_multiprot[tt] ",
 150         "reads a trajectory file and aligns it to a reference structure  ",
 151         "each time frame",
 152         "calling the multiprot program. This allows you to use a reference",
 153         "structure whose sequence is different than that of the protein in the ",
 154         "trajectory, since the alignment is based on the geometry, not sequence.",
 155         "The output of [TT]do_multiprot[tt] includes the rmsd and the number of residues",
 156         "on which it was calculated.",
 157         "[PAR]",
 158         "An aligned trajectory file is generated with the [TT]-ox[tt] option.[PAR]",
 159         "With the [TT]-cr[tt] option, the number of hits in the alignment is given",
 160         "per residue. This number can be between 0 and the number of frames, and",
 161         "indicates the structural conservation of this residue.[PAR]",
 162         "If you do not have the [TT]multiprot[tt] program, get it. [TT]do_multiprot[tt] assumes",
 163         "that the [TT]multiprot[tt] executable is [TT]/usr/local/bin/multiprot[tt]. If this is ",
 164         "not the case, then you should set an environment variable [BB]MULTIPROT[bb]",
 165         "pointing to the [TT]multiprot[tt] executable, e.g.: [PAR]",
 166         "[TT]setenv MULTIPROT /usr/MultiProtInstall/multiprot.Linux[tt][PAR]",
 167         "Note that at the current implementation only binary alignment (your",
 168         "molecule to a reference) is supported. In addition, note that the ",
 169         "by default [TT]multiprot[tt] aligns the two proteins on their C-alpha carbons.",
 170         "and that this depends on the [TT]multiprot[tt] parameters which are not dealt ",
 171         "with here. Thus, the C-alpha carbons is expected to give the same "
 172         "results as choosing the whole protein and will be slightly faster.[PAR]",
 173         "For information about [TT]multiprot[tt], see:",
 174         "http://bioinfo3d.cs.tau.ac.il/MultiProt/.[PAR]"
 175     };
 176     static bool bVerbose;
 177     t_pargs pa[] = {
 178         { "-v",  FALSE, etBOOL, {&bVerbose},
 179           "HIDDENGenerate miles of useless information" }
 180     };
 181
 182     const char *bugs[] = {
 183         "The program is very slow, since multiprot is run externally"
 184     };
 185
 186     t_trxstatus *status;
 187     t_trxstatus *trxout=NULL;
 188     FILE        *tapein,*fo,*frc,*tmpf,*out=NULL,*fres=NULL;
 189     const char  *fnRef;
 190     const char  *fn="2_sol.res";
 191     t_topology  top;
 192     int         ePBC;
 193     t_atoms     *atoms,ratoms,useatoms;
 194     t_trxframe  fr;
 195     t_pdbinfo   p;
 196     int         nres,nres2,nr0;
 197     real        t;
 198     int         i,j,natoms,nratoms,nframe=0,model_nr=-1;
 199     int         cur_res,prev_res;
 200     int         nout;
 201     t_countres  *countres=NULL;
 202     matrix      box,rbox;
 203     int         gnx;
 204     char        *grpnm,*ss_str;
 205     atom_id     *index;
 206     rvec        *xp,*x,*xr;
 207     char        pdbfile[32],refpdb[256],title[256],rtitle[256],filemode[5];
 208     char        out_title[256];
 209     char        multiprot[256],*mptr;
 210     int         ftp;
 211     int         outftp=-1;
 212     real        rmsd;
 213     bool        bTrjout,bCountres;
 214     const char  *TrjoutFile=NULL;
 215     output_env_t oenv;
 216     static rvec translation={0,0,0},rotangles={0,0,0};
 217     gmx_rmpbc_t gpbc=NULL;
 218
 219     t_filenm   fnm[] = {
 220         { efTRX, "-f",   NULL,      ffREAD },
 221         { efTPS, NULL,   NULL,      ffREAD },
 222         { efNDX, NULL,   NULL,      ffOPTRD },
 223         { efSTX, "-r",   NULL     , ffREAD },
 224         { efXVG, "-o",  "rmss",     ffWRITE },
 225         { efXVG, "-rc", "rescount", ffWRITE},
 226         { efXVG, "-cr", "countres", ffOPTWR},
 227         { efTRX, "-ox", "aligned",  ffOPTWR }
 228     };
 229 #define NFILE asize(fnm)
 230
 231     CopyRight(stderr,argv[0]);
 232     parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
 233                       NFILE,fnm, asize(pa),pa, asize(desc),desc,
 234                       asize(bugs),bugs,&oenv
 235         );
 236     fnRef=opt2fn("-r",NFILE,fnm);
 237     bTrjout = opt2bSet("-ox",NFILE,fnm);
 238     bCountres=  opt2bSet("-cr",NFILE,fnm);
 239
 240     if (bTrjout) {
 241         TrjoutFile = opt2fn_null("-ox",NFILE,fnm);
 242     }
 243
 244     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
 245     gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
 246     atoms=&(top.atoms);
 247
 248     ftp=fn2ftp(fnRef);
 249
 250     get_stx_coordnum(fnRef,&nratoms);
 251     init_t_atoms(&ratoms,nratoms,TRUE);
 252     snew(xr,nratoms);
 253     read_stx_conf(fnRef,rtitle,&ratoms,xr,NULL,&ePBC,rbox);
 254
 255     if (bVerbose) {
 256         fprintf(stderr,"Read %d atoms\n",atoms->nr);
 257         fprintf(stderr,"Read %d reference atoms\n",ratoms.nr);
 258     }
 259     if (bCountres) {
 260         snew(countres,ratoms.nres);
 261         j=0;
 262         cur_res=0;
 263         for (i=0;i<ratoms.nr;i++) {
 264             prev_res=cur_res;
 265             cur_res=ratoms.atom[i].resind;
 266             if (cur_res != prev_res) {
 267                 countres[j].resnr=cur_res;
 268                 countres[j].count=0;
 269                 j++;
 270             }
 271         }
 272     }
 273     get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
 274     nres=0;
 275     nr0=-1;
 276     for(i=0; (i<gnx); i++) {
 277         if (atoms->atom[index[i]].resind != nr0) {
 278             nr0=atoms->atom[index[i]].resind;
 279             nres++;
 280         }
 281     }
 282     fprintf(stderr,"There are %d residues in your selected group\n",nres);
 283
 284     strcpy(pdbfile,"ddXXXXXX");
 285     gmx_tmpnam(pdbfile);
 286     if ((tmpf = fopen(pdbfile,"w")) == NULL) {
 287         sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
 288         gmx_tmpnam(pdbfile);
 289         if ((tmpf = fopen(pdbfile,"w")) == NULL) {
 290             gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);
 291         }
 292     }
 293     else {
 294         ffclose(tmpf);
 295     }
 296
 297     if (ftp != efPDB) {
 298         strcpy(refpdb,"ddXXXXXX");
 299         gmx_tmpnam(refpdb);
 300         strcat(refpdb,".pdb");
 301         write_sto_conf(refpdb,rtitle,&ratoms,xr,NULL,ePBC,rbox);
 302     }
 303     else {
 304         strcpy(refpdb,fnRef);
 305     }
 306
 307     if ((mptr=getenv("MULTIPROT")) == NULL) {
 308         mptr="/usr/local/bin/multiprot";
 309     }
 310     if (!gmx_fexist(mptr)) {
 311         gmx_fatal(FARGS,"MULTIPROT executable (%s) does not exist (use setenv MULTIPROT)",
 312                   mptr);
 313     }
 314     sprintf (multiprot,"%s %s %s > /dev/null %s",
 315              mptr, refpdb, pdbfile, "2> /dev/null");
 316
 317     if (bVerbose)
 318         fprintf(stderr,"multiprot cmd='%s'\n",multiprot);
 319
 320     if (!read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_READ_X))
 321         gmx_fatal(FARGS,"Could not read a frame from %s",ftp2fn(efTRX,NFILE,fnm));
 322     natoms = fr.natoms;
 323
 324     if (bTrjout) {
 325         nout=natoms;
 326         /* open file now */
 327         outftp=fn2ftp(TrjoutFile);
 328         if (bVerbose)
 329             fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(outftp));
 330         strcpy(filemode,"w");
 331         switch (outftp) {
 332             case efXTC:
 333             case efG87:
 334             case efTRR:
 335             case efTRJ:
 336                 out=NULL;
 337                 trxout = open_trx(TrjoutFile,filemode);
 338                 break;
 339             case efGRO:
 340             case efG96:
 341             case efPDB:
 342                 /* Make atoms struct for output in GRO or PDB files */
 343                 /* get memory for stuff to go in pdb file */
 344                 init_t_atoms(&useatoms,nout,FALSE);
 345                 sfree(useatoms.resinfo);
 346                 useatoms.resinfo=atoms->resinfo;
 347                 for(i=0;(i<nout);i++) {
 348                     useatoms.atomname[i]=atoms->atomname[i];
 349                     useatoms.atom[i]=atoms->atom[i];
 350                     useatoms.nres=max(useatoms.nres,useatoms.atom[i].resind+1);
 351                 }
 352                 useatoms.nr=nout;
 353                 out=ffopen(TrjoutFile,filemode);
 354                 break;
 355         }
 356         if (outftp == efG87)
 357             fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),"Generated by %s. #atoms=%d, a BOX is"
 358                     " stored in this file.\n",Program(),nout);
 359     }
 360
 361     if (natoms > atoms->nr) {
 362         gmx_fatal(FARGS,"\nTrajectory does not match topology!");
 363     }
 364     if (gnx > natoms) {
 365         gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
 366     }
 367
 368     fo = xvgropen(opt2fn("-o",NFILE,fnm),"RMSD","Time (ps)","RMSD (nm)",oenv);
 369     frc = xvgropen(opt2fn("-rc",NFILE,fnm),"Number of Residues in the alignment","Time (ps)","Residues",oenv);
 370
 371     do {
 372         t = output_env_conv_time(oenv,fr.time);
 373         gmx_rmpbc(gpbc,natoms,fr.box,fr.x);
 374         tapein=ffopen(pdbfile,"w");
 375         write_pdbfile_indexed(tapein,NULL,atoms,fr.x,ePBC,fr.box,' ',-1,gnx,index,NULL,TRUE);
 376         ffclose(tapein);
 377         system(multiprot);
 378         remove(pdbfile);
 379         process_multiprot_output(fn, &rmsd, &nres2,rotangles,translation,bCountres,countres);
 380         fprintf(fo,"%12.7f",t);
 381         fprintf(fo," %12.7f\n",rmsd);
 382         fprintf(frc,"%12.7f",t);
 383         fprintf(frc,"%12d\n",nres2);
 384         if (bTrjout) {
 385             rotate_conf(natoms,fr.x,NULL,rotangles[XX],rotangles[YY],rotangles[ZZ]);
 386             for(i=0; i<natoms; i++) {
 387                 rvec_inc(fr.x[i],translation);
 388             }
 389             switch(outftp) {
 390                 case efTRJ:
 391                 case efTRR:
 392                 case efG87:
 393                 case efXTC:
 394                     write_trxframe(trxout,&fr,NULL);
 395                     break;
 396                 case efGRO:
 397                 case efG96:
 398                 case efPDB:
 399                     sprintf(out_title,"Generated by do_multiprot : %s t= %g %s",
 400                             title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
 401                     switch(outftp) {
 402                         case efGRO:
 403                             write_hconf_p(out,out_title,&useatoms,prec2ndec(fr.prec),
 404                                           fr.x,NULL,fr.box);
 405                             break;
 406                         case efPDB:
 407                             fprintf(out,"REMARK    GENERATED BY DO_MULTIPROT\n");
 408                             sprintf(out_title,"%s t= %g %s",title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
 409                             /* if reading from pdb, we want to keep the original
 410                                model numbering else we write the output frame
 411                                number plus one, because model 0 is not allowed in pdb */
 412                             if (ftp==efPDB && fr.step > model_nr) {
 413                                 model_nr = fr.step;
 414                             }
 415                             else {
 416                                 model_nr++;
 417                             }
 418                             write_pdbfile(out,out_title,&useatoms,fr.x,ePBC,fr.box,' ',model_nr,NULL,TRUE);
 419                             break;
 420                         case efG96:
 421                             fr.title = out_title;
 422                             fr.bTitle = (nframe == 0);
 423                             fr.bAtoms = FALSE;
 424                             fr.bStep = TRUE;
 425                             fr.bTime = TRUE;
 426                             write_g96_conf(out,&fr,-1,NULL);
 427                     }
 428                     break;
 429             }
 430         }
 431         nframe++;
 432     } while(read_next_frame(oenv,status,&fr));
 433     if (bCountres) {
 434         fres=  xvgropen(opt2fn("-cr",NFILE,fnm),"Number of frames in which the residues are aligned to","Residue","Number",oenv);
 435         for (i=0;i<ratoms.nres;i++) {
 436             fprintf(fres,"%10d  %12d\n",countres[i].resnr,countres[i].count);
 437         }
 438         ffclose(fres);
 439     }
 440     ffclose(fo);
 441     ffclose(frc);
 442     fprintf(stderr,"\n");
 443     close_trj(status);
 444     if (trxout != NULL) {
 445         close_trx(trxout);
 446     }
 447     else if (out != NULL) {
 448         ffclose(out);
 449     }
 450     view_all(oenv,NFILE, fnm);
 451     sfree(xr);
 452     if (bCountres) {
 453         sfree(countres);
 454     }
 455     free_t_atoms(&ratoms,TRUE);
 456     if (bTrjout) {
 457         if (outftp==efPDB || outftp==efGRO || outftp==efG96) {
 458             free_t_atoms(&useatoms,TRUE);
 459         }
 460     }
 461     thanx(stderr);
 462     return 0;
 463 }