src/contrib/do_multiprot.c

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