计算所有与 R 中所有基因的相关性：最好的方法是什么？

Question

我有 14000 个基因（列：Gene）和 200 个样本（列：sample1 sample2 ...）

我正在尝试计算 ~14000 个基因的相关性，并将数据集 (test_df) 中的所有基因相关性和所需列附加到新数据框 (df1) 中，并将结果写入文本文件。

当我运行代码时，我得到了（Gene1 和 Gene2）和（Gene1 和 Gene3）之间的相关性。当循环到达 Gene2 时它会中断并且错误提示

cor.test.default(as.matrix(test_df[i, ][, 3:length(test_df)]) 中的错误： 没有足够的有限观察

我每行有 3 到 4 个值，这不应该是这种情况。

请建议任何有效的方法，因为我必须对 14000 个基因进行关联。如何在多个内核上运行此代码以更快地获得结果？

请在下面找到代码和生成的文件。

提前致谢

> test_df <- data.frame(ID=c("ID_3721", "ID_537", "ID_555"), 
                      Gene=c("Gene1","Gene2","Gene3"),
                      sample1=c(11397,78191,44838),
                      sample2=c(33768,33763,7680),
                      sample3=c(74521,33268,72367),
                      sample4=c(51486,11435,28772),
                      sample5=c(73539,21486,0))

> test_df
##       ID  Gene sample1 sample2 sample3 sample4 sample5
##1 ID_3721 Gene1   11397   33768   74521   51486   73539
##2  ID_537 Gene2   78191   33763   33268   11435   21486
##3  ID_555 Gene3   44838    7680   72367   28772       0

for(i in 1:2){
       for(j in i+1:3){

          p.cor <- cor.test(as.matrix(test_df[i,][,3:length(test_df)]), as.matrix(test_df[j,][,3:length(test_df)]), method="pearson")$estimate
          s.cor <- cor.test(as.matrix(test_df[i,][,3:length(test_df)]), as.matrix(test_df[j,][,3:length(test_df)]), method="spearman")$estimate

          df1 <- data.frame(ID1   = test_df[i,1],
                            ID2   = test_df[j,1],
                            Name1 = test_df[i,2],
                            Name2 = test_df[j,2],
                            correlation.p = p.cor
                            correlation.s = s.cor)

         write.table(df1, file="genecorr.txt", row.names=FALSE, sep="\t", append=TRUE, quote=FALSE, col.names = !file.exists("genecorr.txt"))

   }
}

**Error in cor.test.default(as.matrix(test_df[i, ][, 3:length(test_df)]),  : 
  not enough finite observations**

genecorr.txt

ID1     ID2     NAME1   NAME2    correlation.p      correlation.s
ID_3721 ID_537  Gene1   Gene2    -0.136733508500744  -0.1
ID_3721 ID_555  Gene1   Gene3    0.145998550191942    0.3

Answer 1

我建议首先通过以下方式转换您的数据

 dt <- dcast(melt(id.vars=c("ID","Gene"),test_df),variable~Gene)

setDT(dt)

## > dt
##    variable Gene1 Gene2 Gene3
## 1:  sample1 11397 78191 44838
## 2:  sample2 33768 33763  7680
## 3:  sample3 74521 33268 72367
## 4:  sample4 51486 11435 28772
## 5:  sample5 73539 21486     0




nameidx <- combn(names(dt)[-1],2)
 ## > nameidx
 ##      [,1]    [,2]    [,3]   
 ## [1,] "Gene1" "Gene1" "Gene2"
 ## [2,] "Gene2" "Gene3" "Gene3"

请注意使用函数combn 生成名称索引是多么容易。这种方式可以帮助您避免双循环。 如果名称不是唯一的，您可以选择使用 ID 而不是名称

现在只需通过名称 idx 即可

res  <- dt[,lapply(1:ncol(nameidx),
         function(x){ c(pearson=cor.test(get(nameidx[1,x]),
                                    get(nameidx[2,x]),method="pearson")$estimate,
         spearman=cor.test(get(nameidx[1,x]),
                           get(nameidx[2,x]),method="spearman")$estimate)})]

## >  > res
##            V1        V2        V3
## 1: -0.7411691 0.0394641 0.3444608
## 2: -0.6000000 0.1000000 0.3000000

然后我们就可以结束了

 ## > res1 <- setnames(data.table(cbind(t(nameidx),t(res))),c("Name1","Name2","pearson","spearman"))[]
 ## > res1
 ##    Name1 Name2            pearson spearman
 ## 1: Gene1 Gene2 -0.741169112323627     -0.6
 ## 2: Gene1 Gene3 0.0394640960151169      0.1
 ## 3: Gene2 Gene3  0.344460833012615      0.3

Answer 2

首先，对于您的问题，我有一种低效的循环方式：

test_df <- data.frame(ID=c("ID_3721", "ID_537", "ID_555"),
                      Gene=c("Gene1","Gene2","Gene3"),
                      sample1=c(11397,78191,44838),
                      sample2=c(33768,33763,7680),
                      sample3=c(74521,33268,72367),
                      sample4=c(51486,11435,28772),
                      sample5=c(73539,21486,0))

df1<-data.frame(ID1=0,ID2=0,Name1=0,Name2=0,correlation=0)

k<-1

for(i in 1:2){
       for(j in i:3){
       if(i!=j){
          p.cor <- cor.test(as.matrix(test_df[i,][,3:length(test_df)]), as.matrix(test_df[j,][,3:length(test_df)]), method="pearson")$estimate
          s.cor <- cor.test(as.matrix(test_df[i,][,3:length(test_df)]), as.matrix(test_df[j,][,3:length(test_df)]), method="spearman")$estimate


          df1[k,] <- c(as.character(test_df[i,1]),as.character(test_df[j,1]),as.character(test_df[i,2]),as.character(test_df[j,2]),as.character(p.cor))

                            k<-k+1
                            }
   }
}

也许这会快一点

n<-nrow(test_df)

fun<-function(y)cor(x,y)


result<-c()
for(i in 1:(n-1))
{
x<-as.numeric(test_df[i,3:ncol(test_df)])
result<-c(result,apply(test_df[(i+1):nrow(test_df),3:ncol(test_df)],1,fun))
}


m<-rep((n-1):1,(n-1):1)

a<-rep(test_df[,1][-n],(n-1):1)
b<-rep(test_df[,2][-n],(n-1):1)

c<-d<-numeric()
for(i in 2:n)
{
c<-c(c,as.character(test_df[,1][i:n]))
d<-c(d,as.character(test_df[,2][i:n]))
}

df1<-data.frame(ID1=a,ID2=c,Name1=b,Name2=d,correlation=result)

Answer 3

您不需要 for 循环，cor 函数可用于矩阵。默认情况下，它计算列之间的成对相关性，因此根据您的情况，转置矩阵：

rownames(test_df) = test_df[,2]
cor(t(test_df[,-c(1:2)]),method="pearson")
           Gene1      Gene2     Gene3
Gene1  1.0000000 -0.7411691 0.0394641
Gene2 -0.7411691  1.0000000 0.3444608
Gene3  0.0394641  0.3444608 1.0000000

其中一些是多余的，所以我们只取出上面的三角形。我们事先得到了比较的索引：

ind = which(upper.tri(cor(t(test_df[,-c(1:2)]))),arr.ind=TRUE)
     row col
[1,]   1   2
[2,]   1   3
[3,]   2   3

如您所见，这对应于上述矩阵的上三角形。下面我将把矩阵的上三角拉出来和这个向量连接起来。

所以我们加入了 spearman 和 pearson 以及其他信息：

cor_vector = function(M,Method){
res = cor(M,method=Method)
res[upper.tri(res)]
}

data.frame(
test_df[ind[,1],1:2],
test_df[ind[,2],1:2],
pearson = cor_vector(t(test_df[,-c(1:2)]),"pearson"),
spearman = cor_vector(t(test_df[,-c(1:2)]),"spearman")
)

             ID  Gene   ID.1 Gene.1    pearson spearman
Gene1   ID_3721 Gene1 ID_537  Gene2 -0.7411691     -0.6
Gene1.1 ID_3721 Gene1 ID_555  Gene3  0.0394641      0.1
Gene2    ID_537 Gene2 ID_555  Gene3  0.3444608      0.3

但是，我需要提醒您，对于您大小为 14000*200 的矩阵，此计算非常繁琐。如果我进行快速计算，您的输出数据框将是：

choose(14000,2)
[1] 97993000

9000 万行！你确定要存储这么大的data.frame吗？