系统发育树

Question

我正在努力建立一个基于成对基因数据的系统发育树。下面是我的数据子集（test.txt）。树不必基于任何 DNA 序列构建，而只是把它当作文字。

ID  gene1   gene2

1   ADRA1D  ADK
2   ADRA1B  ADK
3   ADRA1A  ADK
4   ADRB1   ASIC1
5   ADRB1   ADK
6   ADRB2   ASIC1
7   ADRB2   ADK
8   AGTR1   ACHE
9   AGTR1   ADK
10  ALOX5   ADRB1
11  ALOX5   ADRB2
12  ALPPL2  ADRB1
13  ALPPL2  ADRB2
14  AMY2A   AGTR1
15  AR  ADORA1
16  AR  ADRA1D
17  AR  ADRA1B
18  AR  ADRA1A
19  AR  ADRA2A
20  AR  ADRA2B

下面是我在 R 中的代码

library(ape)
tab=read.csv("test.txt",sep="\t",header=TRUE)
d=dist(tab,method="euclidean")
fit <- hclust(d, method="ward")
plot(as.phylo(fit))

这里附上我的图

我有一个关于它们如何聚集的问题。因为对

 17 AR  ADRA1B
 18 AR  ADRA1A

和

 2  ADRA1B  ADK
 3  ADRA1A  ADK

应该紧密聚集，因为它们有一个共同的基因。所以17和2应该在一起，18和3应该在一起。

如果我使用这种方法有误（欧几里得距离），我应该使用其他方法吗？

我是否应该将我的数据转换为行和列矩阵，其中gene1是x轴，gene2是y轴，每个单元格填充1或0？（基本上，如果它们成对意味着1，并且如果不是那么 0)

更新代码：

   table=table(tab$gene1, tab$gene2)
   d <- dist(table,method="euclidean")
   fit <- hclust(d, method="ward")
   plot(as.phylo(fit))

但是，在这里我只得到来自gene1而不是gene2列的基因。下图正是我想要的，但也应该有来自gene2列的基因

Answer 1

这个问题的例子有一些解释的余地​​。我的回答只有在每个个体中确实只有两个基因并且每一行描述一个个体时才有效。但是，如果每一行都意味着gene1 与gene2 一起出现，那么我认为肯定不会执行有用的聚类。在这种情况下，我希望有一个额外的列来说明它们常见发生的概率，并且可能首选主成分分析 (PCA) 之类的东西，但我离成为（层次）聚类专家还很遥远。

在您可以使用dist 函数之前，您必须将数据转换为适当的格式：

# convert test data into suitable format
gene.names <- sort(unique(c(tab[,"gene1"],tab[,"gene2"])))
gene.matrix <- cbind(tab[,"ID"],matrix(0L,nrow=nrow(tab),ncol=length(gene.names)))
colnames(gene.matrix) <- c("ID",gene.names)
lapply(seq_len(nrow(tab)),function(x) gene.matrix[x,match(tab[x,c("gene1","gene2")],colnames(gene.matrix))]<<-1)

得到的gene.matrix的形状为：

     ID ACHE ADK ADORA1 ADRA1A ADRA1B ADRA1D ADRA2A ...
[1,]  1    0   1      0      0      0      1      0
[2,]  2    0   1      0      0      1      0      0
[3,]  3    0   1      0      1      0      0      0
[4,]  4    0   0      0      0      0      0      0
...

因此，每一行代表一个观察（=个体），其中第一列标识个体，如果基因存在，则随后的每一列包含1，如果基因缺失，则包含0。在这个矩阵上可以合理地应用dist 函数（ID 列被移除）：

d <- dist(gene.matrix[,-1],method="euclidean")
fit <- hclust(d, method="ward")
plot(as.phylo(fit))

也许，阅读距离度量 euclidean、manhattan 等之间的差异是个好主意。例如，ID=1 和 ID=2 的个体之间的欧几里得距离是：

euclidean_dist = sqrt((0-0)^2 + (1-1)^2 + (0-0)^2 + (0-0)^2 + (0-1)^2 + ...)

而曼哈顿距离是

manhattan_dist = abs(0-0) + abs(1-1) + abs(0-0) + abs(0-0) + abs(0-1) + ...