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

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

首先，将所有数据读入数据帧列表，请参阅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 函数。在循环中，我们正在聚合并获得匹配的 merged 基因名称的相关性。

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)

Answer 2

这是一个替代解决方案。它没有使用嵌套循环，而是使用 expand.grid 来创建组合，然后使用 ‹dplyr› 动词管道来计算主 table.

子集上的相关性

这种方法既有优点也有缺点。最重要的是，它非常适合“整洁数据”方法，并且 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)

首先我们要使基因名称唯一，因为随后的一切都需要每个 table:

的唯一基因

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

此时，我们有一个table，其中包含所有成对组合的索引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() 行，我们留下了所有基因 table 的所有成对组合的完全扩展 table（Value 列原来的两个 table 分别重命名为 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
…

下一行使用相同的组从这些 table 中简单地计算成对相关性。由于整个 table 是长格式的，所以可以很方便地用 ‹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

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

Calculate pairwise spearman's rank correlation from data present in all files in a directory

statistics

r

bioinformatics

dataframe