PageRenderTime 14ms CodeModel.GetById 8ms app.highlight 3ms RepoModel.GetById 1ms app.codeStats 1ms

/cons_test/aa_funcs.R

http://github.com/sbotond/phylosim
R | 99 lines | 74 code | 15 blank | 10 comment | 8 complexity | d3dde47519c637bc6c1402d391d5c816 MD5 | raw file
 1dir<-"aln_aa";
 2system(paste("mkdir",dir));	
 3               
 4# construct substitution process object
 5wag<-WAG();
 6
 7construct_root_sequence<-function(len){
 8        seq<-AminoAcidSequence(length=len);
 9	# attach process
10	attachProcess(seq, wag);
11	return(seq)
12}
13
14# simulate nucleotide data
15simulate_aa<-function(phylo,seq,len,rep){
16		# set all site states to NA
17		clearStates(seq);
18                # sample rates from a discrete gamma model
19                plusGamma(seq,wag,shape=wag.true.gamma.shape,ncat=8);
20                # sample states
21                sampleStates(seq);
22                # construct simulation object
23                sim<-PhyloSim(
24                        phylo=phylo,
25                        root.seq=seq
26                );      
27                # run simulation
28                Simulate(sim,quiet=TRUE);
29                # save alignment
30                fname<-paste(dir,"/aasim_",len,"_",rep,".fas",sep="");
31                saveAlignment(sim,file=fname,skip.internal=TRUE);
32                return(fname);
33}
34
35# estimate nucleotide model parameters
36estimate_aa<-function(phylo, len,rep){
37	old.wd<-getwd();
38	dir.name<-paste(paste("aasim_",len,"_",rep,sep=""));
39	dir.create(dir.name);
40	setwd(dir.name);
41
42	system(paste("(seqret -osf phylipnon -outseq sim.phy ../",dir,"/",dir.name,".fas 2>&1) >/dev/null",sep=""));
43	system("codeml ../ctl/codeml_aa.ctl >/dev/null");
44	lines<-file.slurp("mlc");	
45
46	res.tree<-NA;
47	res.alpha<-NA;
48	
49	i<-1;
50	while(i < length(lines)){
51		if(length(grep(pattern="^tree length =",x=lines[[i]],perl=TRUE,value=FALSE)) > 0){
52			res.tree<-lines[[i+4]];
53			i<-i+4;
54		} 
55		if(length(grep(pattern="^alpha ",x=lines[[i]],perl=TRUE,value=FALSE)) > 0){
56			res.alpha<-as.numeric(strsplit(lines[[i]],"=",fixed=TRUE)[[1]][[3]]);
57			i<-i+1;
58		} else {
59			i<-i+1;
60		}
61
62	}
63
64	cat(res.tree,file="res.nwk");
65	res.tree<-reorder(read.tree("res.nwk"),"pruningwise");
66		
67	setwd(old.wd);
68	system(paste("rm -fr",dir.name));
69	tmp<-list();
70	tmp$len=len;
71	tmp$alpha=res.alpha;
72	tmp$tree=res.tree;
73	return(list(tmp));
74}
75
76plot_dframe<-function(d,p){
77        t<-reorder(p,"pruningwise");
78
79        for(col in names(d)){
80		pl<-my_qplot(d,col);
81
82                true<-NA;
83                if(col=="len"){
84                        next();
85                }
86                else if (col=="alpha"){
87                        true<-wag.true.gamma.shape;
88                }
89                else { 
90                        edge.nr<-as.numeric(strsplit(col,"_")[[1]][[2]]);
91                        true<-t$edge.length[edge.nr];
92                }
93
94		pl<-hline(pl,true);
95                print(pl);
96                Sys.sleep(10);
97        }
98}
99