-
Notifications
You must be signed in to change notification settings - Fork 0
/
make_mr_bayes.R
executable file
·42 lines (31 loc) · 1.24 KB
/
make_mr_bayes.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
library(ape)
library(stringr)
myseqs.fn <- "./output/aligned.fasta" # input file
myseqs.nex.fn <- "./output/tonsa.nex" # output in Nexus format
myseqs.mb <- "./output/tonsa_mb.nex" # output in Nexus with MrBayes block
myseqs <- read.dna(myseqs.fn,format="fasta",as.matrix=FALSE)
myseqs.names <- names(myseqs)
write.nexus.data(as.character(myseqs),myseqs.nex.fn,interleaved=FALSE,gap="-",missing="N")
myseqs.nex <- readLines(myseqs.nex.fn)
# We then write part of the MrBayes block that specifies the model.
#GTR + I + G, from figuerra MS
# set outgroup to Temo_stylifera
# help lset
# For example, "lset nst=6 rates=gamma" would set the model to a general model of DNA substitution (the GTR) with gamma-distributed rate variation across sites.
mbblock1 <- "
begin mrbayes;
set autoclose=yes nowarn=yes;
lset nst=2 rates=invgamma ngammacat=4;
outgroup Temora_stylifera;
"
mbblock3 <- "
mcmc ngen=500000 nruns=2 nchains=2 samplefreq=1000;
sump burninfrac=0.25;
sumt burninfrac=0.25;
quit;
end;
"
myseqs.nexus.withmb <- paste(paste(myseqs.nex,collapse="\n"),
mbblock1,mbblock3,sep="")
write(myseqs.nexus.withmb,file=myseqs.mb)
# We can now run MrBayes with the command mb with the output filename as the single argument.