Вычисление корреляционной матрицы с помощью пакета pspearman с функцией apply () - PullRequest
Новая
       84

Вычисление корреляционной матрицы с помощью пакета pspearman с функцией apply ()

0 голосов
/ 24 августа 2018

Я пытаюсь вычислить корреляцию Спирмена и значение p для фрейма данных. Для лучшего приближения p-значения я должен придерживаться пакета pspearman. Я ожидаю результата, аналогичного функции rcorr(). Но у меня проблема при выполнении pspearman:spearman.test() строка за строкой.

Мой фрейм данных содержит 5000 строк (генов) и 200 столбцов (пятен). И я хочу получить матрицу корреляции и матрицу р-значений для этих 5000 * 5000 пар генов-генов. Корреляция рассчитывается только тогда, когда оба гена не являются NA в более чем двух точках.

Я могу добиться этого с помощью циклов, но это слишком медленно для моего большого набора данных. У меня проблемы при попытке использовать apply(),sapply(),mapply() для повышения скорости.

Вот что я пробовал: 

data = data.frame(matrix(rbinom(10*100000, 50, .5), ncol=200))
dim(data) #5000, 200
rownames(data) <- paste("gene", 1:5000, sep="") 
colnames(data) <- paste("spot",1:200,sep='')

library(pspearman)
spearFunc = function(x,y=data) {
  df = rbind(x,y)
  # Check the number of complete spots.There are no NAs in this set.
  complete = sum(!(is.na(x)) & !(is.na(y)))
  if (complete >=2 ) {
    pspearman::spearman.test(as.numeric(x),as.numeric(y))
    # This function returns a list containing 8 values, like pvalue,correlation
    }}

pair.all1 = mapply(spearFunc,data,data)
dim(pair.all1)
# 8 200, 200 is the number of columns 
pair.all2 = apply(data,1,spearFunc)

Что приводит к ошибке:

Ошибка в pspearman :: spearman.test (as.numeric (x), as.numeric (y)): (список) объект не может быть приведен к типу 'double'

Я надеюсь использовать spearman.test для каждой пары генов с apply (), чтобы сделать 

spearman.test(data[gene1],data[gene1]) 
spearman.test(data[gene1],data[gene2])
....
spearman.test(data[gene1],data[gene5000])
...
spearman.test(data[gene5000],data[gene5000])

Он должен вернуть кадр данных из 8 строк и 25 000 000 столбцов (5000 * 5000 пар генов).

Можно ли использовать apply () внутри apply () для достижения моей цели?

Thx!

1 Ответ

0 голосов
/ 24 августа 2018

Рассмотрите возможность создания парных комбинаций генов из row.names с combn и затем итерации по списку пар через определенную функцию.Обязательно верните структуру NA из логики if, чтобы избежать вывода NULL в матричном выводе.

Однако, имейте в виду, что парные перестановки из 5000 генов (choose(5000, 2)) дают очень высокие значения при12 497 500 элементов!Следовательно, sapply (сам цикл) может не сильно отличаться по производительности от for.Посмотрите на распараллеливание итерации. 

gene_combns <- combn(row.names(data), 2, simplify = FALSE)

spear_func <- function(x) {
  # EXTRACT ROWS BY ROW NAMES  
  row1 <- as.numeric(data[x[1],])
  row2 <- as.numeric(data[x[2],]) 

  # Check the number of complete spots.There are no NAs in this set.
  complete = sum(!(is.na(x)) & !(is.na(y)))

  if (complete >=2 ) {
    pspearman::spearman.test(row1, row2)        
  } else {
    c(statistic=NA, parameter=NA, p.value=NA, estimate=NA, 
      null.value=NA, alternative=NA,   method=NA, data.name=NA)
  }
}

pair.all2 <- sapply(gene_combns, spear_func)

Тестирование

Выше был протестирован с cor.test (точно такие же входные аргументы и список вывода, что и * 1020)* spearman.test , но точнее p-value) с использованием небольшого образца набора данных (50 obs, 20 vars): 

set.seed(82418)
data <- data.frame(matrix(rbinom(10*100000, 50, .5), ncol=200))[1:50, 1:20]
rownames(data) <- paste0("gene", 1:50) 
colnames(data) <- paste0("spot", 1:20)

gene_combns <- combn(row.names(data), 2, simplify = FALSE)
# [[1]]
# [1] "gene1" "gene2"    
# [[2]]
# [1] "gene1" "gene3"    
# [[3]]
# [1] "gene1" "gene4"    
# [[4]]
# [1] "gene1" "gene5"    
# [[5]]
# [1] "gene1" "gene6"    
# [[6]]
# [1] "gene1" "gene7"

test <- sapply(gene_combns, spear_func)  # SAME FUNC BUT WITH cor.test
test[,1:5]

#             [,1]                              [,2]                             
# statistic   885.1386                          1659.598                         
# parameter   NULL                              NULL                             
# p.value     0.1494607                         0.2921304                        
# estimate    0.3344823                         -0.2478179                       
# null.value  0                                 0                                
# alternative "two.sided"                       "two.sided"                      
# method      "Spearman's rank correlation rho" "Spearman's rank correlation rho"
# data.name   "row1 and row2"                   "row1 and row2"                  
#             [,3]                              [,4]                             
# statistic   1554.533                          1212.988                         
# parameter   NULL                              NULL                             
# p.value     0.4767667                         0.7122505                        
# estimate    -0.1688217                        0.08797877                       
# null.value  0                                 0                                
# alternative "two.sided"                       "two.sided"                      
# method      "Spearman's rank correlation rho" "Spearman's rank correlation rho"
# data.name   "row1 and row2"                   "row1 and row2"                  
#             [,5]                             
# statistic   1421.707                         
# parameter   NULL                             
# p.value     0.7726922                        
# estimate    -0.06895299                      
# null.value  0                                
# alternative "two.sided"                      
# method      "Spearman's rank correlation rho"
# data.name   "row1 and row2"
...