如何从成对距离矩阵生成 Newick 树输出

Question

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

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

# 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 我用过的包： 猿, 帕角, phytools, 盖革

此代码似乎有效 -

# 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)

输出树：

Answer 1

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

不过，话虽如此，该任务仍然可以使用phangorn package 完成。例如，您可以根据距离矩阵创建分裂光谱（即矩阵中可能存在的分裂 (more details here - paywalled)。

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）简单地估计树。

Answer 2

正如@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 运行的长度需要根据链收敛进行调整。首先，代码应该很好地了解数据可以推断的关系。