【绘图进阶】之通路与菌的相关性分析热图(六)

原创生信团上海天昊生物

前段时间，我们介绍了几种heatmap图的绘制，很多老师对其中的颜色、标记和分组有兴趣，在后台咨询的比较多。有很多老师咨询微生物多样性数据中菌和菌、菌与功能以及菌与环境因子之间相关性的分析，也有做转录组的老师咨询差异基因或代谢产物的热图绘制。这里，我们以多样性的数据为例，展示数据提取的方法、apply使用方法以及heatmap图的细节调整。

一

数据读取

1. 读取通路数据和物种丰度数据

In [1]:
df_kegg = read.table('kegg_for_stamp_level2.xls',header = T,sep = ' ',row.names = 1,check.names = F)
head(df_kegg)

Out[1]:

In [2]:
df_taxon = read.table('phylum.taxon.Abundance.xls',header = T,sep = ' ',row.names = 1,check.names = F)head(df_taxon)

Out[2]:

2. 物种数据过滤
过滤掉'All'行和'Abundance','superkingdom'列
最大值是 All行Abundance列

In [3]:
df_taxon = subset(df_taxon,Abundance>1000 & Abundance < max(df_taxon$Abundance) ,select = c(-Abundance,-superkingdom))head(df_taxon)

Out[3]:

二

相关性矩阵计算

spearman是非参数的，不要求数据必须是正态分布，也不要求方差具有齐性。

In [4]:
df_cor = apply(df_taxon,1,function(x){ apply(df_kegg,1,function(y)cor(x,y, method = "spearman")) })
head(df_cor)

Out[4]:

In [5]:
df_pvalue = apply(df_taxon,1,function(x){ apply(df_kegg,1,function(y)cor.test(x,y, method = "spearman")$p.value) })head(df_pvalue)

Out[5]:

三

绘制相关性热图

:
In [6]:
library(pheatmap)

1. 普通heatmap图
cluster 调节行和列是否聚类
fontsize可以调节字体大小；

In [7]:
pheatmap(df_cor,cluster_rows = T, cluster_cols = T,fontsize_row = 7, fontsize_col = 7)

Out[7]:

2. 添加显著性标记
In [8]:
display = df_pvalue
display[which(display < 0.001)] = "***"
display[which(display < 0.01)] = "**"
display[which(display > 0.01 & display < 0.05)] = "*"
display[which(display >= 0.05)] = ""
head(display)

Out[8]:

In [9]:
pheatmap(df_cor, fontsize_row = 8, fontsize_col = 5, cluster_rows = T, cluster_cols = T, display_numbers = display)

Out[9]:

3.仅绘制P值显著的结果

In [10]:
df_pvalue = subset(df_pvalue,apply(df_pvalue,1,min)<0.05)
df_pvalue

Out[10]:

In [11]:
df_cor = df_cor[rownames(df_pvalue),]
df_cor

Out[11]:

In [12]:
display = df_pvalue
display[which(display < 0.001)] = "***"
display[which(display < 0.01)] = "**"
display[which(display > 0.01 & display < 0.05)] = "*"
display[which(display >= 0.05)] = ""

修改字体大小，修改cell宽度和高度，添加border颜色，添加title，直接生成pdf文件，设置pdf中图片的宽度和高度。

In [13]:
pheatmap(df_cor, fontsize_row = 8, fontsize_col = 5,cellwidth = 15,cellheight = 40,border_color = 'RED',
cluster_rows = T, cluster_cols = T, display_numbers = display,
main='phylum.taxon_kegg.Correlation_p0.05', file='phylum.taxon_kegg.Correlation_p0.05.pdf',width = 10,height=10)

Out[12]: