R语言计算两组数据相关性及对应p值的方法——psych包和Hmisc包

1. psych包的corr.test()函数 2. Hmisc包的rcorr()函数

在生信分析中,经常会用到R来计算两组的相关性,常用的有psych包的corr.test()函数和Hmisc包的rcorr()函数


一、psych包的corr.test()函数

psych这个包里的corr.test()函数是可以直接计算两个数据集变量之间的相关性的,这个结果里也有显著性检验的p值

输入文件格式见例子中的df1与df2

library(psych)
cor.result <- corr.test(df1,df2,method = method,adjust = adjust)
result.r <- cor.result$r
result.p <- cor.result$p

但是如果有时候数据量非常大,比如有上千上万个,用上述方法计算的时候是非常慢的

这时候我们可以使用另外一个函数,是Hmisc包中的rcorr()函数


二、Hmisc包的rcorr()函数

Hmisc包中的rcorr()函数速度快很多,但是不能计算两个数据集之间变量的相关性。这样的话可以先计算,然后再筛选。

这个函数要求的输入数据是矩阵格式,注意数据最少有5行

rcorr()这个函数的帮助文档也写到了

a numeric matrix with at least 5 rows and at least 2 columns (if y is absent). For print, x is an object produced by rcorr.

另外的一个知识点:如果想要用某个包里的函数,有两种办法,

第一种办法是先使用library()函数加载这个包,然后直接输入函数名;

另外一种办法是不加载,直接使用包名+两个冒号+函数

library(Hmisc)
rcorr(as.matrix(df1),as.matrix(df2))
或者
Hmisc::rcorr(as.matrix(df1),as.matrix(df2))


三、例子

比如在代谢组与转录组联合分析中,会分析差异基因与差异代谢物的相关性,通常的做法是利用表达量数据计算皮尔逊相关系数,然后设置一定的阈值进行筛选,再进行一些图形展示。


以下用一个例子带大家了解一下R计算相关性常用的两种方法:

我先模拟两个数据集,有两组样本A和B,每组重复3次,即A1,A2,A3;B1,B2,B3。

 df1是差异基因在两组样本中的表达量,df2是差异代谢物在两组样本中的表达量。

# df1是差异基因在不同样本中的表达量
df1<-data.frame(Var1=rnorm(10),
                Var2=rnorm(10),
                Var3=rnorm(10),
                Var4=rnorm(10),
                Var5=rnorm(10),
                Var6=rnorm(10))
colnames(df1) <- c("A1","A2","A3","B1","B2","B3")
rownames(df1) <- c("gene1","gene2","gene3","gene4","gene5","gene6","gene7",
                    "gene8","gene9","gene10")
#df2是差异代谢物在不同样本中的表达量
df2<-data.frame(Var7=rnorm(15),
                Var8=rnorm(15),
                Var9=rnorm(15),
                Var10=rnorm(15),
                Var11=rnorm(15),
                Var12=rnorm(15))
colnames(df2) <- c("A1","A2","A3","B1","B2","B3")
rownames(df2) <- c("meta1","meta2","meta3","meta4","meta5","meta6","meta7",
                    "meta8","meta9","meta10","meta11","meta12","meta13",
                    "meta14","meta15")

以下分别是构建的差异基因与差异代谢物的数据 df1和df2

attachments-2023-03-IbkQSlLA6411a1e4cd0ba.pngattachments-2023-03-Q1SrwFt36411a1f2214b8.png

处理数据时R包识别的是两组数据的行要一致,所以进行转置
#转置
gene <- as.data.frame(t(df1)) 
meta <- as.data.frame(t(df2)) 
以下分别是转置后的数据 gene和meta
attachments-2023-03-WduXBOyr6411a2720294e.png
attachments-2023-03-0yRWc4VS6411a2816ab6a.png


接下来我们用这个数据带大家分析一下

1.psych包的corr.test()函数

library(psych)
cor.result <- corr.test(meta,gene,method = "pearson", adjust = "adjust") #计算相关性

计算出来之后是个数据集,我们分别查看一下

r <- cor.result$r     #相关性
r

attachments-2023-03-bKrqNoEs6411a096990ea.png

p <- cor.result$p     #相关性的pvalue值
p

attachments-2023-03-jGlpmzLj6411a0b000f4e.png

结果是很方便观看的宽型数据,我们可以直接输出。

write.csv(r,"correlation.csv")    #输出相关性结果
write.csv(p,"pvalue.csv") #输出pvalue结果


2.Hmisc包的rcorr()函数

Hmisc包中的rcorr()函数速度快很多,但是不能计算两个数据集之间变量的相关性。这样的话可以先计算,然后再筛选。

这个函数要求的输入数据是矩阵格式,注意数据最少有5行

result <- Hmisc::rcorr(as.matrix(df1),as.matrix(df2))

attachments-2023-03-1QpKYUct6411a31fa59fe.png

自定义函数将这个结果转换成一个四列的数据框格式

flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
} #自定义函数
res1 <- flattenCorrMatrix(res.cor$r,res.cor$P) #输出四列结果

这生成的四列结果部分如下,可以看到Hmisc包中的rcorr()函数对每一个数据都互相进行了相关性和p值的计算:

attachments-2023-03-HPSr3p006411a394f36bc.png

所以,我们要用变量名去匹配,只保留我们想要的两组变量之间的相关性和p值结果。

geneid<-data.frame(gene_id=colnames(gene))
metaid<-data.frame(gene_id=colnames(meta))

attachments-2023-03-udn8NLXu6411a415888c0.pngattachments-2023-03-cfR5Ca0V6411a44f60ebe.png

library(dplyr)                                                   #加载dplyr包
res.final<- res1 %>%
merge(geneid,by.x="row",by.y = "gene_id") %>%
merge(metaid,by.x = "column",by.y = "meta_id") #用变量名去匹配,只保留两组变量之间的相关性和p值结果
colnames(res.final) <- c("meta","gene","cor","cor.pvalue") #更改列名

attachments-2023-03-KdaEt2pn6411a5536a6d2.png

至此,我们使用Hmisc包中的rcorr()函数来计算相关性和p值也已经完成了,更适用于大量数据的相关性分析。

这是一个典型的长型数据,如果想要更方便观察,我们也可以转换为宽型数据


关于使用tidyr包进行长宽数据转换的教程链接:R语言用tidyr包进行数据的长宽转换 - 组学大讲堂问答社区 (omicsclass.com)

library(tidyr)
res.cor <- res.final[,1:3] #提取前三列,相关性结果
wide.cor <- spread(res.cor, key = "gene", value = "cor") #转换为宽型数据
res.p <- res.final[,c("meta","gene","cor.pvalue")] #提取meta,gene,pvalue值对应的列,即pvalue结果
wide.p <- spread(res.p, key = "gene", value = "cor.pvalue") #转换为宽型数据

下图分别是两组变量之间的相关性和p值结果:

attachments-2023-03-kjyFtrZF64127f78c13ae.png

attachments-2023-03-aJzxZ5Ar64127ff960710.png



参考:https://www.jianshu.com/p/3cbb31a4b71f


  • 发表于 2023-03-07 13:14
  • 阅读 ( 7638 )
  • 分类:R

0 条评论

请先 登录 后评论
星莓
星莓

生物信息工程师

58 篇文章

作家榜 »

  1. omicsgene 654 文章
  2. 安生水 325 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. CORNERSTONE 72 文章
  7. rzx 67 文章
  8. xun 66 文章