Skip to content

Commit

Permalink
remove compiler warning messages, #4
Browse files Browse the repository at this point in the history
  • Loading branch information
@ziheng-yang
ziheng-yang committed Sep 8, 2020
1 parent a88bd93 commit 881a9f1
Show file tree
Hide file tree
Showing 19 changed files with 3,572 additions and 3,187 deletions.
184 changes: 92 additions & 92 deletions Technical/Simulation/Codon/PositiveSites.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
/* PositiveSites.c
Ziheng Yang, June 2004
Ziheng Yang, June 2004
cc -O2 -DNEB -o PositiveSitesNEB PositiveSites.c -lm
cc -O2 -DBEB -o PositiveSitesBEB PositiveSites.c -lm
Expand All @@ -13,9 +13,9 @@
PositiveSitesBEB <#sites> <#repl> <Evolverf> <Codemlf>
PositiveSitesBS <#sites> <#repl> <Evolverf> <Codemlf> <positive site classes>
This compares siteID from evolverNSsites and mlc from codeml to calculate
the accuracy, power, and false positive rate of codeml inference of sites
under positive selection. The measures are defined as follows (Anisimova et
This compares siteID from evolverNSsites and mlc from codeml to calculate
the accuracy, power, and false positive rate of codeml inference of sites
under positive selection. The measures are defined as follows (Anisimova et
al. 2002; Wong et al. 2004).
codeml inference
Expand All @@ -28,13 +28,13 @@
Power = N++/N+.
FalsePositive = N-+/N-.
The program collects N++ (NmatchB & NmatchC), N+. (NEvolver), and
The program collects N++ (NmatchB & NmatchC), N+. (NEvolver), and
N.+ (NCodemlB & NcodemlC), and then calculates the three measures as above.
Note that codeml inference depends on cutoff P, hence the B (for binned) and
C (for cumulative) difference. All proportions are calculated as the ratio
of averages, taking the ratio after counting sites over replicate data sets.
Note that codeml inference depends on cutoff P, hence the B (for binned) and
C (for cumulative) difference. All proportions are calculated as the ratio
of averages, taking the ratio after counting sites over replicate data sets.
Output is on the screen.
*/

#include <stdio.h>
Expand All @@ -48,149 +48,149 @@
*/

#if (defined BEB)
int model=0, NSsites=3; char *Codemlf="mlc", startCodeml[]="Bayes Empirical", startEvolver[]="replicate ";
int model = 0, NSsites = 3; char* Codemlf = "mlc", startCodeml[] = "Bayes Empirical", startEvolver[] = "replicate ";
#elif(defined NEB)
int model=0, NSsites=3; char *Codemlf="mlc", startCodeml[]="Naive Empirical", startEvolver[]="replicate ";
int model = 0, NSsites = 3; char* Codemlf = "mlc", startCodeml[] = "Naive Empirical", startEvolver[] = "replicate ";
#elif(defined BranchSite)
int model=1, NSsites=3; char *Codemlf="mlc", startCodeml[]="Bayes Empirical", startEvolver[]="replicate ";
int model = 1, NSsites = 3; char* Codemlf = "mlc", startCodeml[] = "Bayes Empirical", startEvolver[] = "replicate ";
int nPositiveClass, PositiveClass[1000];
#endif


