2014-02-17 65 views
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)) 

我的数字是附在这里enter image description here

我有一个问题,他们是如何clustered.Since成对

17 AR ADRA1B 
18 AR ADRA1A 

2 ADRA1B ADK 
3 ADRA1A ADK 

应密切聚集,因为他们有一个共同的gene.so 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)) 

然而,在此我只得到从基因1的基因,而不是基因2如下图所示column.The正是我想要的,但应该从基因2的基因列以及

enter image description here

+2

我有点不解。如何计算字符串的欧式距离?您的示例中的gene1和gene2列是否真的是字符串或因素?如果他们是因子,并且'dist'计算了因子水平上的欧氏距离,我认为没有什么合理的预期。 – Georg

+0

我不确定,但它看起来像是基于特定分类群中存在的基因的* union *进行聚类,而我认为'hclust'会基于每个*的身份进行聚类*基因 - 即,如果分类1具有'基因1 = A','基因2 = B'且分类2具有'基因2 = B','基因2 = A',则它们根本不匹配... –

+0

@Georg那些是字符串,我正在寻找方法从这个表中获得一些关系,以便有一些聚类并构建一棵树。我同意欧几里德不能被使用,我只是想举一个例子,关于我想要的。 – Rgeek

回答

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

也许,它是读取了距离测量euclideanmanhattan等之间的差异是个好主意例如,对于ID=1ID=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) + ... 
+0

我用这个命令:table(tab $ gene1,tab $ gene2)来获得矩阵。我想绘制基因的系统发生关系,而不是成对。我想看看他们是如何分别聚类而不是成对的。 – Rgeek

+0

组成员的数量由聚类算法决定。我认为您的(示例)数据中缺少一些信息以正确执行分层聚类。每行应包含个体中存在的基因的所有字符串。 – Georg

+0

我能够获得我正在寻找的系统发育图。正如我所说我想聚集单个基因(而不是配对),所以不得不稍微修改一下代码。但是,也可以这样做有趣。任何方式来获得表格信息哪些是聚集在一起的(簇号)?所以我可以有一个表中有簇号和它的分配项目,在我们计算出“fit < - hclust(d,method =”ward“ )” – Rgeek