вычисление корреляции всех против всех генов в R: Каков наилучший подход для этого?

Question

У меня 14000 генов (столбец: Gene) и 200 образцов (столбец: sample1 sample2 ...)

Я пытаюсь вычислить корреляции для ~ 14000 генов против всех и добавить все генные корреляции и необходимые столбцы из набора данных (test_df) в новом фрейме данных (df1) и запишите результаты в текстовый файл.

Когда я запускаю код, я получаю корреляции между (Gene1 и Gene2) и (Gene1 и Gene3) , Когда l oop приходит к 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. Этот способ поможет вам избежать двойного l oop. вы можете выбрать go с идентификатором вместо имени, если имя не уникально

Теперь нужно просто пройти через имя 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

Тогда мы можем завершить sh это с

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

Вам не нужно для l oop, функция 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

Как вы можете видеть, это соответствует верхнему треугольнику этой матрицы выше. Ниже я вытащу верхний треугольник матрицы и соединю его с этим вектором.

Итак, мы соединяем копейщика и Пирсона с другой информацией:

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

90 миллионов строк! Вы уверены, что храните такой огромный массив данных.

Answer 3

Во-первых, у меня есть неэффективный способ с циклами для вашей задачи:

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)