https://github.com/wf8/ChromoSSE
1.输入文件的准备.
两个文件, 一个是进化树文件, nwk格式就可以:
(((A: 0.1, B: 0.2): 0.3, C: 0.15): 0.1, D: 0.4);
另一个是染色体数量文件。
文件内容是两列, 第一列是物种名(和进化树的相同); 第二列是该物种的染色体条数; 两列是tab分隔符:
A 18
B 20
C 30
D 24
2.安装/进入revbayes.
https://revbayes.github.io/download
我这边下载的singularity的镜像。
singularity run revbayes-v1.2.1-linux64-singularity.simg
进入到revbayes界面:
如图就是成了
输入example <- 1.0, 如果没报错就没问题了
运行之前确定一下工作路径
所有的输入文件都放在这个路径
如图
#然后在revBayes中输入以下所有代码
phylogeny <- readBranchLengthTrees("input.tree")[1]
#读入第一个输入文件即进化树
max_chromo = 30
#定义最大的祖先核型
chromo_data = readCharacterDataDelimited("data/aristolochia_chromosome_counts.tsv", stateLabels=(max_chromo + 1), type="NaturalNumbers", delimiter="\t", headers=FALSE)
#读入第二个输入文件即物种的核型数量
------------------
gamma ~ dnExponential(10.0)
delta ~ dnExponential(10.0)
rho ~ dnExponential(10.0)
#使用先验均值为0.1的指数先验来模拟染色体发生分支上沿着的多倍性和异倍性事件的发生率
#gamma是染色体获得率，delta是染色体丢失率，rho是多倍化率
mvi = 1
moves[mvi++] = mvScale(gamma, lambda=1, weight=1)
moves[mvi++] = mvScale(delta, lambda=1, weight=1)
moves[mvi++] = mvScale(rho, lambda=1, weight=1)
#Add MCMC moves for each of the rates
Q := fnChromosomes(max_chromo, gamma, delta, rho)
#为染色体进化模型创建速率矩阵
#在这里，使用一个简单的ChromEvol模型，该模型仅包括染色体获得、丢失和多倍化的速率
root_frequencies := simplex(rep(1, max_chromo + 1))
#Here we assume an equal prior probability for the frequency of chromosome numbers at the root of the tree. This does not mean that the frequencies are actually equal, we just give it an equal prior probability.
#Alternatively, we could have treated the root frequencies as a free variable and estimated them from the observed data. This approach will be illustrated in further examples
chromo_ctmc ~ dnPhyloCTMC(Q=Q, tree=phylogeny, rootFreq=root_frequencies,
type="NaturalNumbers")
chromo_ctmc.clamp(chromo_data)
#创建随机节点的染色体进化连续时间马尔可夫链（CTMC）
mymodel = model(phylogeny)
#模型的各项参数已经制定完毕
#把进化树弄进去
monitors[1] = mnScreen(printgen=10)
monitors[2] = mnJointConditionalAncestralState(filename="output/
ChromEvol_simple_anc_states.log", printgen=10, tree=phylogeny, ctmc=chromo_ctmc, type="NaturalNumbers")
monitors[3] = mnModel(filename="output/ChromEvol_simple_model.log", printgen=10)
#没看懂, 但跟着跑就行
mymcmc = mcmc(mymodel, monitors, moves)
mymcmc.run(200)
#在本例中，将链长度设置为200
#但是对于真实的分析，需要运行更多的迭代并检查收敛性
anc_state_trace = readAncestralStateTrace("output/ChromEvol_simple_anc_states.log")
#
ancestralStateTree(phylogeny, anc_state_trace, "output/ChromEvol_simple_final.tree", burnin=0.25, reconstruction="marginal")
#
#运行完毕
#在前面getwd()得到的路径底下多了一个文件夹: Output，里面就是输出文件
q() #退出
4.画图
三个输出文件
# Install the RevGadgets R package from:
# https://github.com/revbayes/RevGadgets
#
library(RevGadgets)
setwd("E:/workspace/ChromoSSE")
# plot the ancestral states tree generated by RevBayes
p = plot_ancestral_states(tree_file="ChromEvol_simple_final.tree",
summary_statistic="MAPChromosome",
include_start_states=FALSE,
tip_label_size=3,
tip_label_offset=0.002,
tip_label_italics=TRUE,
node_label_size=2.5,
node_label_nudge_x=0.0005,
xlim_visible=c(0, 0.105),
shoulder_label_size=2,
alpha=.3,
show_posterior_legend=TRUE,
title="A Basic ChromEvol Analysis")
# using ggtree and ggplot2 functions we can now add all sorts of details
# to our figure:
# annotate some of the clades
#p = p + geom_cladelabel(node=60, label="Isotrema", align=TRUE, offset=0.02, size=2)
#p = p + geom_cladelabel(node=36, label="Aristolochia", align=TRUE, offset=0.02, size=2)
# write pdf
print(p)
ggsave("ChromEvol_simple.pdf")
使用上感觉还是像PastML这种软件
就这样了
