老师，我想要全部显著相关的基因对

0.需求

这是我的直播课学员提的需求，觉得挺有意义的，就帮他实现了一下。

想要找出一个表达矩阵里所有相关性r>0.8且p<0.05的基因对。

不是直接从矩阵或者里看，而是得到若干对基因作为输出结果。

1.编一个表达矩阵

set.seed(10086)
exp = matrix(rnorm(600,sd = 10),nrow = 60)
rownames(exp) = paste0("gene",1:nrow(exp))
colnames(exp) = paste0("sample",1:ncol(exp))
exp[1:10,] = exp[1:10,] + 5
exp[1:4,1:4]

##         sample1    sample2   sample3    sample4
## gene1  10.49789   4.964393 -0.473901 13.6810594
## gene2 -22.44959  19.323965 15.036701 -0.6999084
## gene3  10.66458 -11.270693  3.908851 14.5247626
## gene4   9.85479   3.768001 -4.621285 10.5618159

boxplot(exp)

2.计算基因相关性和p值

cor.test函数不支持矩阵化计算，所以使用corrplot里的函数cor.mtest。

library(corrplot)

## corrplot 0.92 loaded

corm = cor(t(exp)) 
copm = cor.mtest(t(exp))$p
pheatmap::pheatmap(corm)

3.宽变长

两个矩阵分别宽变长。然后合并到一起

library(tidyverse)
corm2 = corm %>% 
  as.data.frame() %>% 
  rownames_to_column("G1") %>% 
  pivot_longer(cols = starts_with("gene"),
               names_to = "G2",
               values_to = "correlation")
copm2 = copm %>% 
  as.data.frame() %>% 
  rownames_to_column("G1") %>% 
  pivot_longer(cols = starts_with("gene"),
               names_to = "G2",
               values_to = "pvalue")
identical(copm2$G1,corm2$G1)
## [1] TRUE
identical(copm2$G2,corm2$G2)
## [1] TRUE
dat = mutate(corm2,pvalue = copm2$pvalue)
head(dat)
## # A tibble: 6 × 4
##   G1    G2    correlation pvalue
##   <chr> <chr>       <dbl>  <dbl>
## 1 gene1 gene1     1       0     
## 2 gene1 gene2    -0.513   0.129 
## 3 gene1 gene3     0.579   0.0795
## 4 gene1 gene4     0.427   0.218 
## 5 gene1 gene5     0.122   0.737 
## 6 gene1 gene6     0.00537 0.988

4.去重复

现在，G1-G2组成的基因对也是有重复的，给你看个例子就明白了。

dat[c(2,61),]
## # A tibble: 2 × 4
##   G1    G2    correlation pvalue
##   <chr> <chr>       <dbl>  <dbl>
## 1 gene1 gene2      -0.513  0.129
## 2 gene2 gene1      -0.513  0.129

所以这样的行只保留一行即可。

直接按照三四列去重也不是不行，但不怕一万就怕万一啊，基因数量多了的话那谁说得准的。。留下这个bug，万一以后踩雷了岂不学术造假。

观察宽变长的规律可以发现，基因的排列是有顺序的。第一列是由行名来的，重复了60次(111222333这样的)，第二列是由列名变来的，重复了六十轮(123123123这样的)。 所以思路就是把原来的矩阵对角线及以上的格子去掉即可。

dat$x = rep(1:60,each = 60)
dat$y = rep(1:60,times = 60)
# x<y或者x>y都行，反正只留一半。
dat2 = filter(dat,x<y)

5. 筛选符合要求的基因对

相关系数大于0.8且p值小于0.05

filter(dat2,abs(correlation)>0.8 & pvalue<0.05)
## # A tibble: 8 × 6
##   G1     G2     correlation    pvalue     x     y
##   <chr>  <chr>        <dbl>     <dbl> <int> <int>
## 1 gene3  gene55      -0.844 0.00213       3    55
## 2 gene14 gene20       0.804 0.00508      14    20
## 3 gene15 gene18      -0.859 0.00146      15    18
## 4 gene15 gene39      -0.832 0.00284      15    39
## 5 gene18 gene42      -0.803 0.00512      18    42
## 6 gene29 gene57       0.853 0.00169      29    57
## 7 gene31 gene60       0.812 0.00431      31    60
## 8 gene33 gene38       0.938 0.0000599    33    38