R语言之相关系数计算篇

R语言之相关系数计算篇

简介：

在环境微生物类的文章中，经常出现计算物种与基因、基因与基因、基因与代谢物之间的相关系数的内容，在这个计算的基础之上再进行相关的可视化。例如相关性热图、网络图等等。文献中常出现的相关系数有Spearman、Pearson两种。

案例：

之间课题组一个师兄想代谢组学中代谢物与基因之间的相关性，共选择了95种代谢物，3313个相关基因，三个实验组一个对照组（每组三个生物学重复，共计12个样本）如下：

想分析每种代谢物和基因之间的关系，需要求解P值和R值，如下：

求解思路：

方法一：cor() 函数与cor.test()函数

这两种函数是R自带的函数，使用简单，操作过程如下：

cor()函数

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1) #导入数据
R_pearson <- cor(metabolite,gene,method = 'pearson') #求解相关系数

结果：

成功求解相关系数，but，该函数不能求解P值，失败*1。

cor.test()函数

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
R_pearson <- cor.test(metabolite,gene,method = 'pearson')

结果：

> R_pearson <- cor.test(metabolite,gene,method = 'pearson')
Error in cor.test.default(metabolite, gene, method = "pearson") : 
  'x'和'y'的长度必需相同

cor.test()函数只能进行同维度的数据计算，但是相较于cor()函数，cor.test()函数可以计算出P值，如下：

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
R_pearson <- cor.test(metabolite$X2.Hydroxypyridine,gene$Cluster.1184.0,method = 'pearson') #只选择其中一列数据进行计算

结果：

But！我们需要计算的是不同维度的数据，失败*2

方法二：Hmisc包中的rcorr()函数以及Psych包中的corr.test()函数

Hmisc包中的rcorr()函数

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
library(Hmisc) #载入包（前提是你安装了这个包）
rcorr_test_pearson <- rcorr(as.matrix(gene),as.matrix(metabolite),type = 'pearson')#这里的编写规则和前两个不一样，数据需要转化为矩阵as.matrix()

结果：

该计算过程需要耗费一丢丢时间，当然不超过30s的样子

该结果包含了相关系数和P值，算是成功求解，但是一看结果：

R_pearson$r

该结果中包含了gene、metabolite本身的相关系数以及P值，原因是啥俺也不知道/(ㄒoㄒ)/~~，失败*3

Psych包中的corr.test()函数

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
library(psych) #载入包（前提是你安装了这个包）
R_pearson <- corr.test(gene,metabolite,method = 'pearson')

结果：

结果就是出不来结果，┭┮﹏┭┮。好像是因为数量较大和corr.test()函数还包括了P值矫正的过程，这里需要指明的是虽然我没有写出矫正的方法但是corr.test()函数会有一个默认的矫正过程，所以大家在使用corr.test()的时候一定要注意这里有默认矫正的过程！！！！

代码一直停不下来，最后我按了ESC中止了代码。

后续：

只运行部分数据进行求解：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
library(psych) #载入包（前提是你安装了这个包）
R_pearson <- corr.test(gene[1:3],metabolite[1:6],method = 'pearson')

后续的结果：

至少说明在数据量比较小的时候，corr.test()函数有着比较好的优势。失败*3

结果：

最后还是决定使用Hmisc包中的rcorr()函数进行求解：

代码：

metabolite <- read.csv('metabolite.csv',row.names = 1)
gene <- read.csv('gene.csv',row.names = 1)
library(Hmisc) #载入包（前提是你安装了这个包）
rcorr_test_pearson <- rcorr(as.matrix(gene),as.matrix(metabolite),type = 'pearson')#这里的编写规则和前两个不一样，数据需要转化为矩阵as.matrix()
rcorr_test_pearson_R <- rcorr_test_pearson$r[1:3313,3314:3408]
rcorr_test_pearson_P <- rcorr_test_pearson$P[1:3313,3314:3408]#后来琢磨了一下，我只选择结果里面我需要的数据即可
write.csv(rcorr_test_pearson_R,'R.csv')#导出数据
write.csv(rcorr_test_pearson_P,'P.csv')

结果：


芜湖！下班！么么哒！