📄 codeml.c
字号:
OMEGA_FIX=com.omega;
if((com.model==0 && com.NSsites==NSdiscrete)
|| (com.model && com.NSsites && com.NSsites!=NSselection
&&com.NSsites!=NSdiscrete&&com.NSsites!=NSbetaw))
error2("\afix_omega?");
if(com.aaDist) error2("fix_omega & aaDist");
}
if (com.model>NSbranch3) error2("seqtype or model.");
if (com.model==NSbranch3 && com.NSsites!=3) error2("model or NSsites.");
/*
if (com.model==NSbranch2 && com.clock==2)
error2("NSbranch & local clock.");
*/
if(com.kappa<0) error2("kappa..");
if(com.runmode==-2 && (com.NSsites||com.alpha||com.aaDist))
error2("err: incorrect model for pairwise comparison.");
if(com.runmode>0 && com.model==2) error2("tree search & model");
if(com.aaDist && com.NSsites!=0 && com.NSsites!=NSdiscrete)
error2("NSsites && aaDist.");
if(com.aaDist==0) {
if(!com.fix_omega || (com.Mgene && com.Mgene>=3)) com.nrate++;
}
else {
if(com.aaDist<=6) com.nrate+=2; /* a & b, PSB2000 */
else if(com.aaDist==FIT1) com.nrate+=4; /* fitness models: */
else if(com.aaDist==FIT2) com.nrate+=5; /* ap, p*, av, v*, b */
if(com.aaDist>=FIT1) FOR(i,2) FOR(j,20) AAchem[i][j]/=AAchem[i][20];
}
if (com.Mgene>=3 && com.nrate==0) error2("Mgene");
if(com.NSsites) {
com.nrate=(com.hkyREV ? 5 : !com.fix_kappa);
if(com.model && com.ngene>1) error2("NSbranchsites with ngene.");
if(com.model && com.NSsites)
if((com.model!=2 && com.model!=3)
|| (com.NSsites!=NSselection&&com.NSsites!=NSdiscrete))
error2("only NSsites=2,3 & model=2,3 are compatible.");
if(com.alpha || com.fix_alpha==0) error2("NSsites & gamma");
if(com.NSsites<=NSselection) com.ncatG=com.NSsites+1;
if(com.NSsites==NSbetaw || com.NSsites==NS02normal) com.ncatG++;
if(com.model==2) com.ncatG=4; /* branchsite models */
if(com.model==3) { /* Joe's models */
if(com.NSsites==NSselection) com.ncatG=3;
if(com.NSsites==NSdiscrete && com.ncatG!=2 && com.ncatG!=3)
error2("use 2 or 3 for ncatG for model=3?");
}
if(com.NSsites==NSdiscrete && !com.model && com.ncatG>3) {
puts("\nncatG may be too big for M3 (discrete)?");
sleep2(3);
}
if(com.NSsites==NSfreqs && com.ncatG!=5)
{ puts("\nncatG changed to 5."); com.ncatG=5; }
if(com.NSsites==NSselection || com.NSsites==NSbetaw)
{ if(!com.fix_omega) com.nrate++; else OMEGA_FIX=com.omega; }
else if(com.NSsites==NSdiscrete) {
com.nrate+=com.ncatG; /* omega's */
if(com.aaDist) {
if (com.aaDist<=6) com.nrate+=com.ncatG; /* a&b PSB2000 */
else { /* fitness models */
com.nrate=!com.fix_kappa+4*com.ncatG;
if(com.aaDist==FIT2) com.nrate+=com.ncatG;
}
}
}
if(com.model==NSbranch3) com.nrate++;
}
}
}
else
error2 ("seqtype..");
if(com.runmode==-2 && com.cleandata==0) {
com.cleandata=1; puts("gaps are removed for pairwise comparison.");
}
if(com.method &&(com.clock||com.rho))
{ com.method=0; puts("\aiteration method reset"); }
if (com.runmode==3 && (com.clock)) error2("runmode+clock");
if (com.aaDist<=6 && (com.seqtype==CODONseq || com.model==FromCodon))
strcpy(com.daafile, daafiles[abs(com.aaDist)]);
else if (com.seqtype==AAseq && com.model>=REVaa_0)
strcpy(com.daafile,"mtmam.dat");
if (com.fix_alpha && com.alpha==0) {
if (com.rho) puts("rho set to 0."); com.fix_rho=1; com.rho=0;
}
if(!com.fix_alpha && com.alpha<=0) { com.alpha=0.5; puts("alpha reset"); }
if(!com.fix_rho && com.rho==0) { com.rho=0.001; puts("init rho reset"); }
if(com.alpha||com.NSsites)
{ if(com.ncatG<2 || com.ncatG>NCATG) error2("ncatG"); }
else if (com.ncatG>1) com.ncatG=1;
#ifdef NSSITES_K1_K2_CLASSES
if((com.NSsites==NSgamma||com.NSsites==NS2gamma||com.NSsites>=NSbetagamma)
&& com.ncatG<10) error2("need more categories for NSsites");
#endif
fin=fopen("in.codeml","r");
return (0);
}
int testx (double x[], int np)
{
int i,k;
double tb[]={.4e-5, 29}, rgeneb[]={0.1,20}, rateb[]={1e-5,999}, omegab=.005;
double alphab[]={0.005,99}, rhob[]={0.01,0.99};
FOR (i,(com.clock?1:com.ntime)) if (x[i]<tb[0] || x[i]>tb[1]) return (-1);
if(com.clock)
for(i=1;i<com.ntime;i++) if (x[i]<1e-6 || x[i]>1-1e-6) return (-1);
if (np==com.ntime) return (0);
FOR (i,com.nrgene)
if (x[com.ntime+i]<rgeneb[0] || x[com.ntime+i]>rgeneb[1]) return (-1);
FOR (i,com.nrate) if (x[com.ntime+com.nrgene+i]>rateb[1]) return (-1);
if (com.seqtype==CODONseq && com.nrate==2) {
if(x[com.ntime+com.nrgene]<rateb[0] || x[com.ntime+com.nrgene+1]<omegab)
return (-1);
}
else
FOR (i,com.nrate) if (x[com.ntime+com.nrgene+i]<rateb[0]) return (-1);
k=com.ntime+com.nrgene+com.nrate;
for (i=0; i<com.nalpha; i++,k++)
if (x[k]<alphab[0] || x[k]>alphab[1]) return (-1);
if (!com.fix_rho) if (x[np-1]<rhob[0] || x[np-1]>rhob[1]) return (-1);
return(0);
}
int SetxInitials (double x[])
{
/* This forces initial values into the boundary of the space
*/
int i, k;
double rateb[]={1e-5,89}; /* what are the negative rate? */
double tb[]={.0001,25}, rgeneb[]={0.1,9};
double alphab[]={0.05,9}, rhob[]={0.01,0.9};
FOR (i,com.ntime)
{ if (x[i]>tb[1]) x[i]=tb[1]; if (x[i]<tb[0]) x[i]=tb[0]; }
FOR (i,com.nrgene) {
if (x[com.ntime+i]>rgeneb[1]) x[com.ntime+i]=rgeneb[1];
if (x[com.ntime+i]<rgeneb[0]) x[com.ntime+i]=rgeneb[0];
}
for (i=0,k=com.ntime+com.nrgene; i<com.nrate; i++,k++) {
if (x[k]>rateb[1]) x[k]=rateb[1];
if (x[k]<rateb[0]) x[k]=rateb[0];
}
for (i=0; i<com.nalpha; i++,k++) {
if (x[k]>alphab[1]) x[k]=alphab[1];
if (x[k]<alphab[0]) x[k]=alphab[0];
}
if (!com.fix_rho)
{ if (x[k]>rhob[1]) x[k]=rhob[1]; if (x[k]<rhob[0]) x[k]=rhob[0]; }
/*
FOR (i,com.np) if (x[i]<1e-4) printf("x[%d]=%.6f\n", i,x[i]);
*/
return(0);
}
int SetxBound (int np, double xb[][2])
{
int i,j,k, K=com.ncatG;
double tb[]={4e-6,50}, rgeneb[]={0.01,99}, rateb[]={1e-4,99};
double alphab[]={0.005,99}, betab[]={.001,99}, omegab[]={.0001,999};
double rhob[]={0.01,0.99}, pb[]={.000001,.999999};
if(com.NSsites) omegab[0]=0.00001;
if(com.clock) {
xb[0][0]=.00005; xb[0][1]=tb[1];
for(i=1;i<tree.nnode-com.ns;i++) /* for node times, clock=1 or 2 */
FOR(j,2) { xb[i][0]=.001; xb[i][1]=.9999; }
for(; i<com.ntime; i++) /* for rates for branches, clock=2 */
FOR(j,2) { xb[i][0]=tb[0]; xb[i][1]=tb[1]; }
}
else
FOR (i,com.ntime) FOR (j,2) xb[i][j]=tb[j];
for(i=com.ntime;i<np;i++) { xb[i][0]=rateb[0]; xb[i][1]=rateb[1]; }
FOR (i,com.nrgene) FOR (j,2) xb[com.ntime+i][j]=rgeneb[j];
FOR (i,com.nrate) FOR (j,2) xb[com.ntime+com.nrgene+i][j]=rateb[j];
k=com.ntime+com.nrgene+!com.fix_kappa;
if (com.NSsites) { /* p's before w's in xb[] */
switch(com.NSsites) {
case(NSneutral): xb[k][0]=pb[0]; xb[k++][1]=pb[1]; break; /* p0 */
case(NSselection): /* for p0, p1, w2 */
FOR(j,2) { xb[k][0]=-99; xb[k++][1]=99; }
if(!com.fix_omega) { xb[k][0]=omegab[0]; xb[k++][1]=omegab[1]; }
break;
case(NSdiscrete): /* pK[] & rK[] */
if(com.model) K=3;
if(com.nparK==0)
FOR(j,K-1) { xb[k][0]=-99; xb[k++][1]=99; }
FOR(j,K) { xb[k][0]=omegab[0]; xb[k++][1]=omegab[1]; }
if(com.nparK)
FOR(j,K*(K-1)) { xb[k][0]=-99; xb[k++][1]=99; }
if(com.seqtype==CODONseq&&com.aaDist)
FOR(j,K) { xb[k][0]=omegab[0]; xb[k++][1]=omegab[1]; }
break;
case(NSfreqs): /* p0...pK */
FOR(j,K-1) { xb[k][0]=-99; xb[k++][1]=99; }
break;
case(NSgamma):
FOR(j,2) { xb[k][0]=alphab[0]; xb[k++][1]=alphab[1]; } break;
case(NS2gamma): /* p0, alpha1,beta1,alpha2=beta2 */
xb[k][0]=pb[0]; xb[k++][1]=pb[1];
FOR(j,3) { xb[k][0]=alphab[0]; xb[k++][1]=alphab[1]; }
break;
case(NSbeta): /* p_beta,q_beta */
FOR(j,2) { xb[k][0]=betab[0]; xb[k++][1]=betab[1]; }
break;
case(NSbetaw): /* p0, p_beta, q_beta, w */
xb[k][0]=pb[0]; xb[k++][1]=pb[1]; /* p0 */
FOR(j,2) { xb[k][0]=betab[0]; xb[k++][1]=betab[1]; } /* p & q */
if(!com.fix_omega) { xb[k][0]=omegab[0]; xb[k++][1]=omegab[1]; }
break;
case(NSbetagamma): /* p0, p_beta, q_beta, alpha, beta */
xb[k][0]=pb[0]; xb[k++][1]=pb[1]; /* p0 */
FOR(j,4) { xb[k][0]=betab[0]; xb[k++][1]=betab[1]; } /* p&q, a&b */
break;
case(NSbeta1gamma): /* p0, p_beta, q_beta, alpha, beta */
xb[k][0]=pb[0]; xb[k++][1]=pb[1]; /* p0 */
FOR(j,4) { xb[k][0]=betab[0]; xb[k++][1]=betab[1]; } /* p&q, a&b */
break;
case(NSbeta1normal): /* p0, p_beta, q_beta, mu, s */
xb[k][0]=pb[0]; xb[k++][1]=pb[1]; /* p0 */
FOR(j,4) { xb[k][0]=betab[0]; xb[k++][1]=betab[1]; } /* p&q, mu&s */
xb[k-2][0]=1; xb[k-2][1]=9; /* mu */
break;
case(NS02normal): /* p0, p1, mu2, s1, s2 */
FOR(j,2) { xb[k][0]=pb[0]; xb[k++][1]=pb[1]; } /* p0 & p1, */
FOR(j,3) { xb[k][0]=.0001; xb[k++][1]=29; } /* mu2,s1,s2 */
break;
case(NS3normal): /* p0, p1, mu2, s0, s1, s2 */
FOR(j,2) { xb[k][0]=-99; xb[k++][1]=99; } /* p0 & p1, tranformed */
FOR(j,4) { xb[k][0]=.0001; xb[k++][1]=29; } /* mu2,s0,s1,s2 */
break;
}
}
else if((com.seqtype==CODONseq||com.model==FromCodon)&&com.model!=AAClasses)
{ if(!com.fix_omega) { xb[k][0]=omegab[0]; xb[k][1]=omegab[1]; } }
if(com.seqtype==CODONseq&&com.model)
FOR(j,N_OMEGA-com.fix_omega)
{ xb[k+j][0]=omegab[0]; xb[k+j][1]=omegab[1]; }
if (com.aaDist<0 && (com.seqtype==1||com.model==FromCodon)) {
/* linear relationship between d_ij and w_ij */
if(com.nrate != !com.fix_kappa+1+(com.seqtype==1)) error2("in Setxbound");
xb[com.ntime+com.nrgene+!com.fix_kappa][1]=1; /* 0<b<1 */
}
k=com.ntime+com.nrgene+com.nrate;
for (i=0;i<com.nalpha;i++,k++) FOR (j,2) xb[k][j]=alphab[j];
if (!com.fix_rho) FOR (j,2) xb[np-1][j]=rhob[j];
/*
if(noisy>=3 && np<20) {
puts("\nBounds:\n");
FOR(i,np) printf(" %10.6f", xb[i][0]); FPN(F0);
FOR(i,np) printf(" %10.6f", xb[i][1]); FPN(F0);
}
*/
return(0);
}
void getpcodonClass(double x[], double pcodonClass[])
{
/* This uses pcodon0[], paa0[], and x[] to calculate pcodonclass[] and
com.pi[] for the fitness models.
pcodon0[] has the codon frequencies observed (codonFreq=3) or expected
(codonFreq=2 or 1 or 0) rootally. Under the fitness models, the expected
codon frequencies pcodonClass[] differs among site classes and from the
rootal pi[] (pcodon0[]).
This is called by SetParameters().
*/
int i,iclass,iaa, k, nclass=(com.NSsites==0?1:com.ncatG);
double paaClass[20], *w,fit;
if(com.seqtype!=1 || com.aaDist<FIT1) error2("getpcodonClass");
k=com.ntime+com.nrgene+!com.fix_kappa+nclass-1;
FOR(iclass, nclass) {
w=x+k+iclass*(4+(com.aaDist==FIT2));
FOR(iaa,20) {
fit=-w[0]*square(AAchem[0][iaa]-w[1])
-w[2]*square(AAchem[1][iaa]-w[3]);
paaClass[iaa]=exp(2*fit);
}
abyx(1/sum(paaClass,20), paaClass, 20);
FOR(i,com.ncode) {
iaa=GeneticCode[com.icode][FROM61[i]];
pcodonClass[iclass*64+i]=pcodon0[i]/paa0[iaa]*paaClass[iaa];
}
if(fabs(1-sum(pcodonClass+iclass*64,com.ncode))>1e-5) error2("pcodon!=1");
/*
fprintf(frst,"\nSite class %d: ",iclass+1);
matout (frst,paaClass,2, 10);
matout (frst,pcodonClass+iclass*64,16,4);
*/
}
if(nclass==1) FOR(i,com.ncode) com.pi[i]=pcodonClass[i];
}
int GetInitials(double x[], int* fromfile)
{
/* This caculates the number of parameters (com.np) and get initial values.
Perhaps try to restruct the code and make two sections for amino ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -