技术贴 | R语言:组学关联分析和pheatmap可视化

2021
03/22

+
分享
评论
微生态
A-
A+

基本方法是用R语言psych包corr.test函数进行两组数据的相关分析,结果经格式化后用pheatmap可视化得热图。


本文由微科盟阿童木根据实践经验而整理,希望对大家有帮助。

原创微文,欢迎转发转载。


举例展示R语言组学关联分析的方法。宏基因组数据以KO-样品丰度表为例。代谢组数据以metabolite-样品丰度表为例。基本方法是用R语言psych包corr.test函数进行两组数据的相关分析,结果经格式化后用pheatmap可视化得热图。


一、模拟输入

 

1. KO丰度表

代码:

ko_abun = as.data.frame(matrix(abs(round(rnorm(20010010))), 1020))
colnames(ko_abun) = paste("KO"1:20, sep="_")
rownames(ko_abun) = paste("sample"1:10, sep="_")
ko_abun


 图 1


2. metabolite丰度表

代码:

metabo_abun = as.data.frame(matrix(abs(round(rnorm(20020010))), 1020))
colnames(metabo_abun) = paste("met"1:20, sep="_")
rownames(metabo_abun) = paste("sample"1:10, sep="_")
metabo_abun



  图 2


二、相关分析函数


思路

  1. 参数一:other -> KO或其他组学丰度表

  2. 参数二:metabo -> 代谢物丰度表

  3. 参数三:route -> 输出目录【提前创建】

  4. corr.test进行两组数据相关分析

  5. 用stringr split将ko-metabolite结果列拆成两列

  6. 结果保留r_value p_value

  7. 显著相关标记“*”


library(psych)
library(stringr)

correlate = function(other, metabo, route)
{
    #  读取方式:check.name=F, row.names=1, header=T
    # 计算相关性:
    result=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="fdr", alpha=.05, ci=TRUE, minlength=50), short=FALSE, digits=5))

    # 整理结果
    pair=rownames(result)  # 行名
    result2=data.frame(pair, result[, c(24)])  # 提取信息

    # P值排序
    result3=data.frame(result2[order(result2[,"raw.p"], decreasing=F),])

    # 格式化结果【将细菌代谢物拆成两列】
    result4=data.frame(str_split_fixed(result3$pair, "-"2), result3[, c(23)])
colnames(result4)=c("feature_1""feature_2""r_value""p_value")

    # 保存提取的结果
    write.table(result4, file=paste(route, "Correlation_result.txt", sep="/"), sep="\t", row.names=F, quote=F)
}

三、相关性分析

代码

dir.create("Result")  # 创建结果目录
correlate(ko_abun, metabo_abun, "Result")


结果Correlation_result.txt如下,行数为20X20



图3


四、pheatmap可视化函数

思路:

  1. 参数一:infile ->Correlation_result.txt相关性分析结果

  2. 参数二:route -> Route输出路径

  3. 用reshape2 dcast函数把feature_1 feature_2 p_value做成matrix,作为pheatmap输入文件

  4. 用reshape2 dcast函数把feature_1 feature_2 r_value做成matrix,作为pheatmap display

 

代码:

library(reshape2)
library(pheatmap)

correlate_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="\t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]

    # 用"*"代替<=0.05的p值,用""代替>0.05的相对丰度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            data_mark[i,j]=ifelse(data_p[i,j] <= 0.05"*""")
        }
    }

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.png", sep="/"))
}


五、pheatmap绘制热图

代码:

correlate_pheatmap("Correlation_result.txt", "Result")


结果目录,新增结果图(一个png一个pdf):



图4

 

打开pdf,如下:

图5

 

随机模拟的数据,没有显著的不奇怪。如果是做项目一般两组数据的相关分析都可以得到一些显著相关的结果,如下:



图6



本文由作者自行上传,并且作者对本文图文涉及知识产权负全部责任。如有侵权请及时联系(邮箱:nanxingjun@hmkx.cn
关键词:
R语言,关联,技术,分析,函数,参数

人点赞

收藏

人收藏

打赏

打赏

我有话说

0条评论

0/500

评论字数超出限制

表情
评论

为你推荐

推荐课程


社群

精彩视频

您的申请提交成功

确定 取消
剩余5
×

打赏金额

认可我就打赏我~

1元 5元 10元 20元 50元 其它

打赏

打赏作者

认可我就打赏我~

×

扫描二维码

立即打赏给Ta吧!

温馨提示:仅支持微信支付!