codeml.c

一个神经网络工具箱

      }

      ReadSeq((com.verbose?fout:NULL),fseq); /*may change seqtype*/
      i=(com.ns*2-1)*sizeof(struct TREEN);
      if((nodes=(struct TREEN*)malloc(i))==NULL) error2("oom nodes");
      /* BootstrapSeq("boot");  exit(0); */

      /* if(com.ngene>1 && com.Mgene==1) printSeqsMgenes (); */
      if(com.ngene>1 && com.aaDist>=FIT1)  /* because of pcodon0[] */
         { error2("ngene for fitness models");}

      pmodel=(com.seqtype==CODONseq?NSbranchmodels[com.model]:aamodels[com.model]);
      fprintf(fout,"%s (in %s)  ",seqtypestr[com.seqtype-1],VerStr);
      fprintf(fout,"  %s   Model: %s ",com.seqf,pmodel);
      if(com.clock) fprintf(fout," %s ",clockstr[com.clock]);
      if(com.seqtype==CODONseq||com.model==FromCodon) {
         if (com.fix_kappa) fprintf(fout, " kappa = %.3f fixed\n", com.kappa);
         if (com.fix_omega) fprintf(fout, " omega = %.3f fixed\n", com.omega);
      }
      if (com.seqtype==AAseq && (com.model==Empirical||com.model==Empirical_F))
         fprintf (fout, "(%s) ", com.daafile);
      if(com.seqtype==AAseq&&com.nrate) fprintf(fout,"(nrate:%d) ",com.nrate);
      if (com.alpha && com.rho) fprintf (fout, "Auto-");
      if (com.alpha) fprintf (fout,"dGamma (ncatG=%d) ", com.ncatG);
      if (com.ngene>1)
         fprintf (fout, " (%d genes: %s)  ", com.ngene, Mgenestr[com.Mgene]);

      if (com.alpha==0) com.nalpha=0;
      else              com.nalpha=(com.nalpha?com.ngene:!com.fix_alpha);
      if (com.Mgene==1) com.nalpha=!com.fix_alpha;
      if(com.nalpha>1 && (!com.alpha || com.ngene==1 || com.fix_alpha))
         error2("Malpha");
      if (com.nalpha>1 && com.rho) error2("Malpha or rho");
      if (com.nalpha>1) fprintf (fout,"(%d gamma)", com.nalpha);
     
      if (com.Mgene && com.ngene==1) error2("Mgene for one gene.");
      if (com.seqtype==CODONseq) {
         fprintf (fout, "\nCodon frequencies: %s\n", codonfreqs[com.codonf]);
         if(com.alpha) 
            fputs("Warning: Gamma model for codons.  See documentation.",fout);
      }
      if ((com.seqtype==CODONseq||com.model==FromCodon) 
         && (com.aaDist &&com.aaDist<10))
         fprintf(fout,"%s, %s\n",com.daafile,(com.aaDist>0?"geometric":"linear"));

      if (com.NSsites) {
         fprintf(fout,"Site-class models");
         if (nnsmodels==1) {
            fprintf(fout," %s",NSsitesmodels[com.NSsites]);
            if(com.NSsites>=NSdiscrete)fprintf(fout," (%d categories)",com.ncatG);
         }
         if(com.nparK) fprintf(fout," & HMM");
         FPN(fout);
         if(com.aaDist)
            fprintf(fout,"\nFitness models: aaDist: %d\n",com.aaDist);
      }

      if(com.clock==2 && 
        (rateBranch=(int*)realloc(rateBranch,com.ns*2*sizeof(int)))==NULL) 
           error2("oom rateBranch");
      com.sspace=max2(800000,3*com.ncode*com.ncode*(int)sizeof(double));
      if(com.NSsites) 
         com.sspace=max2(com.sspace,2*com.ncode*com.ncode+4*com.npatt)*(int)sizeof(double);
