1

我想从遗传数据中产生系统发育树。我在 R 和 python 中发现了一些看起来很棒的树图包，例如 R 中的 ggtree。但是这些需要已经采用树格式的数据输入，例如 Newick。

我认为大多数人从 vcf 文件开始并生成 FASTA 文件，但我的出发点是基因型表 - 我使用单倍体生物体，因此每个位置都是 0（参考）或 1（非参考）。据此，我使用 R 中的 dist() 计算成对遗传距离。 5 个样本 AE 的示例数据，成对距离超过十个变异位置：

# Generate dataframe with example genotypes
Variant <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
A <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0)
B <- c(1, 1, 0, 0, 1, 1, 0, 0, 1, 1)
C <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 1)
D <- c(1, 0, 1, 1, 0, 0, 1, 1, 0, 1)
E <- c(1, 0, 0, 0, 1, 1, 0, 0, 1, 1)
df = data.frame(Variant, A, B, C, D, E)
df
#  Remove first column with variant names
df$Variant <- NULL
# Transpose the columns and rows
df_t = t(df)
# Compute pairwise distance (Euclidean)
pdist = dist(df_t, method = "euclidean", diag = TRUE, upper = TRUE, p = 2)
pdist

我想从 pdist 生成一个分层树输出文件，例如 Newick 格式，以便我可以将其插入 ggtree 之类的包中以绘制漂亮的树，例如圆形系统发育图、用协变量着色等。我尝试过搜索但不确定从哪儿开始。

编辑/更新 这个网站很有帮助 http://www.phytools.org/Cordoba2017/ex/2/Intro-to-phylogenies.html 我使用了包：ape、phangorn、phytools、geiger

这段代码似乎工作 -

# Produce dendrogram
hclust = hclust(pdist)
# Check dendrogram looks sensible
plot(hclust)
class(hclust) # check that class is hclust
# Save to Newick file
my_tree <- as.phylo(hclust) 
write.tree(phy=my_tree, file="ExampleTree.newick") # Writes a Newick file
# Produce tree
plot(unroot(my_tree),type="unrooted",cex=1.5,
     use.edge.length=TRUE,lab4ut="axial",
     edge.width=2,
     no.margin=TRUE)

输出树： 来自示例数据的无根距离树

4

2 回答 2

2

这是一项不平凡的任务。要从距离矩阵构建一棵树（如在分叉树中），您将需要使用系统发育算法，并且最好不要从距离矩阵中构建（请注意，使用欧几里得距离作为二元矩阵也可能存在缺点）。

但是，也就是说，仍然可以使用phangornpackage完成任务。例如，您可以从距离矩阵创建分裂光谱（即矩阵中可能存在的分裂（更多细节在这里 - 付费专区）。

require(phangorn)
# Generate dataframe with example genotypes
Variant <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
A <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0)
B <- c(1, 1, 0, 0, 1, 1, 0, 0, 1, 1)
C <- c(0, 0, 1, 1, 0, 0, 1, 1, 0, 1)
D <- c(1, 0, 1, 1, 0, 0, 1, 1, 0, 1)
E <- c(1, 0, 0, 0, 1, 1, 0, 0, 1, 1)
df = data.frame(Variant, A, B, C, D, E)
df
#  Remove first column with variant names
df$Variant <- NULL
# Transpose the columns and rows
df_t = t(df)
# Compute pairwise distance (Euclidean)
pdist = dist(df_t, method = "euclidean", diag = TRUE, upper = TRUE, p = 2)

# calculate the Hadamard distance spectrum
distances <- distanceHadamard(as.matrix(pdist))
# representing the distances
lento(distances)
# Plotting the distances as a tree (a network actually)
plot(as.networx(distances), "2D")

请注意，在同一个包neighborNet中也可以使用，但手册强调此功能是实验性的。我建议联系包作者以获取更多信息。

然后，您可以将您的网络转换为"phylo"可供ape并且可能ggtree通过强制使用的方式：

# Converting into a phylo object
phylo <- as.phylo(distances)

但同样，请注意，这个生成的树在系统发育意义上可能是不正确的（即假设经过修改下降），我强烈建议使用基于模型的方法（例如使用MrBayes或BEAST2）简单地估计树。

于 2018-10-08T07:06:48.610 回答
0

正如@thomas-guillerme 所提到的，二进制数据可以有效地用于使用 MrBayes 构建系统发育树。输入文件应包含二进制data块和mrbayes命令。

#nexus
begin data;
dimensions ntax = 5 nchar = 10;
format datatype = restriction;
matrix
A 0011001100
B 1100110011
C 0011001101
D 1011001101
E 1000110011;
end;

begin mrbayes;
lset coding = variable;
mcmc ngen = 1000000 samplefreq = 1000;
sump burnin = 200;
sumt burnin = 200;
end;

运行的长度mcmc需要根据链收敛进行调整。首先，代码应该很好地了解数据可以推断的关系。

于 2018-10-11T15:09:58.230 回答