Корреляция один ко многим в матрице в R

Попытка вычислить корреляцию конкретного гена (здесь его ген 1) со всеми другими (35999 кандидатов) в матрице

         samp1   samp2   samp3   samp4
Gene 1  3.7891  2.4487  1.1939  0.6013
Gene 2  1.4484  3.2316  2.841   1.9545
Gene 3  0.4505  2.6062  2.0729  0.6403
.
.
Gene 36000  1.8828  5.2633  2.7552  1.7335

Я использовал следующий код

library(Hmisc)
A <- read.table("C:/Users/Desktop/exp.txt", header=T, sep="\t")
cor <- rcorr(as.matrix(A), type="pearson")


write.csv(cor$r,'C:/Users/Desktop/pCC VALUES.csv')
write.csv(cor$P,'C:/Users/Desktop/p VALUES.csv')

Но приведенный выше код для расчета матрицы 36000*36000. Однако я хотел бы получить соотношение один ко многим, где интересующий ген всегда является первым геном. Вот его ген 1. Это сэкономит время обработки. Один из методов, очевидно, заключается в извлечении 36000 пар моего интереса из вывода. Я хотел бы знать, возможен ли какой-либо другой метод определения корреляции моего гена со всеми остальными без расчета многих на многие.

Редактировать:

Я ищу вывод, как в формате,

Gene 1 Gene 2   pcc  p-value
Gene 1 Gene 3   pcc  p-value
.
.
Gene 1 Gene 36000 pcc p-value
end

1 ответ

Решение

Если я правильно вас понял (корреляция первого ряда с каждым другим, взятым по одному за раз), то что-то в этом духе может помочь вам начать:

dat <- as.matrix(read.table(text = "samp1;samp2;samp3;samp4
Gene 1;3.7891;2.4487;1.1939;0.6013
Gene 2;1.4484;3.2316;2.841;1.9545
Gene 3;0.4505;2.6062;2.0729;0.6403
Gene 4;0.4705;2.4062;1.0729;0.6003
Gene 5;1.8828;5.2633;2.7552;1.7335", sep=";"))

corr_list <- list()

for (i in 2:nrow(dat)) {
  r <- cor.test(dat[1,], dat[i,])
  corr_list[[paste("Genes 1 &", i)]] <- c(r$estimate, p.val=r$p.value)
}


# Results
corr_list

$`Genes 1 & 2`
       cor      p.val 
-0.3070573  0.6929427 

$`Genes 1 & 3`
       cor      p.val 
-0.1417635  0.8582365 

$`Genes 1 & 4`
       cor      p.val 
0.04777015 0.95222985 

$`Genes 1 & 5`
      cor     p.val 
0.1425788 0.8574212 

Вы также можете поместить результаты в data.frame, если это более удобно:

corr_list <- data.frame(Gene1=numeric(), Gene2=numeric(), cor=numeric(), p.value=numeric())

for (i in 2:nrow(dat)) {
  r <- cor.test(dat[1,], dat[i,])
  corr_list[i-1,] <- c(1, i, r$estimate, r$p.value)
}

corr_list

  Gene1 Gene2         cor   p.value
1     1     2 -0.30705735 0.6929427
2     1     3 -0.14176355 0.8582365
3     1     4  0.04777015 0.9522299
4     1     5  0.14257884 0.8574212
Другие вопросы по тегам