/*
      i=com.ns*(com.ns-1)/2;
      com.sspace=max2(com.sspace,
        (int)sizeof(double)*((com.ns*2-2)*(com.ns*2-2+4+i)+i));
*/
      if((com.space=(double*)realloc(com.space,com.sspace))==NULL) {
         printf("\nfailed to get %d bytes for space", com.sspace);
         error2("oom space");
      }
      SeqDistance=(double*)realloc(SeqDistance, i*sizeof(double));
      ancestor=(int*)realloc(ancestor, i*sizeof(int));
      if(SeqDistance==NULL||ancestor==NULL) error2("oom distance&ancestor");

      if(com.seqtype==AAseq) {
         Initialize (fout);
         if (com.model==FromCodon||com.model==AAClasses) 
            AA2Codonf(com.pi, com.fb61);  /* get codon freqs from aa freqs */ 
      }
      else {  /* codon sequences */
         if(com.sspace<max2(com.ngene+1,com.ns)*(64+12+4)*(int)sizeof(double)) {
            com.sspace=max2(com.ngene+1,com.ns)*(64+12+4)*sizeof(double);
            if((com.space=(double*)realloc(com.space,com.sspace))==NULL)
               error2("oom space for InitializeCodon");
         }
         if(InitializeCodon(fout,com.space)) error2("giving up on stop codons");
         if(com.Mgene==3) FOR (i,com.ngene) xtoy(com.pi,com.piG[i],com.ncode);
#if (0)
         /* read in site pattern frequencies as data.  
            Note that codon frequencies are incorrect, and it is better to do this 
            right after PatternWeight(). Also have to change com.ls and 
            com.pose */
         com.pose=(int*) realloc(com.pose, 188*sizeof(int));
         for(i=0,com.ls=0; i<com.npatt; i++) {
            fscanf(fseq,"%lf", &com.fpatt[i]);
            for(k=0; k<(int)com.fpatt[i]; k++) com.pose[com.ls++]=i;
         }
         matout2(F0,com.fpatt,1,com.npatt,4,0); 
         getchar();
#endif
      }
      if(com.seqtype==CODONseq)
         DistanceMatNG86(fout,fpair[0],fpair[1],fpair[2],0);
      else
         DistanceMatAA(fout);

      if(com.ndata==1) fclose(fseq);
      fflush(fout);

      if(com.seqtype==AAseq && com.model==Poisson && !com.print) 
         PatternJC69like (fout);
      if(com.alpha || com.NSsites) {
         s2=com.npatt*com.ncatG*sizeof(double);
         if((com.fhK=(double*)realloc(com.fhK,s2))==NULL) error2("oom fhK");
      }

      if(com.runmode==-2 && com.Mgene!=1) {
         if(com.seqtype==CODONseq) 
            PairwiseCodon(fout,fpair[3],fpair[4],fpair[5],com.space);  
         else
            PairwiseAA(fout, fpair[0]);  
      }
      else {
         if(!com.cleandata) {
            slkl0=com.ns*com.ncode*com.npatt*sizeof(double);
            com.lkl0=(double*)malloc(slkl0);
         }
         com.slkl1= 2 *com.ncode*com.npatt*sizeof(double);
         /* com.slkl1= (com.ns-1)*com.ncode*com.npatt*sizeof(double); */
         com.lkl=(double*)realloc(com.lkl,com.slkl1);

         printf("\n%9ld bytes for distance",com.ns*(com.ns-1)/2*sizeof(double));
         printf("\n%9d bytes for lkl0\n%9d bytes for lkl1\n",slkl0,com.slkl1);
         printf ("%9d bytes for fhK\n%9d bytes for space\n",s2,com.sspace);
         if((!com.cleandata&&com.lkl0==NULL) || com.lkl==NULL)
            error2("oom");

    /*
SlidingWindow(fout,com.space);
exit(0);
*/

         if (nnsmodels>1) {
            printf("\nNSsites batch run.\nUsing ncatG values as in YNGP2000.\n");
            for(insmodel=0; insmodel<nnsmodels; insmodel++) {
               com.NSsites=nsmodels[insmodel];
               if(com.NSsites<=NSselection)  com.ncatG=com.NSsites+1;
               else if (com.NSsites==NSdiscrete) com.ncatG=3;
               else if (com.NSsites==NSfreqs) com.ncatG=5;
               else if (com.NSsites==NSbetaw||com.NSsites==NS02normal) 
                     com.ncatG=ncatG0+1;
               else  com.ncatG=ncatG0;

               com.nrate=com.nkappa=(com.hkyREV?5:!com.fix_kappa);
               if(com.NSsites==0 || com.NSsites==NSselection || com.NSsites==NSbetaw)
                  com.nrate+=!com.fix_omega;
               if(com.NSsites==NSdiscrete)   com.nrate+=com.ncatG;

               printf("\n\nModel %d: %s\n",com.NSsites, NSsitesmodels[com.NSsites]);
               fprintf(fout,"\n\nModel %d: %s",com.NSsites,NSsitesmodels[com.NSsites]);
               fprintf(frst,"\n\nModel %d: %s",com.NSsites,NSsitesmodels[com.NSsites]);
               fprintf(frub,"\n\nModel %d: %s",com.NSsites,NSsitesmodels[com.NSsites]);
               if(com.NSsites) fprintf(fout," (%d categories)",com.ncatG);
               FPN(fout);
               Forestry (fout);

               printf("\nTime used: %s\n", printtime(timestr));
               fprintf(fout,"\nTime used: %s\n", printtime(timestr));
            }
            exit(0);
         }

         if (com.Mgene==1)        MultipleGenes(fout, fpair, com.space);
         else if (com.runmode==0) Forestry(fout);
         else if (com.runmode==3) StepwiseAddition(fout, com.space);
         else if (com.runmode>=4) Perturbation(fout,(com.runmode==4),com.space);
         else                     StarDecomposition(fout, com.space);
      }
      FPN(frst);  fflush(frst);  fflush(frst1);
      free(nodes);
   }  /* for (idata) */
   if (noisy) putchar ('\a');
   fclose(frst);   if(fin)fclose(fin);
   k=0;
   if(com.seqtype==1) k=(com.runmode==-2?6:3);
   else if (com.runmode==-2) k=1;
   FOR(i,k) fclose(fpair[i]);

   if(com.ndata>1 && fseq) fclose(fseq);  
   free(PMat);
   printf("\nTime used: %s\n", printtime(timestr));
   return (0);
}


