可能是因为某个基因集相关的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)<br style="visibility: visible;">[1] 20 500 <br style="visibility: visible;">> dim(dat_lnc)<br style="visibility: visible;">[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_4m6a_1 1.46 -1.22 -1.26 -0.58m6a_2 -0.14 0.46 2.27 -0.62m6a_3 -1.96 -1.37 -0.32 -0.90m6a_4 -0.83 0.58 0.67 0.82> round(dat_lnc[1:4,1:4],2) sample_1 sample_2 sample_3 sample_4lnc_1 0.76 0.91 0.17 0.40lnc_2 0.47 -1.89 -0.10 0.62lnc_3 -0.57 -0.34 -1.07 -1.25lnc_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$rdata.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 variablescor_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_m6Am2=dat_lnccor_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循环嵌套得到的两个矩阵进行整合后筛选统计学显著的结果也非常简单。

转载自《生信技能树》，如有侵权，请联系删除。