0

我正在努力建立一个基于成对基因数据的系统发育树。下面是我的数据子集（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列中获得基因

在此处输入图像描述

4

1 回答 1

2

这个问题的例子有一些解释的余地​​。我的回答只有在每个个体中确实只有两个基因并且每一行描述一个个体时才有效。但是，我认为，如果每一行都意味着肯定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) + ...
于 2014-02-17T22:51:40.490 回答