Skip to content
Commits on Source (5)
Expand all Hide whitespace changes
Inline Side-by-side
debian/changelog
View file @ bc993320
paml (4.9i+dfsg-1) UNRELEASED; urgency=medium
* New upstream version
* debhelper-compat 12
* Standards-Version: 4.4.0
-- Andreas Tille <tille@debian.org> Thu, 01 Aug 2019 19:19:13 +0200
paml (4.9h+dfsg-1) unstable; urgency=medium
* New upstream version
......
debian/control
View file @ bc993320
......@@ -5,8 +5,8 @@ Uploaders: Pjotr Prins <pjotr.debian@thebird.nl>,
Andreas Tille <tille@debian.org>
Section: science
Priority: optional
Build-Depends: debhelper (>= 11~)
Standards-Version: 4.1.4
Build-Depends: debhelper-compat (= 12)
Standards-Version: 4.4.0
Vcs-Browser: https://salsa.debian.org/med-team/paml
Vcs-Git: https://salsa.debian.org/med-team/paml.git
Homepage: http://abacus.gene.ucl.ac.uk/software/paml.html
......
doc/pamlHistory.txt
View file @ bc993320
......@@ -8,6 +8,41 @@ discussion site
https://groups.google.com/forum/#!forum/pamlsoftware.
Version 4.9i, March 2019
(*) mcmctree: I have added an option (duplication = 1) for dating a
tree with both speciations and gene duplications, so that some nodes
on the tree share divergence times. Nodes sharing ages are identified
using labels in the tree file: #1, #2, .... I have yet to update the
document about specification of the model.
(*) mcmctree: The TipDate option was written for one locus or
partition and never worked for more than two loci/partitions. I have
edited the code so that it works for multiple partitions, some of
which may be molecular and the others morphological.
(*) codeml: The option estFreq = 0 when codonFreq = 6 (FMutSel0) and 7
(FMutSel) is not working in versions 4.9g and 4.9h. This is fixed
now. This option uses the observed codon or amino acid frequencies
for the mutation-selection models of codon usage. Instead the program
estimates the frequencies using maximum likeihood, which is what the
option estFreq = 1 does. Look at the README file in the
examples/mtCDNAape/ folder.
(*) codeml clade model D: The bounds for the w (dN/dS) ratios in the
first site classes are set tp (0.0001, 0.5) for w0 and (0.5, 1.5) for
w1, in versions 4.9b,c,d,e,f,g, since I added the BEB calculation for
clade model D in 4.9b. The motivation for the bounds is that site
class 0 represents strong purifying selection with a small w0, while
site class 1 should include sites under weak purifying selection with
a larger w1. However the bounds are arbitrary. In some datasets, the
MLEs are found to be at the bounds, making the interpretation awkward.
I have changed the bounds to the following:
w0b[]={0.0001, 1.0}, w1b[]={0.01, 1.5}.
This means that the user should swap the estimates of w0 and w2 if w0 > w1.
Version 4.9h, March 2018
(*) mcmctree: gamma-Dirichlet versus conditional i.i.d. priors for
......
doc/pamlHistory.txt~ 0 → 100644
View file @ bc993320
This diff is collapsed.
examples/DatingSoftBound/mtCDNApri.trees
View file @ bc993320
7 1
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla),
(orangutan, sumatran)) '>.12<.16', gibbon);
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);
//end of file
((((human, (chimpanzee, bonobo)) 'B(.06, .08)', gorilla),
(orangutan, sumatran)) 'B(.12, .16)', gibbon);
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);
((((human, (chimpanzee, bonobo)) 'B(.06, .08)', gorilla), (orangutan, sumatran)) 'B(.12, .16)', gibbon);
examples/TipDate.FluH1/baseml.ctl
View file @ bc993320
......@@ -18,7 +18,7 @@
fix_kappa = 0 * 0: estimate kappa; 1: fix kappa at value below
kappa = 2 * initial or fixed kappas
fix_alpha = 0 * 0: estimate alpha; 1: fix alpha at value below
fix_alpha = 1 * 0: estimate alpha; 1: fix alpha at value below
alpha = 0.25 * initial or fixed alpha, 0:infinity (constant rate)
ncatG = 5 * # of categories in the dG, AdG, or nparK models of rates
fix_rho = 1 * 0: estimate rho; 1: fix rho at value below
......
examples/TipDate.FluH1/mcmctreeClock1.ctl
View file @ bc993320
seed = -1357
seed = 1357
seqfile = H1.txt
* seqfile = ../H1.CP123.txt
treefile = H1.tre
......@@ -7,7 +7,7 @@
seqtype = 0 * 0: nucleotides; 1:codons; 2:AAs
noisy = 3
usedata = 0 in.BV.HKYG5 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
usedata = 2 in.BV.HKYG5 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
ndata = 1 *
clock = 1 * 1: global clock; 2: independent rates; 3: correlated rates
......
examples/TipDate.HIV2/FigTree.tre 0 → 100644
View file @ bc993320
#NEXUS
BEGIN TREES;
UTREE 1 = (((SL92eSM1992: 0.381307, (SL92fSM1992: 0.154397, SL92aSM1992: 0.154397) [&95%={0.0347455, 0.522379}]: 0.226910) [&95%={0.0904505, 0.907043}]: 0.504653, (((Lib1sm1988: 0.142534, FO784h1989: 0.152534) [&95%={0.0748715, 0.520538}]: 0.202566, SIVMNE1982: 0.285100) [&95%={0.152587, 0.820146}]: 0.357506, (SIVSTM1989: 0.615921, ((((((JA1990: 0.060047, ON1990: 0.060047) [&95%={0.0520286, 0.251389}]: 0.066570, EHO1988: 0.106617) [&95%={0.0769271, 0.358953}]: 0.073194, UC1h1988: 0.179811) [&95%={0.105804, 0.465565}]: 0.087350, D205h1986: 0.247161) [&95%={0.152152, 0.584476}]: 0.111152, 2238h1989: 0.388313) [&95%={0.223053, 0.723296}]: 0.156536, (((D194h1987: 0.111214, UC2h1988: 0.121214) [&95%={0.0835921, 0.423991}]: 0.141317, BEN1987: 0.252531) [&95%={0.124754, 0.605509}]: 0.202257, (CAM2CG1987: 0.390471, (ROD1985: 0.316609, ((P07h1995: 0.259167, (P09h1995: 0.168184, NIHZ1986: 0.078184) [&95%={0.0925188, 0.329415}]: 0.090983) [&95%={0.120916, 0.445566}]: 0.110770, (P08h1995: 0.318209, (P010h1995: 0.272484, ((P01h1995: 0.188251, ((P02h1995: 0.071261, P06h1995: 0.071261) [&95%={0.0031347, 0.190992}]: 0.077251, (P05h1995: 0.102844, ALI1989: 0.042844) [&95%={0.0612933, 0.200128}]: 0.045668) [&95%={0.0744984, 0.27015}]: 0.039739) [&95%={0.0946369, 0.324612}]: 0.043111, (P04h1995: 0.107655, P03h1995: 0.107655) [&95%={0.0046477, 0.272953}]: 0.123707) [&95%={0.121765, 0.381337}]: 0.041122) [&95%={0.149241, 0.434521}]: 0.045725) [&95%={0.182706, 0.493541}]: 0.051729) [&95%={0.223594, 0.558168}]: 0.046672) [&95%={0.255413, 0.624248}]: 0.053862) [&95%={0.292723, 0.698665}]: 0.064317) [&95%={0.337601, 0.788576}]: 0.070062) [&95%={0.379475, 0.895678}]: 0.071072) [&95%={0.41195, 1.01957}]: 0.096685) [&95%={0.44965, 1.207}]: 0.143354) [&95%={0.48519, 1.51936}]: 0.312910, (SL2h1993: 0.423337, (SL92bSM1992: 0.160871, SL92cSM1992: 0.160871) [&95%={0.0348426, 0.554373}]: 0.252466) [&95%={0.0928551, 1.06356}]: 0.785533) [&95%={0.51691, 1.97612}];
END;
[Note for FigTree: Under Time Scale, set Offset = 1995.0, Scale factor = -100.0
Untick Scale Bar, & tick Tip Labels, Node Bars, Scale Axis, Reverse Axis, Show Grid.]
src/Makefile.MSVC → src/Makefile.MSC
View file @ bc993320
PRGS = baseml codeml basemlg mcmctree pamp evolver yn00 chi2
CC = cl # cc, gcc, cl
CFLAGS = -O2 -Ot
CFLAGS = -Ox
DEFINE = #
LIBS = #-lm -lM
......
src/baseml.c
View file @ bc993320
......@@ -30,7 +30,7 @@ int Forestry(FILE *fout);
void DetailOutput(FILE *fout, double x[], double var[]);
int GetOptions(char *ctlf);
int GetInitials(double x[], int *fromfile);
int GetInitialsTime(double x[]);
int GetInitialsTimes(double x[]);
int SetxInitials(int np, double x[], double xb[][2]);
int SetParameters(double x[]);
int SetPGene(int igene, int _pi, int _UVRoot, int _alpha, double x[]);
......@@ -52,7 +52,7 @@ int GroupDistances();
struct CommonInfo {
unsigned char *z[NS];
char *spname[NS], seqf[512], outf[512], treef[512], cleandata;
char *spname[NS], seqf[2048], outf[2048], treef[2048], cleandata;
char oldconP[NNODE]; /* update conP for nodes to save computation (0 yes; 1 no) */
int seqtype, ns, ls, ngene, posG[NGENE + 1], lgene[NGENE], *pose, npatt, readpattern;
int np, ntime, nrgene, nrate, nalpha, npi, nhomo, ncatG, ncode, Mgene;
......@@ -72,14 +72,16 @@ struct CommonInfo {
int fix_omega;
double omega;
} com;
struct TREEB {
int nbranch, nnode, root, branches[NBRANCH][2];
double lnL;
} tree;
struct TREEN {
int father, nson, sons[MAXNSONS], ibranch, ipop;
double branch, age, *pkappa, pi[4], *conP, label;
char *nodeStr, fossil;
double branch, age, *pkappa, pi[4], *conP, label, label2;
char fossil, *name, *annotation;
} *nodes, **gnodes, nodes_t[2 * NS - 1];
......@@ -92,7 +94,7 @@ struct SPECIESTREE {
struct TREESPN {
char name[LSPNAME + 1], fossil, usefossil; /* fossil: 0, 1, 2, 3 */
int father, nson, sons[2];
double age, pfossil[7]; /* lower and upper bounds or alpha & beta */
double label, age, pfossil[7]; /* lower and upper bounds or alpha & beta */
double *lnrates; /* log rates for loci */
} nodes[2 * NS - 1];
} stree;
......@@ -129,16 +131,16 @@ int N_PMatUVRoot = 0;
#include "treesub.c"
#include "treespace.c"
int main(int argc, char *argv[])
int main (int argc, char *argv[])
{
FILE *fout, *fseq = NULL, *fpair[6];
char pairfs[1][32] = { "2base.t" };
char rstf[512] = "rst", ctlf[512] = "baseml.ctl", timestr[64];
char rstf[2048] = "rst", ctlf[2048] = "baseml.ctl", timestr[64];
char *Mgenestr[] = { "diff. rate", "separate data", "diff. rate & pi",
"diff. rate & kappa", "diff. rate & pi & kappa" };
int getdistance = 1, i, idata;
size_t s2 = 0;
if (argc > 2 && !strcmp(argv[argc - 1], "--stdout-no-buf"))
setvbuf(stdout, NULL, _IONBF, 0);
......@@ -266,7 +268,13 @@ int main(int argc, char *argv[])
fprintf(fout, "\n\nPrinting out site pattern counts\n");
printPatterns(fout);
}
if (com.verbose == 2) {
if (com.verbose >= 2) {
fprintf(fout, "\nSite-to-pattern map: ");
for (i = 0; i < com.ls; i++)
fprintf(fout, " %2d", com.pose[i] + 1);
fprintf(fout, "\n");
fprintf(fout, "\n**** Alignment mutated for JC69-like models ****\n");
fprintf(fout, "\n%6d %6d\n", com.ns, com.ls);
for (i = 0; i < com.ns; i++) {
fprintf(fout, "\n%-30s ", com.spname[i]);
......@@ -410,9 +418,11 @@ int Forestry(FILE *fout)
fprintf(flnf, "%6d %6d %6d\n", ntree, com.ls, com.npatt);
for (itree = 0; ntree == -1 || itree < ntree; itree++, iteration = 1) {
if (ReadTreeN(ftree, &haslength, &i, 0, 1)) {
if (ReadTreeN(ftree, &haslength, 0, 1)) {
puts("\nend of tree file."); break;
}
ProcessNodeAnnotation(&i);
if (noisy) printf("\nTREE # %2d: ", itree + 1);
fprintf(fout, "\nTREE # %2d: ", itree + 1);
fprintf(frub, "\n\nTREE # %2d\n", itree + 1);
......@@ -490,7 +500,9 @@ int Forestry(FILE *fout)
else status = 0;
}
if (itree == 0) { lnL0 = lnLbest = lnL; btree = 0; }
if (itree == 0) {
lnL0 = lnLbest = lnL; btree = 0;
}
else if (lnL < lnLbest) {
lnLbest = lnL; btree = itree;
}
......@@ -623,12 +635,12 @@ int Forestry(FILE *fout)
void PrintBestTree(FILE *fout, FILE *ftree, int btree)
{
int itree, ntree, i, j;
int itree, ntree, i;
rewind(ftree);
GetTreeFileType(ftree, &ntree, &i, 0);
for (itree = 0; ntree == -1 || itree < ntree; itree++) {
if (ReadTreeN(ftree, &i, &j, 0, 1)) {
if (ReadTreeN(ftree, &i, 0, 1)) {
puts("\nend of tree file."); break;
}
if (itree == btree)
......@@ -1443,7 +1455,7 @@ int ConditionalPNode(int inode, int igene, double x[])
if (com.clock < 5) {
if (com.clock) t *= GetBranchRate(igene, (int)nodes[ison].label, x, NULL);
else t *= com.rgene[igene];
else t *= com.rgene[igene];
}
GetPMatBranch(PMat, x, t, ison);
......@@ -1473,28 +1485,27 @@ int ConditionalPNode(int inode, int igene, double x[])
return (0);
}
int PMatCijk(double P[], double t)
{
/* P(t)ij = SUM Cijk * exp{Root*t}
*/
/* P(t)ij = SUM Cijk * exp{Root*t}
*/
int i, j, k, n = com.ncode, nr = nR;
double expt[5], pij;
if (t < -.001 && noisy>3)
printf("\nt = %.5f in PMatCijk", t);
if (t < 1e-200) { identity(P, n); return(0); }
for (k = 1; k < nr; k++) expt[k] = exp(t*Root[k]);
for (i = 0; i < n; i++) for (j = 0; j < n; j++) {
for (k = 1, pij = Cijk[i*n*nr + j*nr + 0]; k < nr; k++)
pij += Cijk[i*n*nr + j*nr + k] * expt[k];
P[i*n + j] = (pij > 0 ? pij : 0);
double exptm1[5] = {0};
memset(P, 0, n*n * sizeof(double));
for (k = 1; k < nr; k++) exptm1[k] = expm1(t*Root[k]);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
for (k = 0; k < nr; k++)
P[i*n + j] += Cijk[i*n*nr + j*nr + k] * exptm1[k];
}
P[i*n + i] ++;
}
return (0);
}
#ifdef UNDEFINED
int CollapsSite(FILE *fout, int nsep, int ns, int *ncat, int SiteCat[])
......
src/basemlg.c
View file @ bc993320
......@@ -40,7 +40,7 @@ int lfunG_d (double x[], double *lnL, double dl[], int np);
int lfunG_dd (double x[], double *lnL, double dl[], double ddl[], int np);
struct CommonInfo {
unsigned char *z[NS], *spname[NS], seqf[96],outf[96],treef[96];
unsigned char *z[NS], *spname[NS], seqf[2048],outf[2048],treef[2048];
int seqtype, ns, ls, ngene, posG[NGENE+1],lgene[NGENE],*pose,npatt, readpattern;
int clock,fix_alpha,fix_kappa,fix_rgene,Malpha,print,verbose;
int model, runmode, cleandata, ndata;
......@@ -58,8 +58,8 @@ struct TREEB {
struct TREEN {
int father, nson, sons[MAXNSONS], ibranch;
double branch, age, label, *conP;
char *nodeStr, fossil;
double branch, age, label, label2, *conP;
char *annotation, fossil;
} nodes[2*NS-1];
static int nR=4, CijkIs0[64];
......@@ -83,7 +83,7 @@ int main (int argc, char *argv[])
{
FILE *fout, *fseq=NULL, *fpair[6];
char pairfs[1][32]={"2base.t"};
char ctlf[96]="baseml.ctl";
char ctlf[2048]="baseml.ctl";
double space[50000];
noisy=2; com.cleandata=1; /* works with clean data only */
......@@ -158,7 +158,7 @@ int Forestry (FILE *fout, double space[])
fprintf(frub,"\n\nTREE # %2d", itree+1);
fprintf(frate,"\nTREE # %2d\n", itree+1);
if(ReadTreeN(ftree,&haslength,&i,0,1))
if (ReadTreeN(ftree, &haslength, 0, 1))
{ puts("err or end of tree file."); break; }
com.ntime = com.clock ? tree.nnode-com.ns: tree.nbranch;
......@@ -248,7 +248,7 @@ int SetBranch (double x[])
int testx (double x[], int np)
{
int i,k;
double tb[]={1e-5,4}, rgeneb[]={.01,20}, kappab[]={0,80}, alphab[]={.01,99};
double tb[]={1e-5,99}, rgeneb[]={.01,20}, kappab[]={0,80}, alphab[]={.01,999};
if (SetBranch(x)) return (-1);
FOR (i,com.ntime) if (x[i]<tb[0] || x[i]>tb[1]) return (-1);
......@@ -266,8 +266,8 @@ int SetxBound (int np, double xb[][2])
/* sets lower and upper bounds for variables during iteration
*/
int i,j, k=com.ntime+com.nrgene+com.nrate;
double tb[]={.4e-5, 20}, rgeneb[]={1e-4,99}, rateb[]={1e-5,999};
double alphab[]={0.005, 99}, pb[2]={1e-5,1-1e-5};
double tb[]={.4e-5, 99}, rgeneb[]={1e-4,99}, rateb[]={1e-5,999};
double alphab[]={0.005, 999}, pb[2]={1e-5,1-1e-5};
if(com.clock) {
xb[0][0]=tb[0]; xb[0][1]=tb[1];
......
src/codeml.c
View file @ bc993320
......@@ -106,7 +106,7 @@ double CDFLogis(double x, double m, double s);
struct common_info {
unsigned char *z[NS];
char *spname[NS], seqf[512], outf[512], treef[512], daafile[512], cleandata;
char *spname[NS], seqf[2048], outf[2048], treef[2048], daafile[2048], cleandata;
char oldconP[NNODE]; /* update conP for nodes? to save computation */
int seqtype, ns, ls, ngene, posG[NGENE + 1], lgene[NGENE], npatt, *pose, readpattern;
int runmode, clock, verbose, print, codonf, aaDist, model, NSsites;
......@@ -132,14 +132,16 @@ struct common_info {
& eigenQcodon. Try to remove them? */
double *pomega, pkappa[5], *ppi;
} com;
struct TREEB {
int nbranch, nnode, root, branches[NBRANCH][2];
double lnL;
} tree;
struct TREEN {
int father, nson, sons[MAXNSONS], ibranch, ipop;
double branch, age, omega, *conP, label;
char *nodeStr, fossil, usefossil;
double branch, age, omega, *conP, label, label2;
char fossil, usefossil, *name, *annotation;
} *nodes, **gnodes, nodes_t[2 * NS - 1];
......@@ -152,7 +154,7 @@ struct SPECIESTREE {
struct TREESPN {
char name[LSPNAME + 1], fossil, usefossil; /* fossil: 0, 1, 2, 3 */
int father, nson, sons[2];
double age, pfossil[7]; /* lower and upper bounds or alpha & beta */
double label, age, pfossil[7]; /* lower and upper bounds or alpha & beta */
double *lnrates; /* log rates for loci */
} nodes[2 * NS - 1];
} stree;
......@@ -245,13 +247,13 @@ double lnLmodel;
int main(int argc, char *argv[])
{
FILE *fseq = NULL, *fpair[6] = {NULL};
FILE *fseq = NULL, *fpair[6] = { NULL };
char pairfs[6][32] = { "2NG.dS","2NG.dN","2NG.t", "2ML.dS","2ML.dN","2ML.t" };
char ctlf[96] = "codeml.ctl", *pmodel, timestr[64];
char *seqtypestr[3] = { "CODONML", "AAML", "CODON2AAML" };
char *Mgenestr[] = { "diff. rate", "separate data", "diff. rate & pi",
"diff. rate & k&w", "diff. rate & pi & k&w" };
int getdistance = 0, i,j, k, idata, nc, nUVR, cleandata0;
int getdistance = 0, i, j, k, idata, nc, nUVR, cleandata0;
size_t s2;
starttimer();
......@@ -521,7 +523,7 @@ int main(int argc, char *argv[])
else {
if (com.seqtype == AAseq && com.model == FromCodon0) {
SetMapAmbiguity(com.seqtype, 1);
com.seqtype = CODONseq;
com.seqtype = CODONseq;
/* com.model = 0; */
com.cleandata = 0;
com.ncode = Nsensecodon;
......@@ -597,7 +599,7 @@ int main(int argc, char *argv[])
} /* for (idata) */
fclose(frst);
for (i = 0; i < 6; i++) if(fpair[i]) fclose(fpair[i]);
for (i = 0; i < 6; i++) if (fpair[i]) fclose(fpair[i]);
if (com.ndata > 1 && fseq) fclose(fseq);
fclose(fout); fclose(frub);
if (finitials) fclose(finitials);
......@@ -643,10 +645,11 @@ int Forestry(FILE *fout)
}
for (itree = 0; ntree == -1 || itree < ntree; itree++, iteration = 1) {
if (ReadTreeN(ftree, &haslength, &i, 0, 1)) {
if (ReadTreeN(ftree, &haslength, 0, 1)) {
puts("end of tree file."); break;
}
if (com.seqtype == 1 && com.model)
ProcessNodeAnnotation(&k);
printf("\nTREE # %2d\n", itree + 1);
fprintf(fout, "\n\nTREE # %2d: ", itree + 1);
fprintf(flnf, "\n\n%2d\n", itree + 1);
......@@ -730,8 +733,8 @@ int Forestry(FILE *fout)
printf("Out..\nlnL = %12.6f\n", -lnL);
printf("%d lfun, %d eigenQcodon, %d P(t)\n", NFunCall, NEigenQ, NPMatUVRoot);
if (itree == 0) {
lnL0 = lnL;
if (itree == 0) {
lnL0 = lnL;
for (i = 0; i < np - com.ntime; i++) xcom[i] = x[com.ntime + i];
}
else if (!j)
......@@ -903,15 +906,13 @@ double *PointOmega(double xcom[], int igene, int inode, int isiteclass)
com.NSsites and igene. This sometimes points to com.omega or com.rK[].
This is called by SetParameters(), DetailOutput(), etc.
Ziheng Notes: 8 August 2003.
Difficulties in using this with lfunt() etc.
Trying to remove global variables com.pomega and com.pkappa through
PointOmega and PointKappa, but was unsuccessful when too many changes were
made at the same time. Perhaps look at this later. Note that some
variables are passed over the head from lfunt() etc. to eigenQcodon().
Ziheng Notes: 8 August 2003.
*/
int k = com.nrgene + com.nkappa, backfore;
int nka = (com.hkyREV ? 5 : 1), nw = (com.aaDist == AAClasses ? com.nOmegaType : 1);
......@@ -1308,7 +1309,7 @@ void DetailOutput(FILE *fout, double x[], double var[])
if (com.seqtype == CODONseq && com.ngene == 1 && (com.model == 0 || com.model == FromCodon0 || com.NSsites == 0)) {
tdSdNb = (double*)malloc(tree.nnode * 3 * sizeof(double));
if (tdSdNb == NULL) error2("oom DetailOutput");
if (com.model>=NSbranchB && com.model <= NSbranch3 && com.aaDist != AAClasses) { /* branch models */
if (com.model >= NSbranchB && com.model <= NSbranch3 && com.aaDist != AAClasses) { /* branch models */
fprintf(fout, "\nw (dN/dS) for branches: ");
k = com.ntime + com.nrgene + com.nkappa + com.npi;
for (i = 0; i < com.nOmega - 1; i++)
......@@ -1472,13 +1473,13 @@ int getnrate(int firsttime)
if (com.aaDist) com.nrate++;
if (!com.fix_omega) com.nrate++;
if (com.codonf) {
com.codonf = 0; puts("CodonFreq=0 reset for model=6.");
com.codonf = 0; puts("CodonFreq=0 reset for model=5.");
}
break;
case (FromCodon):
com.nrate = com.nkappa = (com.hkyREV ? 5 : !com.fix_kappa);
if (com.aaDist) com.nrate++;
if(com.fix_omega) error2("fix_omega = 1? omega is inestimable!");
if (com.fix_omega) error2("fix_omega = 1? omega is not estimable!");
com.omega = -1;
if (com.codonf) {
com.codonf = 0; puts("CodonFreq=0 reset for model=6.");
......@@ -1531,20 +1532,22 @@ int getnrate(int firsttime)
/* pi_T, pi_C, pi_A are counted as frequency parameters pi. */
if (com.codonf == 0)
com.npi = 0;
if (com.codonf == FMutSel0) /* FMutSel0: pi_TCA plus 20 AA freqs. */
com.npi = 3 + (com.npi ? 20 - 1 : 0);
else if (com.codonf == FMutSel) /* FMutSel: pi_TCA plus 60 codon freqs. */
com.npi = 3 + (com.npi ? com.ncode - 1 : 0);
else if (com.npi) {
if (com.codonf == F1x4 || com.codonf == F1x4MG) com.npi = 3;
else if (com.codonf == F3x4 || com.codonf == F3x4MG) com.npi = 9;
else if (com.codonf == Fcodon) com.npi = com.ncode - 1;
if (firsttime) {
if (com.codonf == FMutSel0) /* FMutSel0: pi_TCA plus 20 AA freqs. */
com.npi = 3 + (com.npi ? 20 - 1 : 0);
else if (com.codonf == FMutSel) /* FMutSel: pi_TCA plus 60 codon freqs. */
com.npi = 3 + (com.npi ? com.ncode - 1 : 0);
else if (com.npi) {
if (com.codonf == F1x4 || com.codonf == F1x4MG) com.npi = 3;
else if (com.codonf == F3x4 || com.codonf == F3x4MG) com.npi = 9;
else if (com.codonf == Fcodon) com.npi = com.ncode - 1;
}
}
com.nrate += com.npi;
if (com.aaDist != AAClasses) {
if (com.fix_kappa > 1) error2("fix_kappa>1, not tested."); /** ???? */
if (com.model > 0 && com.model!=FromCodon0 && (com.alpha || !com.fix_alpha))
if (com.model > 0 && com.model != FromCodon0 && (com.alpha || !com.fix_alpha))
error2("dN/dS ratios among branches not implemented for gamma");
if (com.model > 0 && com.clock)
error2("model and clock don't work together");
......@@ -1576,8 +1579,8 @@ int getnrate(int firsttime)
puts("\aGamma codon model: are you sure this is the model you want to use? ");
if (com.aaDist == 0) {
if ((!com.fix_omega || (com.Mgene && com.Mgene >= 3))
&& !com.NSsites && (com.model==0 || com.model == FromCodon0))
if ((!com.fix_omega || (com.Mgene && com.Mgene >= 3))
&& !com.NSsites && (com.model == 0 || com.model == FromCodon0))
com.nrate++;
}
else {
......@@ -1590,7 +1593,7 @@ int getnrate(int firsttime)
}
if (com.NSsites) {
if (com.NSsites == NSfreqs && com.ncatG != 5) {
if (com.NSsites == NSfreqs && com.ncatG != 5) {
puts("\nncatG changed to 5."); com.ncatG = 5;
}
if (com.model && com.NSsites)
......@@ -1602,7 +1605,7 @@ int getnrate(int firsttime)
case (NSpselection):
case (NSM2aRel):
com.ncatG = 3; break;
case (NSbetaw):
case (NSbetaw):
case (NS02normal):
if (firsttime) com.ncatG++; break;
}
......@@ -1846,7 +1849,7 @@ int SetxBound(int np, double xb[][2])
double alphab[] = { 0.02,49 }, betab[] = { 0.005,99 };
double rhob[] = { 0.01,0.99 }, pb[] = { .00001,.99999 };
double wb[] = { 0.0001,999 };
double w0b[] = { 0.0001, 0.5 }, w1b[] = { 0.5, 1.5 }; /* clade model D */
double w0b[] = { 0.0001, 1.0 }, w1b[] = { 0.01, 1.5 }; /* clade model D */
SetxBoundTimes(xb);
for (i = com.ntime; i < np; i++) for (j = 0; j < 2; j++) xb[i][j] = rateb[j];
......@@ -2067,7 +2070,7 @@ int GetInitialsCodon(double x[])
x[k++] = com.pf3x4[j * 4 + i] / (com.pf3x4[j * 4 + 3] + .02*rndu());
}
if (com.NSsites == 0 && (com.model == 0 || com.model == FromCodon0)) {
if (!com.aaDist) {
if (!com.aaDist) {
if (!com.fix_omega) x[k++] = com.omega;
}
else if (com.aaDist == AAClasses)
......@@ -2101,7 +2104,7 @@ int GetInitialsCodon(double x[])
}
}
if (com.model>=NSbranchB && com.model <= NSbranch3) { /* branch models */
if (com.model >= NSbranchB && com.model <= NSbranch3) { /* branch models */
if (com.model == NSbranchB) {
com.nbtype = tree.nbranch;
for (i = 0; i < tree.nbranch; i++)
......@@ -2128,7 +2131,7 @@ int GetInitialsCodon(double x[])
/* branch-site and clade models
com.nOmega=2 different w's at a site (three w's in the model: w0,w1,w2) */
if (com.model>=NSbranch2 && com.model <= NSbranch3 && com.NSsites) {
if (com.model >= NSbranch2 && com.model <= NSbranch3 && com.NSsites) {
if (com.model == NSbranch2) { /* branch-site models A & B */
com.ncatG = 4; K = 3;
if (com.NSsites == NSdiscrete)
......@@ -2730,7 +2733,7 @@ int SetParameters(double x[])
k = com.ntime + com.nrgene + com.nkappa + com.npi;
if (com.nrate) {
if ((com.model==0 || com.model == FromCodon0) && !com.aaDist && !com.fix_omega && !com.NSsites)
if ((com.model == 0 || com.model == FromCodon0) && !com.aaDist && !com.fix_omega && !com.NSsites)
com.omega = x[k];
if (com.seqtype == AAseq)
eigenQaa(NULL, Root, U, V, x + com.ntime + com.nrgene);
......@@ -2756,7 +2759,7 @@ int SetParameters(double x[])
*/
/* branch models */
if (com.seqtype == CODONseq && com.model>=NSbranchB && com.model <= NSbranch3 && com.NSsites == 0 && com.aaDist == 0) {
if (com.seqtype == CODONseq && com.model >= NSbranchB && com.model <= NSbranch3 && com.NSsites == 0 && com.aaDist == 0) {
for (j = 0; j < tree.nnode; j++) {
if (j == tree.root) continue;
if (com.fix_omega && (int)nodes[j].label == com.nOmega - 1)
......@@ -3199,7 +3202,7 @@ int eigenQcodon(int mode, double blength, double *S, double *dS, double *dN,
NEigenQ++;
if (blength >= 0 && (S == NULL || dS == NULL || dN == NULL)) error2("eigenQcodon");
memset(Q, 0, n*n*sizeof(double));
memset(Q, 0, n*n * sizeof(double));
for (i = 1; i < n; i++) {
ic1 = FROM61[i]; from[0] = ic1 / 16; from[1] = (ic1 / 4) % 4; from[2] = ic1 % 4;
for (j = 0; j < i; j++) {
......@@ -3415,16 +3418,16 @@ int eigenQaa(FILE *fout, double Root[], double U[], double V[], double rate[])
int Qcodon2aa(double Qc[], double pic[], double Qaa[], double piaa[])
{
/* Qc -> Qaa
This routine constructs the rate matrix for amino acid replacement from
the rate matrix for codon substitution, by congregating states in the
Markov chain. Both processes are time reversible, and only the
symmetrical part of the rate matrix are constructed. Codon frequencies
pic[] are used. They are constructed by assigning equal frequencies for
synonymous codons in the routine AA2Codonf().
Qaa(a_i, a_j) = sum_i sum_j (piC[i]*piC[j]]*Qc[i][j]) / (piAA[i]*piAA[j])
*/
/* Qc -> Qaa
This routine constructs the rate matrix for amino acid replacement from
the rate matrix for codon substitution, by congregating states in the
Markov chain. Both processes are time reversible, and only the
symmetrical part of the rate matrix are constructed. Codon frequencies
pic[] are used. They are constructed by assigning equal frequencies for
synonymous codons in the routine AA2Codonf().
Qaa(a_i, a_j) = sum_i sum_j (piC[i]*piC[j]]*Qc[i][j]) / (piAA[i]*piAA[j])
*/
int i, j, aai, aaj, nc = Nsensecodon, naa = 20;
double ti, tij;
......@@ -3851,7 +3854,7 @@ int AA2Codonf(double faa[20], double fcodon[])
*/
int ic, iaa, i, NCsyn[20];
for (iaa = 0; iaa < 20; iaa++)
for (iaa = 0; iaa < 20; iaa++)
NCsyn[iaa] = 0;
for (ic = 0; ic < 64; ic++) {
iaa = GeneticCode[com.icode][ic];
......@@ -5003,9 +5006,9 @@ int PairwiseAA(FILE *fout, FILE*f2AA)
char GetAASiteSpecies(int species, int sitepatt)
{
/* this returns the amino acid encoded by the codon at sitepatt in species.
Returns '*' if more than two amino acids or '-' if codon is --- or ***.
*/
/* this returns the amino acid encoded by the codon at sitepatt in species.
Returns '*' if more than two amino acids or '-' if codon is --- or ***.
*/
int n = com.ncode, c, naa, k;
char aa = 0, newaa = 0;
......@@ -5157,7 +5160,6 @@ int lfunNSsites_rate(FILE* frst, double x[], int np)
char sitelabel[96], *colors[5] = { "darkblue", "lightblue", "purple", "pinkred", "red" };
char *colorvalues[5] = { "[2,2,120]", "[133,57,240]", "[186,60,200]", "[200,60,160]", "[250,5,5]" };
if (com.nparK) error2("lfunNSsites_rate to be done for HMM.");
if ((meanw = (double*)malloc(com.npatt * sizeof(double))) == NULL)
......@@ -6302,10 +6304,10 @@ int lfunNSsites_M2M8(FILE* frst, double x[], int np)
printf("\nBEBing (dim = %d). This may take several minutes.", dim);
if (com.NSsites == 8) {
if (com.NSsites == 8) {
paras[0] = "p0"; paras[1] = "p"; paras[2] = "q"; paras[3] = "ws";
}
else if (!M2a) {
else if (!M2a) {
paras[0] = "p0"; paras[1] = "p1"; paras[2] = "w2";
}
else {
......@@ -6722,8 +6724,8 @@ int lfunNSsites_ACD(FILE* frst, double x[], int np)
8 March 2016, added BEB for clade model D. deleted ncatG=2 for D so that ncatG=3 for A C D.
*/
int n1d = 10, debug = 0, site = 0;
double w0b[] = { 0.0, 0.5 }; /* bounds for w0 */
double w1b[] = { 0.5, 1.5 }; /* bounds for w1 (for model D only) */
double w0b[] = { 0.0, 1.0 }; /* bounds for w0 */
double w1b[] = { 0.01, 1.5 }; /* bounds for w1 (for model D only) */
double wsb[] = { 1, 11 }; /* bounds for w2 in model A */
double w2b[] = { 0, 3 }; /* bounds for w2 w3 w4... in clade models C & D */
enum ModelsACD { mA, mC, mD };
......
src/evolver.c
View file @ bc993320
/* evolver.c
Copyright, Ziheng Yang, April 1995.
cl -Ot -O2 evolver.c tools.c
cl -Ot -O2 -DCodonNSbranches -FeevolverNSbranches.exe evolver.c tools.c
cl -Ot -O2 -DCodonNSsites -FeevolverNSsites.exe evolver.c tools.c
cl -Ot -O2 -DCodonNSbranchsites -FeevolverNSbranchsites.exe evolver.c tools.c
cl -Ox evolver.c tools.c
cl -Ox -DCodonNSbranches -FeevolverNSbranches.exe evolver.c tools.c
cl -Ox -DCodonNSsites -FeevolverNSsites.exe evolver.c tools.c
cl -Ox -DCodonNSbranchsites -FeevolverNSbranchsites.exe evolver.c tools.c
cc -fast -o evolver evolver.c tools.c -lm
cc -O3 -o evolver evolver.c tools.c -lm
cc -O4 -DCodonNSbranches -o evolverNSbranches evolver.c tools.c -lm
cc -O4 -DCodonNSsites -o evolverNSsites evolver.c tools.c -lm
cc -O4 -DCodonNSbranchsites -o evolverNSbranchsites evolver.c tools.c -lm
......@@ -36,14 +36,14 @@
struct CommonInfo {
unsigned char *z[2 * NS - 1];
char spname[NS][LSPNAME + 1], daafile[512], cleandata, readpattern;
unsigned char *spname[NS], daafile[512], cleandata, readpattern;
int ns, ls, npatt, np, ntime, ncode, clock, rooted, model, icode;
int seqtype, *pose, ncatG, NSsites;
int ngene, lgene[1], posG[1 + 1]; /* not used */
double piG[1][4], rgene[1]; /* not used */
double *fpatt, kappa, omega, alpha, pi[64], *conP, daa[20 * 20];
double freqK[NCATG], rK[NCATG];
char *siteID; /* used if ncatG>1 */
char *siteID; /* used if ncatG>1 */
double *siterates; /* rates for gamma or omega for site or branch-site models */
double *omegaBS, *QfactorBS; /* omega IDs for branch-site models */
} com;
......@@ -52,8 +52,8 @@ struct TREEB {
} tree;
struct TREEN {
int father, nson, sons[MAXNSONS], ibranch;
double branch, age, omega, label, *conP;
char *nodeStr, fossil;
double branch, age, omega, label, label2, *conP;
char *annotation, fossil;
} *nodes;
extern char BASEs[];
......@@ -89,7 +89,7 @@ double PMat[NCODE*NCODE], U[NCODE*NCODE], V[NCODE*NCODE], Root[NCODE];
static double Qfactor = -1, Qrates[5]; /* Qrates[] hold kappa's for nucleotides */
int main(int argc, char*argv[])
int main(int argc, char *argv[])
{
char *MCctlf = NULL, outf[512] = "evolver.out", treefile[512] = "mcmc.txt", mastertreefile[512] = "\0";
int i, option = -1, ntree = 1, rooted, BD = 0, gotoption = 0, pick1tree = 0;
......@@ -107,7 +107,7 @@ int main(int argc, char*argv[])
if (argc == 1)
printf("Results for options 1-4 & 8 go into %s\n", outf);
else if (option != 5 && option != 6 && option != 7 && option != 9) {
puts("Usage: \n\tevolver \n\tevolver option# MyDataFile"); exit(-1);
puts("Usage: \n\tevolver \n\tevolver option datafile"); exit(-1);
}
if (option >= 4 && option <= 6)
MCctlf = argv[2];
......@@ -117,6 +117,11 @@ int main(int argc, char*argv[])
if (argc > 4) sscanf(argv[4], "%d", &pick1tree);
}
for (i = 0; i < NS; i++) {
com.spname[i] = (char*)malloc(LSPNAME * sizeof(char*));
if (com.spname[i] == NULL) error2("oom");
}
#if defined (CodonNSbranches)
option = 6; com.model = 1;
MCctlf = (argc == 3 ? argv[2] : "MCcodonNSbranches.dat");
......@@ -219,8 +224,9 @@ int main(int argc, char*argv[])
Simulate(MCctlf ? MCctlf : MCctlf0[option - 5]);
}
else if (option == 9) {
CladeSupport(fout, treefile, 1, mastertreefile, pick1tree);
/* CladeMrBayesProbabilities("/papers/BPPJC3sB/Karol.trees"); */
/* CladeMrBayesProbabilities("\data\Lakner2008SB/M1510YY03.nex.t"); */
}
return(0);
}
......@@ -248,11 +254,11 @@ int between_f_and_x(void)
void LabelClades(FILE *fout)
{
/* This reads in a tree and scan species names to check whether they form a
paraphyletic group and then label the clade.
It assumes that the tree is unrooted, and so goes through two rounds to check
whether the remaining seqs form a monophyletic clade.
*/
/* This reads in a tree and scan species names to check whether they form a
paraphyletic group and then label the clade.
It assumes that the tree is unrooted, and so goes through two rounds to check
whether the remaining seqs form a monophyletic clade.
*/
FILE *ftree;
int unrooted = 1, iclade, sizeclade, mrca, paraphyl, is, imrca, i, j, k, lasts, haslength;
char key[96] = "A", treef[64] = "/A/F/flu/HA.all.prankcodon.tre", *p, chosen[NS], *endstr = "end";
......@@ -271,12 +277,12 @@ void LabelClades(FILE *fout)
i = (com.ns * 2 - 1) * sizeof(struct TREEN);
if ((nodes = (struct TREEN*)malloc(i)) == NULL) error2("oom");
for (i = 0; i < com.ns * 2 - 1; i++) nodes[i].nodeStr = NULL;
for (i = 0; i < com.ns * 2 - 1; i++) nodes[i].annotation = NULL;
for (i = 0; i < com.ns - 1; i++) {
anc[i] = (int*)malloc((com.ns / SI + 1) * sizeof(int));
if (anc[i] == NULL) error2("oom");
}
ReadTreeN(ftree, &haslength, &j, 1, 0);
ReadTreeN(ftree, &haslength, 1, 0);
fclose(ftree);
if (debug) { OutTreeN(F0, 1, PrNodeNum); FPN(F0); }
......@@ -287,8 +293,7 @@ void LabelClades(FILE *fout)
break;
for (i = 0; i < com.ns; i++)
chosen[i] = '\0';
k = strlen(key);
for (i = 0; i < com.ns; i++) {
if ((p = strstr(com.spname[i], key))
......@@ -400,7 +405,7 @@ void TreeDistanceDistribution(FILE* fout)
puts("\ntree #: mean prop of tree pairs with 0 1 2 ... shared bipartitions\n");
fputs("\ntree #: prop of tree pairs with 0 1 2 ... shared bipartitions\n", fout);
for (i = 0; i < ntree; i++) {
ReadTreeN(ftree, &j, &k, 0, 1);
ReadTreeN(ftree, &j, 0, 1);
nib[i] = tree.nbranch - com.ns;
Tree2Partition(partition + i*lpart);
}
......@@ -462,7 +467,7 @@ void TreeDistances(FILE* fout)
if (k == 2) { /* pairwise comparisons */
fputs("Number of identical bi-partitions in pairwise comparisons\n", fout);
for (i = 0; i < ntree; i++) {
ReadTreeN(ftree, &j, &k, 0, 1);
ReadTreeN(ftree, &j, 0, 1);
nib[i] = tree.nbranch - com.ns;
Tree2Partition(partition + i*lpart);
}
......@@ -477,7 +482,7 @@ void TreeDistances(FILE* fout)
}
}
else { /* first vs. others */
ReadTreeN(ftree, &j, &k, 0, 1);
ReadTreeN(ftree, &j, 0, 1);
nib[0] = tree.nbranch - com.ns;
if (nib[0] == 0) error2("1st tree is a star tree..");
Tree2Partition(partition);
......@@ -496,7 +501,7 @@ void TreeDistances(FILE* fout)
fputs("\nCorrect internal branches compared with the 1st tree:\n", fout);
FOR(k, nib[0]) nIBsame[k] = 0;
for (i = 1, mp = vp = 0; i < ntree; i++, FPN(fout)) {
ReadTreeN(ftree, &j, &k, 0, 1);
ReadTreeN(ftree, &j, 0, 1);
nib[1] = tree.nbranch - com.ns;
Tree2Partition(partition + lpart);
nsame = NSameBranch(partition, partition + lpart, nib[0], nib[1], IBsame);
......@@ -823,7 +828,7 @@ void Simulate(char *ctlf)
if (fixtree) {
fscanf(fin, "%lf", &tlength); fgets(line, lline, fin);
if (ReadTreeN(fin, &i, &j, 1, 1)) /* might overwrite spname */
if (ReadTreeN(fin, &i, 1, 1)) /* might overwrite spname */
error2("err tree..");
if (i == 0) error2("use : to specify branch lengths in tree");
......@@ -897,9 +902,9 @@ void Simulate(char *ctlf)
for (i = 0; i < tree.nnode; i++)
blengthBS[i] = nodes[i].branch;
for (k = 0; k < com.ncatG; k++) {
ReadTreeN(fin, &i, &j, 0, 1);
ReadTreeN(fin, &i, 0, 1);
if (i) error2("do not include branch lengths except in the first tree.");
if (!j) error2("Use # to specify omega's for branches");
/* if (!j) error2("Use # to specify omega's for branches"); */
for (i = 0; i < tree.nnode; i++) com.omegaBS[i*com.ncatG + k] = nodes[i].label;
}
for (i = 0; i < tree.nnode; i++)
......@@ -907,7 +912,7 @@ void Simulate(char *ctlf)
nodes[i].branch = blengthBS[i]; nodes[i].label = nodes[i].omega = 0;
}
for (i = 0; i < tree.nnode; i++) { /* print out omega as node labels. */
nodes[i].nodeStr = pc = (char*)malloc(20 * com.ncatG * sizeof(char));
nodes[i].annotation = pc = (char*)malloc(20 * com.ncatG * sizeof(char));
sprintf(pc, "'[%.2f", com.omegaBS[i*com.ncatG + 0]);
for (k = 1, pc += strlen(pc); k < com.ncatG; k++, pc += strlen(pc))
sprintf(pc, ", %.2f", com.omegaBS[i*com.ncatG + k]);
......@@ -1100,8 +1105,7 @@ void Simulate(char *ctlf)
PatternWeightJC69like();
}
if (format == 2 || format == 3) fprintf(fseq, "\n\n[Replicate # %d]\n", ir + 1);
printSeqs(fseq, NULL, NULL, format); /* printsma not usable as it codes into 0,1,...,60. */
printSeqs(fseq, com.z, com.spname, com.ns, com.ls, com.npatt, com.fpatt, NULL, NULL, format); /* printsma not usable as it codes into 0,1,...,60. */
if ((format == 2 || format == 3) && !fixtree) {
fprintf(fseq, "\nbegin tree;\n tree true_tree = [&U] ");
OutTreeN(fseq, 1, 1); fputs(";\n", fseq);
......@@ -1169,7 +1173,7 @@ void Simulate(char *ctlf)
for (j = 0; j < com.ns * 2 - 1; j++) free(com.z[j]);
free(space);
if (com.model && com.NSsites) /* branch-site model */
for (i = 0; i < tree.nnode; i++) free(nodes[i].nodeStr);
for (i = 0; i < tree.nnode; i++) free(nodes[i].annotation);
free(nodes);
if (com.alpha || com.ncatG) {
free(com.siterates); com.siterates = NULL;
......@@ -1218,10 +1222,10 @@ int GetSpnamesFromMB(FILE *fmb, char line[], int lline)
return(0);
}
char *GrepLine(FILE*fin, char*query, char* line, int lline)
char *GrepLine(FILE *fin, char *query, char *line, int lline)
{
/* This greps infile to search for query[], and returns NULL or line[].
*/
/* This greps infile to search for query[], and returns NULL or line[].
*/
char *p = NULL;
rewind(fin);
......@@ -1235,17 +1239,16 @@ char *GrepLine(FILE*fin, char*query, char* line, int lline)
void CladeMrBayesProbabilities(char treefile[])
{
/* This reads a tree from treefile and then scans a set of MrBayes output files
(mbfiles) to retrieve posterior probabilities for every clade in that tree.
It first scans the first mb output file to get the species names.
/* This reads a tree from treefile and then scans a set of MrBayes output files
(mbfiles) to retrieve posterior probabilities for every clade in that tree.
It first scans the first mb output file to get the species names.
Sample mb output:
6 -- ...........................************* 8001 1.000 0.005 (0.000)
7 -- ....................******************** 8001 1.000 0.006 (0.000)
Sample mb output:
6 -- ...........................************* 8001 1.000 0.005 (0.000)
7 -- ....................******************** 8001 1.000 0.006 (0.000)
Note 4 Jan 2014: This uses parti2B[], and is broken after i rewrote
Tree2Partition().
*/
Note 4 Jan 2014: This uses parti2B[], and is broken after i rewrote Tree2Partition().
*/
int lline = 100000, i, j, k, nib, inode, parti2B[NS];
char line[100000], *partition, *p;
char symbol[2] = ".*", cladestr[NS + 1] = { 0 };
......@@ -1266,13 +1269,13 @@ void CladeMrBayesProbabilities(char treefile[])
if (i && i != com.ns) error2("do you mean to specify ns in the tree file?");
i = (com.ns * 2 - 1) * sizeof(struct TREEN);
if ((nodes = (struct TREEN*)malloc(i)) == NULL) error2("oom");
ReadTreeN(ftree, &i, &j, 0, 1);
ReadTreeN(ftree, &i, 0, 1);
FPN(F0); OutTreeN(F0, 0, 0); FPN(F0); FPN(F0);
nib = tree.nbranch - com.ns;
for (i = 0; i < tree.nnode; i++) {
nodes[i].nodeStr = NULL;
if (i > com.ns) nodes[i].nodeStr = (char*)malloc(100 * sizeof(char));
nodes[i].annotation = NULL;
if (i > com.ns) nodes[i].annotation = (char*)malloc(100 * sizeof(char));
}
partition = (char*)malloc(nib*com.ns * sizeof(char));
......@@ -1294,19 +1297,19 @@ void CladeMrBayesProbabilities(char treefile[])
for (k = 0; k < nmbfiles; k++) {
if (GrepLine(fmb[k], cladestr, line, lline)) {
p = strstr(line, cladestr);
sscanf(p + com.ns, "%lf%lf\0", &t, &Pclade[i*nmbfiles + k]);
sscanf(p + com.ns, "%lf%lf", &t, &Pclade[i*nmbfiles + k]);
}
}
for (k = 0; k < nmbfiles; k++) printf("%6.2f", Pclade[i*nmbfiles + k]);
FPN(F0);
for (k = 0, p = nodes[inode].nodeStr; k < nmbfiles; k++) {
for (k = 0, p = nodes[inode].annotation; k < nmbfiles; k++) {
sprintf(p, "%3.0f%s", Pclade[i*nmbfiles + k] * 100, (k < nmbfiles - 1 ? "/" : ""));
p += 4;
}
}
FPN(F0); OutTreeN(F0, 1, PrLabel); FPN(F0);
for (i = 0; i < tree.nnode; i++) free(nodes[i].nodeStr);
for (i = 0; i < tree.nnode; i++) free(nodes[i].annotation);
free(nodes); free(partition); free(Pclade);
fclose(ftree);
for (k = 0; k < nmbfiles; k++) fclose(fmb[k]);
......
src/mcmctree.c
View file @ bc993320
This diff is collapsed.
src/paml.h
View file @ bc993320
......@@ -13,6 +13,7 @@
#include <limits.h>
#include <float.h>
#include <time.h>
#include <assert.h>
#define square(a) ((a)*(a))
#define FOR(i,n) for(i=0; i<n; i++)
......@@ -138,7 +139,7 @@ int f_mono_di (FILE *fout, char z[], int ls, int iring, double fb1[], double fb2
int PickExtreme (FILE *fout, char z[], int ls, int iring, int lfrag, int ffrag[]);
int print1seq (FILE*fout, unsigned char *z, int ls, int pose[]);
void printSeqs(FILE *fout, int *pose, char keep[], int format);
void printSeqs(FILE *fout, unsigned char *z[], unsigned char *spnames[], int ns, int ls, int npatt, double fpatt[], int *pose, char keep[], int format);
int printPatterns(FILE *fout);
void printSeqsMgenes (void);
int printsma (FILE*fout, char*spname[], unsigned char*z[], int ns, int l, int lline, int gap, int seqtype,
......@@ -150,10 +151,8 @@ double SeqDivergence (double x[], int model, double alpha, double *kapa);
int symtest (double freq[], int n, int nobs, double space[], double *chisym,
double* chihom);
int dSdNNG1986(char *z1, char *z2, int lc, int icode, int transfed, double *dS, double *dN, double *Ssites, double *Nsites);
int difcodonNG(char codon1[], char codon2[], double *SynSite,double *AsynSite,
double *SynDif, double *AsynDif, int transfed, int icode);
int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3],
double vdiff[3], int transfed, int icode);
int difcodonNG(char codon1[], char codon2[], double *SynSite,double *AsynSite, double *SynDif, double *AsynDif, int transfed, int icode);
int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3], double vdiff[3], int transfed, int icode);
int testTransP (double P[], int n);
double testDetailedBalance (double P[], double pi[], int n);
void pijJC69 (double pij[2], double t);
......@@ -163,12 +162,10 @@ int PMatTN93 (double P[], double a1t, double a2t, double bt, double pi[]);
int PMatUVRoot (double P[],double t,int n,double U[],double V[],double Root[]);
int PMatCijk (double PMat[], double t);
int PMatQRev(double P[], double pi[], double t, int n, double space[]);
int EvolveHKY85 (char source[], char target[], int ls, double t,
double rates[], double pi[], double kapa, int isHKY85);
int EvolveHKY85 (char source[], char target[], int ls, double t, double rates[], double pi[], double kapa, int isHKY85);
double DistanceIJ (int is, int js, int model, double alpha, double *kappa);
int DistanceMatNuc (FILE *fout, FILE*f2base, int model, double alpha);
int EigenQGTRbase (FILE* fout, double kappa[], double pi[],
int *nR, double Root[], double Cijk[]);
int EigenQGTRbase (FILE* fout, double kappa[], double pi[], int *nR, double Root[], double Cijk[]);
int DistanceMatNG86 (FILE *fout, FILE*fds, FILE*fdn, FILE*dt, double alpha);
int setmark_61_64 (void);
......@@ -236,12 +233,11 @@ int eigenRealSym(double A[], int n, double Root[], double work[]);
int MeanVar (double x[], int n, double *mean, double *var);
int variance (double x[], int n, int nx, double mx[], double vx[]);
int correl (double x[], double y[], int n, double *mx, double *my,
double *vxx, double *vxy, double *vyy, double *r);
int correl (double x[], double y[], int n, double *mx, double *my, double *vxx, double *vxy, double *vyy, double *r);
int comparefloat (const void *a, const void *b);
int comparedouble (const void *a, const void *b);
double Eff_IntegratedCorrelationTime (double x[], int n, double *mx, double *varx, double *rho1);
double Eff_IntegratedCorrelationTime2(double x[], int n, int nbatch, double *mx, double *varx);
double Eff_IntegratedCorrelationTime (double x[], int n, double *mx, double *vx, double *rho1);
double Eff_IntegratedCorrelationTime2(double x[], int n, int nbatch, double *mx, double *vx);
int HPDinterval(double x[], int n, double HPD[2], double alpha);
int DescriptiveStatistics(FILE *fout, char infile[], int nbin, int propternary, int SkipColumns);
......@@ -249,12 +245,9 @@ int DescriptiveStatisticsSimple (FILE *fout, char infile[], int SkipColumns);
int splitline(char line[], int nfields, int fields[]);
int scanfile (FILE*fin, int *nrecords, int *nx, int *HasHeader, char line[], int ifields[]);
double bound (int nx, double x0[], double p[], double x[],
int (*testx) (double x[], int nx));
int gradient (int n, double x[], double f0, double g[],
double (* fun)(double x[],int n), double space[], int Central);
int Hessian (int nx, double x[], double f, double g[], double H[],
double (*fun)(double x[], int n), double space[]);
double bound (int nx, double x0[], double p[], double x[], int (*testx) (double x[], int nx));
int gradient (int n, double x[], double f0, double g[], double (* fun)(double x[],int n), double space[], int Central);
int Hessian (int nx, double x[], double f, double g[], double H[], double (*fun)(double x[], int n), double space[]);
int HessianSKT2004 (double xmle[], double lnLm, double g[], double H[]);
int H_end (double x0[], double x1[], double f0, double f1, double e1, double e2, int n);
......@@ -264,7 +257,6 @@ double LineSearch2 (double(*fun)(double x[],int n), double *f, double x0[],
void xtoFreq(double x[], double freq[], int n);
int SetxBound (int np, double xb[][2]);
int ming1 (FILE *fout, double *f, double (* fun)(double x[], int n),
int (*dfun) (double x[], double *f, double dx[], int n),
......@@ -288,19 +280,19 @@ int nls2 (FILE *fout, double *sx, double * x0, int nx,
int (*testx) (double x[], int nx),
int ny, double e);
int ResetFinetuneSteps(FILE *fout, double Pjump[], double finetune[], int nsteps);
int ResetStepLengths(FILE *fout, double Pjump[], double finetune[], int nsteps);
/* tree structure functions in treesub.c */
void NodeToBranch (void);
void BranchToNode (void);
void ClearNode (int inode);
int ReadTreeN (FILE *ftree, int *haslength, int *haslabel, int copyname, int popline);
int ReadTreeN (FILE *ftree, int *haslength, int copyname, int popline);
int ReadTreeB (FILE *ftree, int popline);
int OutTreeN (FILE *fout, int spnames, int printopt);
int OutTreeB (FILE *fout);
int DeRoot (void);
int GetSubTreeN (int keep[], int space[]);
void printtree (int timebranches);
void PrintTree (int timebranches);
void PointconPnodes (void);
int SetBranch (double x[]);
int DistanceMat (FILE *fout, int ischeme, double alfa, double *kapa);
......@@ -317,8 +309,7 @@ int PopEmptyLines (FILE* fseq, int lline, char line[]);
int hasbase (char *str);
int blankline (char *str);
void BranchLengthBD(int rooted, double birth, double death, double sample,
double mut);
void BranchLengthBD(int rooted, double birth, double death, double sample, double mut);
int RandomLHistory (int rooted, double space[]);
void Tree2Partition (char partition[]);
......@@ -359,8 +350,7 @@ int GetSubSeqs(int nsnew);
/* functions for evolving sequences */
int GenerateSeq (void);
int Rates4Sites (double rates[],double alfa,int ncatG,int ls, int cdf,
double space[]);
int Rates4Sites (double rates[],double alfa,int ncatG,int ls, int cdf, double space[]);
void Evolve (int inode);
void EvolveJC (int inode);
......@@ -396,8 +386,6 @@ enum {PrBranch=1, PrNodeNum=2, PrLabel=4, PrNodeStr=8, PrAge=16, PrOmega=32} Out
#define PAML_RELEASE 0
#define FullSeqNames 0 /* 1: numbers at the beginning of sequence name are part of name */
#define pamlVerStr "paml version 4.9h, March 2018"
#define pamlVerStr "paml version 4.9i, September 2018"
#endif
src/pamp.c
View file @ bc993320
This diff is collapsed.
src/tools.c
View file @ bc993320
This diff is collapsed.