PageRenderTime 39ms CodeModel.GetById 31ms app.highlight 6ms RepoModel.GetById 1ms app.codeStats 0ms

/cons_test/codon_funcs.R

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