int main (int argc, char* argv[])
int main(int argc, char* argv[])
{
int nbin=21, noisy=0, nr=1, ls=100;
double PCut[]={.525, .55, .575, .6, .625, .65, .675, .7, .725, .75, .775, .8, .825, .85, .875, .9, .925, .95, .975, .99, 1}, PCut0=0.5;
/*
int nbin=11, noisy=0, nr=1, ls=100;
double PCut[]={0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.99, 1}, PCut0=0.5;
*/
int ir, ib, i,j,k, ch, nmatch0;
double NmatchB[50], NCodemlB[50], NEvolver=0;
int nbin = 21, noisy = 0, nr = 1, ls = 100;
double PCut[] = { .525, .55, .575, .6, .625, .65, .675, .7, .725, .75, .775, .8, .825, .85, .875, .9, .925, .95, .975, .99, 1 }, PCut0 = 0.5;
/*
int nbin=11, noisy=0, nr=1, ls=100;
double PCut[]={0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.99, 1}, PCut0=0.5;
*/
int ir, ib, i, j, k, ch, nmatch0;
double NmatchB[50], NCodemlB[50], NEvolver = 0;
double NmatchC[50], NCodemlC[50];
int *siteEvolver, *siteCodeml, *ibin, nsiteEvolver, nsiteCodeml, lline=1000;
int* siteEvolver, * siteCodeml, * ibin, nsiteEvolver, nsiteCodeml, lline = 1000;
int nPositiveClass, PositiveClass[100];
char *Evolverf="siterates.txt", line[1001];
char* Evolverf = "siterates.txt", line[1001];
double p, AccuracyB, AccuracyC, Power, FalsePositive;
FILE *fEvolver, *fCodeml;
FILE* fEvolver, * fCodeml;

if(model)
if (model)
puts("Usage:\n\tPositiveSitesBR <#sites> <#repl> <Evolverf> <Codemlf> <positive site classes>");
else {
puts("Usage:\n\tPositiveSitesBEB <#sites> <#repl>");
puts("Usage:\n\tPositiveSitesBEB <#sites> <#repl> <Evolverf> <Codemlf>");
}

if(argc<3) exit(-1);
if(argc>1) sscanf(argv[1],"%d", &ls);
if(argc>2) sscanf(argv[2],"%d", &nr);
if(argc>3) Evolverf=argv[3];
if(argc>4) Codemlf=argv[4];
printf("%d sites, %d replicate.\n", ls,nr);
if (argc < 3) exit(-1);
if (argc > 1) sscanf(argv[1], "%d", &ls);
if (argc > 2) sscanf(argv[2], "%d", &nr);
if (argc > 3) Evolverf = argv[3];
if (argc > 4) Codemlf = argv[4];
printf("%d sites, %d replicate.\n", ls, nr);
printf("codeml file %s starts with '%s' for each replicate.\n", Codemlf, startCodeml);

if(nr>20) { printf("Hit Enter to continue."); getchar(); }
if (nr > 20) { printf("Hit Enter to continue."); getchar(); }

fEvolver=fopen(Evolverf,"r"); fCodeml=fopen(Codemlf,"r");
if(!fEvolver || !fCodeml) { puts("file error"); exit(-1); }
fEvolver = fopen(Evolverf, "r"); fCodeml = fopen(Codemlf, "r");
if (!fEvolver || !fCodeml) { puts("file error"); exit(-1); }

siteEvolver=(int*)malloc(ls*3*sizeof(int));
if(siteEvolver==NULL) { puts("oom"); exit(-1); }
siteCodeml=siteEvolver+ls; ibin=siteCodeml+ls;
siteEvolver = (int*)malloc(ls * 3 * sizeof(int));
if (siteEvolver == NULL) { puts("oom"); exit(-1); }
siteCodeml = siteEvolver + ls; ibin = siteCodeml + ls;

for(i=0; i<nbin; i++) {
NmatchB[i]=NCodemlB[i]=0;
NmatchC[i]=NCodemlC[i]=0;
for (i = 0; i < nbin; i++) {
NmatchB[i] = NCodemlB[i] = 0;
NmatchC[i] = NCodemlC[i] = 0;
}

if(model && NSsites) { /* branch-site model */
nPositiveClass = argc-5;
for(i=0; i<nPositiveClass; i++)
sscanf(argv[5+i], "%d", &PositiveClass[i]);
if (model && NSsites) { /* branch-site model */
nPositiveClass = argc - 5;
for (i = 0; i < nPositiveClass; i++)
sscanf(argv[5 + i], "%d", &PositiveClass[i]);
printf("%d site classes are under positive selection: ", nPositiveClass);
for(i=0; i<nPositiveClass; i++) printf(" %2d", PositiveClass[i]);
for (i = 0; i < nPositiveClass; i++) printf(" %2d", PositiveClass[i]);
printf("\n");
}
for(ir=0; ir<nr; ir++) {
for (ir = 0; ir < nr; ir++) {
/* Read true sites from evovler siteID */
for( ; ; ) {
if(fgets(line, lline, fEvolver)==NULL) break;
if(strstr(line, startEvolver)) break;
for (; ; ) {
if (fgets(line, lline, fEvolver) == NULL) break;
if (strstr(line, startEvolver)) break;
}

if(NSsites && !model) { /* site models */
if(!strchr(line,':')) { puts("did not find ':' in line."); exit(-1); }
sscanf(strchr(line,':')+1, "%d", &nsiteEvolver);
if(nsiteEvolver>ls) { puts("Too many sites. ls wrong?"); exit(-1); }
for(i=0; i<nsiteEvolver; i++)
if (NSsites && !model) { /* site models */
if (!strchr(line, ':')) { puts("did not find ':' in line."); exit(-1); }
sscanf(strchr(line, ':') + 1, "%d", &nsiteEvolver);
if (nsiteEvolver > ls) { puts("Too many sites. ls wrong?"); exit(-1); }
for (i = 0; i < nsiteEvolver; i++)
fscanf(fEvolver, "%d", &siteEvolver[i]);
}
else { /* branch-site models */
for(i=0,nsiteEvolver=0; i<ls; i++) {
for (i = 0, nsiteEvolver = 0; i < ls; i++) {
fscanf(fEvolver, "%d", &k);
for(j=0; j<nPositiveClass; j++)
if(k==PositiveClass[j]) {
siteEvolver[nsiteEvolver++]=i+1;
for (j = 0; j < nPositiveClass; j++)
if (k == PositiveClass[j]) {
siteEvolver[nsiteEvolver++] = i + 1;
break;
}
}
}

NEvolver+=nsiteEvolver;
if(noisy) {
NEvolver += nsiteEvolver;
if (noisy) {
printf("\n\n%d sites from Evolver:\n", nsiteEvolver);
for(i=0;i<nsiteEvolver; i++) printf(" %3d", siteEvolver[i]);
for (i = 0; i < nsiteEvolver; i++) printf(" %3d", siteEvolver[i]);
}

/* Read inferred sites and probs from codeml mlc or rst, bin probs */
for( ; ; ) {
if(fgets(line, lline, fCodeml)==NULL) break;
if(strstr(line, startCodeml)) break;
for (; ; ) {
if (fgets(line, lline, fCodeml) == NULL) break;
if (strstr(line, startCodeml)) break;
}
for(i=0; i<3; i++) fgets(line, lline, fCodeml);
for(i=nsiteCodeml=0; i<ls; i++,nsiteCodeml++) {
if(fscanf(fCodeml, "%d %c%lf", &siteCodeml[i], &ch, &p)!=3) break;
for (i = 0; i < 3; i++) fgets(line, lline, fCodeml);
for (i = nsiteCodeml = 0; i < ls; i++, nsiteCodeml++) {
if (fscanf(fCodeml, "%d %d%lf", &siteCodeml[i], &ch, &p) != 3) break;
fgets(line, lline, fCodeml);
for(j=0; j<nbin-1; j++) if(p<=PCut[j]) break;
ibin[i]=j;
for (j = 0; j < nbin - 1; j++) if (p <= PCut[j]) break;
ibin[i] = j;
}
if(noisy) {
if (noisy) {
printf("\n%d sites from codeml at 50%%:\n", nsiteCodeml);
for(i=0; i<nsiteCodeml; i++) printf("%4d", siteCodeml[i]);
for (i = 0; i < nsiteCodeml; i++) printf("%4d", siteCodeml[i]);
}

/* count matches by going through codeml sites */
for(i=0,nmatch0=0; i<nsiteCodeml; i++) {
ib=ibin[i];
for (i = 0, nmatch0 = 0; i < nsiteCodeml; i++) {
ib = ibin[i];
NCodemlB[ib]++;
for(j=0; j<=ib; j++) NCodemlC[j]++;
for(j=0; j<nsiteEvolver; j++)
if(siteCodeml[i]==siteEvolver[j]) break;
if(j<nsiteEvolver) { /* a match */
for (j = 0; j <= ib; j++) NCodemlC[j]++;
for (j = 0; j < nsiteEvolver; j++)
if (siteCodeml[i] == siteEvolver[j]) break;
if (j < nsiteEvolver) { /* a match */
nmatch0++;
NmatchB[ib]++;
for(j=0; j<=ib; j++)
for (j = 0; j <= ib; j++)
NmatchC[j]++;
}
}

printf("\nReplicate %3d: %3d evolver sites, %3d codeml sites at 50%%, %3d matches",
ir+1,nsiteEvolver,nsiteCodeml,nmatch0);
printf("\nReplicate %3d: %3d evolver sites, %3d codeml sites at 50%%, %3d matches",
ir + 1, nsiteEvolver, nsiteCodeml, nmatch0);
}

printf("\n\n%6s%22s%10s%17s%10s%10s\n\n",
printf("\n\n%6s%22s%10s%17s%10s%10s\n\n",
"P", "AccuracyBin", "Pcut", "AccuracyCum", "Power", "FalsePos");
for(j=0; j<nbin; j++) {
AccuracyB = (NmatchB[j] ? NmatchB[j]/NCodemlB[j] : 0);
AccuracyC = (NmatchC[j] ? NmatchC[j]/NCodemlC[j] : 0);
Power = (NmatchC[j] ? NmatchC[j]/NEvolver : 0);
FalsePositive = NCodemlC[j]-NmatchC[j];
if(FalsePositive) FalsePositive/=(ls*nr-NEvolver);

p = (j==0 ? PCut0 : PCut[j-1]);
printf("%5.3f - %5.3f: %7.3f (%5.0f) ", p, PCut[j], AccuracyB, NCodemlB[j]);
printf( " >%4.3f: %7.3f (%5.0f) %7.3f %7.3f\n", p, AccuracyC, NCodemlC[j], Power, FalsePositive);
for (j = 0; j < nbin; j++) {
AccuracyB = (NmatchB[j] ? NmatchB[j] / NCodemlB[j] : 0);
AccuracyC = (NmatchC[j] ? NmatchC[j] / NCodemlC[j] : 0);
Power = (NmatchC[j] ? NmatchC[j] / NEvolver : 0);
FalsePositive = NCodemlC[j] - NmatchC[j];
if (FalsePositive) FalsePositive /= (ls * nr - NEvolver);

p = (j == 0 ? PCut0 : PCut[j - 1]);
printf("%5.3f - %5.3f: %7.3f (%5.0f) ", p, PCut[j], AccuracyB, NCodemlB[j]);
printf(" >%4.3f: %7.3f (%5.0f) %7.3f %7.3f\n", p, AccuracyC, NCodemlC[j], Power, FalsePositive);
}
printf( "\nTrue positive sites from evolver: %5.0f out of %5d total sites\n", NEvolver,ls*nr);
printf("\nTrue positive sites from evolver: %5.0f out of %5d total sites\n", NEvolver, ls * nr);
fclose(fEvolver); fclose(fCodeml);
}
6 changes: 3 additions & 3 deletions examples/TipDate.HIV2/mcmctree.ctl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
outfile = out

seqtype = 0 * 0: nucleotides; 1:codons; 2:AAs
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 * 3 for mt2G3P and 6 for 4G6P
clock = 1 * 1: global clock; 2: independent rates; 3: correlated rates
Expand Down Expand Up @@ -33,5 +33,5 @@

print = 1
burnin = 20000
sampfreq = 2
nsample = 100000
sampfreq = 10
nsample = 500000
1 change: 1 addition & 0 deletions src/BFdriver.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
using bpp and mcmctree.
cc -o BFdriver -O3 BFdriver.c tools.c -lm
cl -Ox BFdriver.c tools.c
BFdriver <controlfilename> <npoints> <scriptname.sh>
BFdriver mcmctree.ctl 16 tmp.sh
Expand Down
60 changes: 40 additions & 20 deletions src/Makefile
View file
Original file line number Diff line number Diff line change
@@ -1,27 +1,47 @@
PRGS = baseml codeml basemlg mcmctree pamp evolver yn00 chi2
PRGS = baseml codeml basemlg mcmctree pamp evolver yn00 chi2
CC = cc # cc, gcc, cl

CFLAGS = -O3
#CFLAGS = -fast
CFLAGS = -O3 -Wall -Wno-unused-result
#CC = icc
#CFLAGS = -fast -Wall

LIBS = -lm # -lM

all : $(PRGS)

baseml : baseml.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ baseml.c tools.c $(LIBS)
basemlg : basemlg.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ basemlg.c tools.c $(LIBS)
codeml : codeml.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ codeml.c tools.c $(LIBS)
evolver : evolver.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ evolver.c tools.c $(LIBS)
pamp : pamp.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ pamp.c tools.c $(LIBS)
mcmctree : mcmctree.c tools.c treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ mcmctree.c tools.c $(LIBS)
$(CC) $(CFLAGS) -o infinitesites -D INFINITESITES mcmctree.c tools.c $(LIBS)
yn00: yn00.c tools.c paml.h
$(CC) $(CFLAGS) -o $@ yn00.c tools.c $(LIBS)
chi2 : chi2.c
baseml : baseml.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ baseml.o tools.o $(LIBS)
basemlg : basemlg.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ basemlg.o tools.o $(LIBS)
codeml : codeml.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ codeml.o tools.o $(LIBS)
evolver : evolver.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ evolver.o tools.o $(LIBS)
pamp : pamp.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ pamp.o tools.o $(LIBS)
mcmctree : mcmctree.o tools.o treesub.c treespace.c paml.h
$(CC) $(CFLAGS) -o $@ mcmctree.c tools.o $(LIBS)
$(CC) $(CFLAGS) -o infinitesites -D INFINITESITES mcmctree.c tools.o $(LIBS)
yn00: yn00.o tools.o paml.h
$(CC) $(CFLAGS) -o $@ yn00.o tools.o $(LIBS)
chi2 : chi2.o
$(CC) $(CFLAGS) -o $@ chi2.c $(LIBS)

tools.o : paml.h tools.c
$(CC) $(CFLAGS) -c tools.c
baseml.o : paml.h baseml.c treesub.c treespace.c
$(CC) $(CFLAGS) -c baseml.c
basemlg.o : paml.h basemlg.c treesub.c
$(CC) $(CFLAGS) -c basemlg.c
codeml.o : paml.h codeml.c treesub.c treespace.c
$(CC) $(CFLAGS) -c codeml.c
evolver.o: evolver.c treesub.c treespace.c
$(CC) $(CFLAGS) -c evolver.c
mcmctree.o : paml.h mcmctree.c treesub.c treespace.c
$(CC) $(CFLAGS) -c mcmctree.c
pamp.o : paml.h pamp.c treesub.c treespace.c
$(CC) $(CFLAGS) -c pamp.c
yn00.o : paml.h yn00.c
$(CC) $(CFLAGS) -c yn00.c

clean :
-rm *.o $(PRGS)
Loading

0 comments on commit 881a9f1

Please sign in to comment.