/* x[]: t[ntime]; rgene[ngene-1]; kappa; p[](NSsites); omega[]; 
        { alpha(for NSsites) !! alpha, rho || rK[], fK[] || rK[], MK[] }
*/

int Forestry (FILE *fout)
{
   static int times=0;
   FILE *ftree, *frate=NULL;
   int  status=0, i,j=0,k, itree, ntree, np, iteration=1;
   int pauptree=0, haslength;
   double x[NP],xb[NP][2], lnL=0,lnL0=0, e=1e-8,tl=0, *var=NULL, nchange=-1;

   if(com.clock==3) GetSeqTimes();

   if ((ftree=fopen(com.treef,"r"))==NULL) {
      printf("\ntree file %s not found.\n", com.treef);
      exit(-1);
   }
   GetTreeFileType(ftree, &ntree, &pauptree, 0);
   if (com.alpha)
      frate=(FILE*)fopen(ratef,"w");
   if (ntree>10 && com.print) puts("\nlarge lnf file");
   FOR(i,NP-NBRANCH) xcom[i]=0.1;
   flnf=fopen("lnf","w+");
   fprintf(flnf,"%6d %6d %6d\n", ntree, com.ls, com.npatt);

   if(com.seqtype==1 && com.aaDist>=FIT1) {
      xtoy(com.pi,pcodon0,64);
      zero(paa0,20);
      FOR(i,com.ncode) paa0[GeneticCode[com.icode][FROM61[i]]]+=pcodon0[i];
      pcodonClass=(double*)malloc(com.ncatG*64*sizeof(double));
      if(pcodonClass==NULL) error2("oom pcodonClass");
   }
   if(!com.cleandata) InitPartialLikelihood();

   for(itree=0; ntree==-1||itree<ntree; iteration=1,itree++) {
      if((pauptree && PaupTreeRubbish(ftree)) || 
         ReadaTreeN(ftree,&haslength, &i,1)) 
            { puts("err or end of tree file."); break; }
      printf("\nTREE # %2d\n", itree+1);
      fprintf(fout,"\nTREE # %2d:  ", itree+1);
      fprintf(flnf,"\n\n%2d\n", itree+1);
      if(com.print) fprintf (frst,"\n\nTREE # %2d\n", itree+1);
      fprintf(frub,"\n\nTREE #%2d\n", itree+1);

      if (com.fix_blength==2 && !haslength) error2("no blengths in tree");
      if (com.fix_blength>0 && !haslength) com.fix_blength=0;
      if (times++==0 && com.fix_blength>0 && haslength) {
         if(com.clock) puts("\nBranch lengths in tree are ignored");
         else {
            if(com.fix_blength==2) puts("\nBranch lengths in tree are fixed.");
            else if(com.fix_blength==1) puts("\nBranch lengths in tree used as initials.");
            if(com.fix_blength==1) {
               FOR(i,tree.nnode) 
                  if((x[nodes[i].ibranch]=nodes[i].branch)<0) 
                     x[nodes[i].ibranch]=1e-5;
            }
         }
      }
      LASTROUND=0;
      if(com.cleandata) nchange=MPScore(com.space);
      if(com.ns<40) { OutaTreeN(F0,0,0); printf("   MP score: %.0f",nchange); }
      OutaTreeN(fout,0,0); fprintf(fout,"   MP score: %.0f",nchange);

      if(!com.clock && nodes[tree.root].nson<=2 && com.ns>2) {
         puts("\nThis is a rooted tree, without clock.  Check.");
         fputs("\nThis is a rooted tree.  Please check!",fout);
      }
      fflush(fout),  fflush(flnf);

      GetInitials(x, &i);
      if(i==-1) iteration=0;
      if(iteration) SetxInitials (x); /* start within the feasible region */

      if((np=com.np)>NP || np-com.ntime>NP-NBRANCH) error2("raise NP");
      if((i=spaceming2(np))>com.sspace) {
         com.sspace=i;
         printf ("\nspace adjusted to %9d bytes\n",com.sspace);
         if((com.space=(double*)realloc(com.space,com.sspace))==NULL) {
            printf("\ntrying to get %d bytes for ming2", com.sspace);
            error2("oom space");
         }
      }
      printf("\nntime & nrate & np:%6d%6d%6d\n",com.ntime,com.nrate,com.np);
      if(itree && !fin) for(i=0;i<np-com.ntime;i++) x[com.ntime+i]=xcom[i];
/*
      if(com.seqtype==CODONseq && com.NSsites==0 && com.model){
         printf("\n%d dN/dS ratios for branches assumed:\n",N_OMEGA);
         FOR(i,tree.nbranch) printf("%4d",OmegaBranch[i]); FPN(F0);
         fprintf (fout, "\n%d dN/dS ratios for branches assumed:\n",N_OMEGA);
         FOR(i,tree.nbranch) fprintf(fout,"%4d",OmegaBranch[i]); FPN(fout);
      }
*/
      if (com.clock==2) {
         printf("\n%d rates for branches assumed:\n",N_rateBranch);
         FOR (i,tree.nbranch) printf("%3d", rateBranch[i]); FPN(F0);
         FPN(fout); FOR(i,tree.nbranch)fprintf(fout,"%3d",rateBranch[i]); 
         FPN(fout);
      }

      if(noisy>=2&&com.npi0)
         printf("\n\a%d zero frequencies.\n", com.npi0);

      PointLklnodes ();
      lnL = com.plfun (x,np);

      if(noisy) {
         printf("\nnp =%6d", np);
         if(noisy>2 && np<100) matout(F0,x,1,np);
         printf("\nlnL0 = %12.6f\n",-lnL);

/*
         gradient (np, x, lnL, com.space, com.plfun, space+np, 1);
         FOR(i,np) printf("%12.6f", com.space[i]);  FPN(F0);
*/
      }

      if(iteration) {
         SetxBound(np,xb);
         SetxInitials (x);
         if(com.method)
            j=minB(noisy>2?frub:NULL, &lnL,x,xb, e*100, com.space);
         else