根据目录中所有文件中存在的数据计算成对 spearman 等级相关性

Question

我正在尝试计算 Spearman 的排名相关性，其中每个实验的数据（带有名称和排名的 tsv）存储在目录中的单独文件中。

以下是输入文件的格式：

#header not present
#geneName   value
ENSMUSG00000026179.14   14.5648627685587
ENSMUSG00000026179.14   0.652158034413075
ENSMUSG00000026179.14   0.652158034413075
ENSMUSG00000026179.14   1.852158034413075
ENSMUSG00000026176.13   4.13033421794948
ENSMUSG00000026176.13   4.13033421794948
ENSMUSG00000026176.13   15.4344068144428
ENSMUSG00000026176.13   15.4344068144428
ENSMUSG00000026176.13   6.9563523670728
...

我的问题是键（基因名称）是重复的，每个实验文件都包含不同但重叠的基因名称集。我需要的是每对基因名称的交集，同时执行相关性和删除重复项，可能类似于以下伪代码：

# Find correlation for all possible pairs of input(i.e. files in directory)
files = list_Of_files("directory")
for(i in files) {
    for(k in files) {
    CommonGenes <- intersect (i,k)
    tempi <- removeRepetitive(i, CommonGenes) #Keep the gene with highest value and remove all other repeating genes. Also, keep only common genes.
    tempk <- removeRepetitive(k, CommonGenes) #Keep the gene with highest value and remove all other repeating genes. Also, keep only common genes. 
    correlationArray[] <- spearman(tempi, tempk) #Perform correlation for only the common genes
}
}

最后，我想使用 corrplot 或 qtlcharts 绘制相关矩阵。

Answer 1

这是另一种解决方案。它没有嵌套循环，而是使用expand.grid 创建组合，然后使用‹dplyr›动词管道来计算主表子集的相关性。

这种方法既有优点也有缺点。最重要的是，它非常适合“整洁数据”方法和there are some who advocate to work in tidy data as much as possible。实际代码和zx8754差不多。

library(dplyr)

genes = sprintf('ENSMUSG%011d', 1 : 50)
my_dfs = replicate(4, tibble(Gene = sample(genes, 20, replace = TRUE), Value = runif(20)),
                   simplify = FALSE)

首先，我们要使基因名称唯一，因为随后每张表都需要唯一的基因：

my_dfs = lapply(my_dfs, function (x) summarize(group_by(x, Gene), Value = max(Value)))

现在我们可以创建这个列表的所有排列：

combinations = bind_cols(expand.grid(i = seq_along(my_dfs), j = seq_along(my_dfs)),
                         expand.grid(x = my_dfs, y = my_dfs))

此时，我们有一个表，其中包含所有成对组合 i、j 的索引，以及组合本身作为列表列：

# A tibble: 16 x 4
       i     j                 x                 y
   <int> <int>            <list>            <list>
 1     1     1 <tibble [17 x 2]> <tibble [17 x 2]>
 2     2     1 <tibble [18 x 2]> <tibble [17 x 2]>
 3     3     1 <tibble [19 x 2]> <tibble [17 x 2]>
…

我们现在按索引分组，并按基因名称加入每个组中的单个列表列：

correlations = combinations %>%
    group_by(i, j) %>%
    do(inner_join(.$x[[1]], .$y[[1]], by = 'Gene')) %>%
    print() %>%
    summarize(Cor = cor(Value.x, Value.y, method = 'spearman'))

中场休息：在print() 行，我们留下了所有基因表的所有成对组合的完全扩展表（两个原始表的Value 列已重命名为Value.x 和Value.y ，分别）：

# A tibble: 182 x 5
# Groups:   i, j [16]
       i     j               Gene    Value.x    Value.y
   <int> <int>              <chr>      <dbl>      <dbl>
 1     1     1 ENSMUSG00000000014 0.93470523 0.93470523
 2     1     1 ENSMUSG00000000019 0.21214252 0.21214252
 3     1     1 ENSMUSG00000000028 0.65167377 0.65167377
 4     1     1 ENSMUSG00000000043 0.12555510 0.12555510
 5     1     1 ENSMUSG00000000010 0.26722067 0.26722067
 6     1     1 ENSMUSG00000000041 0.38611409 0.38611409
 7     1     1 ENSMUSG00000000042 0.01339033 0.01339033
…

下一行使用相同的组从这些表中简单地计算成对相关性。由于整个表格是长格式的，可​​以方便地用‹ggplot2›绘制：

library(ggplot2)

ggplot(correlations) +
    aes(i, j, color = Cor) +
    geom_tile() +
    scale_color_gradient2()

...但是如果您需要将其作为平方相关矩阵，那就再简单不过了：

corr_mat = with(correlations, matrix(Cor, nrow = max(i)))

      [,1]  [,2]  [,3]  [,4]
[1,]  1.00  1.00 -0.20 -0.26
[2,]  1.00  1.00 -0.43 -0.50
[3,] -0.20 -0.43  1.00 -0.90
[4,] -0.26 -0.50 -0.90  1.00

Answer 2

首先，将所有数据读入数据帧列表，更多信息请参见this post，这里我们只是创建一个虚拟数据。

library(dplyr)

# dummy data
set.seed(1)
myDfs <- list(
  data.frame(geneName = sample(LETTERS[1:4], 15, replace = TRUE), value = runif(15)),
  data.frame(geneName = sample(LETTERS[1:4], 15, replace = TRUE), value = runif(15)),
  data.frame(geneName = sample(LETTERS[1:4], 15, replace = TRUE), value = runif(15)),
  data.frame(geneName = sample(LETTERS[1:4], 15, replace = TRUE), value = runif(15)),
  data.frame(geneName = sample(LETTERS[1:4], 15, replace = TRUE), value = runif(15))
)

然后，就像你的两个嵌套 for 循环一样，我们这里有两个嵌套的 apply 函数。在循环中，我们正在聚合并获得匹配的合并基因名称的相关性。

res <- sapply(myDfs, function(i){
  # group by gene, get max value
  imax <- i %>% group_by(geneName) %>% summarise(i_Max = max(value))
  sapply(myDfs, function(j){
    # group by gene, get max value
    jmax <- j %>% group_by(geneName) %>% summarise(j_Max = max(value))
    # get overlapping genes
    ij <- merge(imax, jmax, by = "geneName")
    # return correlation
    cor(ij$i_Max, ij$j_Max, method = "spearman")
  })
})

res 会有相关矩阵。

res

#      [,1] [,2] [,3] [,4] [,5]
# [1,]  1.0 -0.2  1.0  0.4 -0.4
# [2,] -0.2  1.0 -0.2  0.8  0.0
# [3,]  1.0 -0.2  1.0  0.4 -0.4
# [4,]  0.4  0.8  0.4  1.0 -0.4
# [5,] -0.4  0.0 -0.4 -0.4  1.0

对于相关图有many alternatives to choose from。这里作为一个例子，我们使用 corrplot：

corrplot::corrplot(res)