相关性分析返回相关性系数的同时返回p值
这个分析需求已经不是第一次有人问我了，可能是因为某个基因集相关的lncRNA的数据分析策略深入人心吧。越来越多的人选择了它相关性分析。如果是2万多个蛋白质编码基因和2万多个lncRNA基因的相关性，计算量就有点可怕，不过几十个m6a基因或者小班焦亡基因去跟其它基因进行相关性计算，基本上还是绝大部分小伙伴可以hold住的。
首先创造模拟数据；
dat_m6A = matrix(rnorm(10000),nrow = 20)
dim(dat_m6A)
dat_lnc = matrix(rnorm(15000*ncol(dat_m6A)),nrow = 15000)
dim(dat_lnc)
可以看到是20个m6a基因，以及 1.5万个lncRNA的表达量矩阵，而且样品数量是500个；
> dim(dat_m6A)
[1] 20 500
> dim(dat_lnc)
[1] 15000 500
接下来，我们就开始对 dat_m6A 和 dat_lnc 两个矩阵的不同基因，进行相关性分析。
首先对它们进行简单的行列命名，如下所示：
colnames(dat_m6A) = paste0('sample_',1:ncol(dat_m6A))
rownames(dat_m6A) = paste0('m6a_',1:nrow(dat_m6A))
colnames(dat_lnc) = paste0('sample_',1:ncol(dat_lnc))
rownames(dat_lnc) = paste0('lnc_',1:nrow(dat_lnc))
> round(dat_m6A[1:4,1:4],2)
sample_1 sample_2 sample_3 sample_4
m6a_1 1.46 -1.22 -1.26 -0.58
m6a_2 -0.14 0.46 2.27 -0.62
m6a_3 -1.96 -1.37 -0.32 -0.90
m6a_4 -0.83 0.58 0.67 0.82
> round(dat_lnc[1:4,1:4],2)
sample_1 sample_2 sample_3 sample_4
lnc_1 0.76 0.91 0.17 0.40
lnc_2 0.47 -1.89 -0.10 0.62
lnc_3 -0.57 -0.34 -1.07 -1.25
lnc_4 -1.47 0.02 -1.33 -0.73
因为，这两个矩阵，都是完全随机的，所以后续进行相关性分析，理论上R值和p值都表现不好。
最简单的是 corr.test 函数：
而 corr.test 函数 来自于 psych 这个包：
## do corr.test
data.corr <- corr.test(dat_m6A, dat_lnc,
method="pearson", adjust="fdr")
data.r <- data.corr$r
data.p <- data.corr$p
虽然 corr.test 函数 方便，一下子得到了r值矩阵和配套的p值矩阵。很容易进行筛选。但是它运行速度实在是太慢了！
起码在我写完这个教程的时候，它还没有运行完毕。而且我注意到它的帮助文档里面写的是：For very large matrices (> 200 x 200), there is a noticeable speed improvement if confidence intervals are not found.
也仅仅是说，它认为 matrices (> 200 x 200), 是 very large ，简直是搞笑啊，对生物信息学数据分析完全不适用啊它这个评价标准。
R基础包stats里面的cor函数是效率最高的，但是缺p值
dat=rbind(dat_m6A,dat_lnc)
# 这个 cor 函数运行非常快
# Unfortunately, the function cor() returns only the correlation coefficients between variables
cor_dat=cor(t(dat))
cor_dat[1:4,1:4]
cor_df=cor_dat[1:nrow(dat_m6A),
(nrow(dat_m6A)+1):(nrow(dat_m6A)+nrow(dat_lnc))]
round(cor_df[1:4,1:4], 2)
library(reshape2)
cor_df=melt(cor_df)
head(cor_df)
colnames(cor_df)=c('m6A' ,'lncRNA' ,'cor')
R_thre=0.2 # 因为是模拟数据，所以迫不得已，设置了 0.2
cor_df=cor_df[abs(cor_df$cor) > R_thre,]
cor_df$R = ifelse(cor_df$cor > 0,'postive','negative')
table(cor_df$R)
table(as.character(cor_df$m6A))
因为是模拟数据，所以迫不得已，设置了R的阈值是 0.2 ，如下所示：
> cor_df
m6A lncRNA cor R
93795 m6a_15 lnc_4690 -0.2197739 negative
141146 m6a_6 lnc_7058 -0.2058085 negative
239808 m6a_8 lnc_11991 -0.2090980 negative
298828 m6a_8 lnc_14942 0.2139784 postive
可以看到，完全随机的数值，也是可以达到约0.2的相关性哦，不过，这里没有给出p对应的p值，并不能说是统计学显著的相关性哦。
另外，因为我这个教程里面并没有设置随机数种子，所以呢，你也基本上不可能得到跟我一模一样的结果哦。
两个apply循环嵌套
这个问题是粉丝提问，我让对方发给我了代码，我看了看， 虽然对方已经是很灵活应用了apply函数，以及unlist函数，而且还可以自己创造函数，比如下面的cor_2_matrix函数，但是代码冗余太多了。
可能是对 R基础包stats里面的cor函数 不熟悉，以为它只能是对两个向量进行相关性计算，其实它可以直接对一个表达量矩阵进行相关性计算。
m1=dat_m6A
m2=dat_lnc
cor_2_matrix <- function(m1,m2){
apply(m2 , 1, function(x){
# x = m2[2,]
unlist(apply(m1, 1,function(y){
# y = m1[1,]
cor(as.numeric(x),
as.numeric(y))
}))
})
}
rdf=cor_2_matrix( m1 , m2 )
这个速度其实超级快了。但是粉丝也是不知道如何设置p值，我给他进行了简单的修改，如下所示；
corP_2_matrix <- function(m1,m2){
apply(m2 , 1, function(x){
# x = m2[2,]
unlist(apply(m1, 1,function(y){
# y = m1[1,]
cor.test(as.numeric(x),
as.numeric(y))$p.value
}))
})
}
pdf=corP_2_matrix( m1 , m2 )
速度也是一般般，而且比较麻烦的就是两次得到了两个矩阵，需要整合！可以看到，同样的，因为是模拟数据，所以基本上相关性都很弱，而且p值不太可能是小于0.05的， 很难有统计学显著性。
> rdf[1:4,1:4]
lnc_1 lnc_2 lnc_3 lnc_4
m6a_1 0.03162697 0.0160654358 -0.041608606 -0.02066662
m6a_2 -0.02402029 -0.0002989616 0.025571886 0.02516795
m6a_3 -0.03827284 0.0072451576 -0.001705459 0.03600183
m6a_4 -0.01619382 -0.0743613971 0.028690169 0.01480454
> pdf[1:4,1:4]
lnc_1 lnc_2 lnc_3 lnc_4
m6a_1 0.4804326 0.72007536 0.3531637 0.6447899
m6a_2 0.5920671 0.99467954 0.5683611 0.5744888
m6a_3 0.3931174 0.87161827 0.9696560 0.4218182
m6a_4 0.7179333 0.09673156 0.5221350 0.7412269
这两个矩阵如何整合，筛选统计学显著的结果，就是后话了。
最辣鸡的两个for循环嵌套
当我把这个问题发在讨论群，让学员们尝试解决，发现绝大部分小伙伴给出来的都是最辣鸡的两个for循环嵌套，运行效率本身就堪忧，而且极度的不美观。
## 最垃圾的代码，两个for循环
if(F){
outTab=data.frame()
for(i in row.names(m6A_exp)){
for(j in row.names(lnc_exp)){
x=as.numeric(m6A_exp[i,])
y=as.numeric(lnc_exp[j,])
corT=cor.test(x,y)
cor=corT$estimate
pvalue=corT$p.value
if((cor>corFilter) & (pvalue<pvalueFilter)){
outTab=rbind(outTab,cbind(m6A=i,lncRNA=j,
cor,pvalue,Regulation="postive"))
}
if((cor< -corFilter) & (pvalue<pvalueFilter)){
outTab=rbind(outTab,cbind(m6A=i,lncRNA=j,
cor,pvalue,Regulation="negative"))
}
}
}
}
这个是最垃圾的代码，两个for循环，速度超级慢，只不过里面添加好了筛选标准而已。
我们前面的两个apply循环嵌套得到的两个矩阵进行整合后筛选统计学显著的结果也非常简单。
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你