Calculer la corrélation de rang de Spearman à partir de données présentes dans tous les fichiers d'un répertoire

Question

J'essaie de calculer la corrélation de rang de Spearman, où les données (tsv avec le nom et le rang) pour chaque expérience sont stockées dans des fichiers distincts dans un répertoire.

Voici le format des fichiers d'entrée:

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

Mon problème est que les clés (noms de gènes) sont répétitives, et chaque fichier d'expérience contient jeu différent mais qui se chevauchent des noms de gènes. Ce que je besoin est une intersection de noms de gènes pour chaque paire tout en effectuant la corrélation et la suppression des doublons, probablement quelque chose comme ce code pseudo:

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

En fin de compte, je veux tracer la matrice de corrélation en utilisant corrplot ou qtlcharts.

Answer 1

Tout d'abord, lisez toutes les données dans une liste de données, voir this post pour plus d'informations, ici nous créons juste une donnée fictive.

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

Alors, comme vos deux imbriqués pour boucles, ce que nous avons ici est deux imbriqués appliquent fonctions. Au sein des boucles, nous sommes en train d'agréger et d'obtenir une corrélation sur les noms de gènes appariés.

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 aura la matrice de corrélation.

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 

Pour courbe de corrélation il y a many alternatives to choose from. Voici à titre d'exemple, nous utilisons corrplot:

corrplot::corrplot(res) 

Answer 2

est ici une solution de rechange. Plutôt que d'avoir une boucle imbriquée, il utilise expand.grid pour créer les combinaisons, puis utilise un pipeline de verbes pour calculer les corrélations sur un sous-ensemble de la table maître.

Cette approche présente à la fois des avantages et des inconvénients. Tout d'abord, il s'intègre bien dans l'approche "données soignées", et there are some who advocate to work in tidy data as much as possible. Le code actuel est à peu près aussi long que celui de 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) 

Tout d'abord nous voulons rendre les noms de gène unique, car tout nécessite ensuite des gènes uniques par table:

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

Maintenant, nous pouvons créer toutes les permutations de cette liste:

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

À ce stade, nous avons un tableau avec les indices de toutes les combinaisons par paires i, j, ainsi que les combinaisons elles-mêmes comme colonnes de la liste:

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

Nous maintenant groupe par les indices et rejoindre les colonnes de la liste unique dans chaque groupe par des noms de gènes:

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

Intermission : à la ligne print() nous nous retrouvons avec une table entièrement étendue de toutes les combinaisons par paires de toutes les tables de gènes (les colonnes Value des deux tables originales ont été renommées en Value.x et Value.y, respectivement):

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

La ligne suivante calcule de manière triviale des corrélations par paires à partir de ces tables, en utilisant les mêmes groupes. Étant donné que toute la table est au format long, il peut être facilement tracée avec :

library(ggplot2) 

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

... mais si vous avez besoin cela comme une matrice de corrélation carrée au lieu, rien de plus facile:

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