【原】技術(shù)貼 | R語言：組學(xué)關(guān)聯(lián)分析和pheatmap可視化

微生態(tài) 2021-04-13

展開全文

本文由阿童木根據(jù)實(shí)踐經(jīng)驗(yàn)而整理，希望對(duì)大家有幫助,。

原創(chuàng)微文,，歡迎轉(zhuǎn)發(fā)轉(zhuǎn)載,。

舉例展示R語言組學(xué)關(guān)聯(lián)分析的方法。宏基因組數(shù)據(jù)以KO-樣品豐度表為例。代謝組數(shù)據(jù)以metabolite-樣品豐度表為例,?；痉椒ㄊ怯肦語言psych包c(diǎn)orr.test函數(shù)進(jìn)行兩組數(shù)據(jù)的相關(guān)分析,，結(jié)果經(jīng)格式化后用pheatmap可視化得熱圖。

一,、模擬輸入

1. KO豐度表

代碼：

ko_abun = as.data.frame(matrix(abs(round(rnorm(200, 100, 10))), 10, 20))
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(200, 200, 10))), 10, 20))
colnames(metabo_abun) = paste("met", 1:20, sep="_")
rownames(metabo_abun) = paste("sample", 1:10, sep="_")
metabo_abun

圖 2

二,、相關(guān)分析函數(shù)

思路

參數(shù)一：other -> KO或其他組學(xué)豐度表
參數(shù)二：metabo -> 代謝物豐度表
參數(shù)三：route -> 輸出目錄【提前創(chuàng)建】
corr.test進(jìn)行兩組數(shù)據(jù)相關(guān)分析
用stringr split將ko-metabolite結(jié)果列拆成兩列
結(jié)果保留r_value p_value
顯著相關(guān)標(biāo)記“*”

library(psych)
library(stringr)

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

    # 整理結(jié)果
    pair=rownames(result)  # 行名
    result2=data.frame(pair, result[, c(2, 4)])  # 提取信息

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

    # 格式化結(jié)果【將細(xì)菌代謝物拆成兩列】
    result4=data.frame(str_split_fixed(result3$pair, "-", 2), result3[, c(2, 3)])
colnames(result4)=c("feature_1", "feature_2", "r_value", "p_value")

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

三、相關(guān)性分析

代碼

dir.create("Result")  # 創(chuàng)建結(jié)果目錄
correlate(ko_abun, metabo_abun, "Result")

結(jié)果Correlation_result.txt如下,，行數(shù)為20X20

圖3

四,、pheatmap可視化函數(shù)

思路：

參數(shù)一：infile ->Correlation_result.txt相關(guān)性分析結(jié)果
參數(shù)二：route -> Route輸出路徑
用reshape2 dcast函數(shù)把feature_1 feature_2 p_value做成matrix，作為pheatmap輸入文件
用reshape2 dcast函數(shù)把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的相對(duì)豐度
    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="/"))
}