引言
上一期内容我们介绍了基于R语言对微生物丰度进行计算并进行可视化,今天这期内容主要介绍如何使用热图(Heatmap)来展示物种丰度的计算结果!
正文——代码
工作目录设置及R包的加载
rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\物种丰度计算及可视化')#设置工作路径
#安装包
install.packages("pheatmap")
#加载R包
library (pheatmap)
读取数据
df1 <- read.table(file="Genus.txt",sep="\t",header=T,check.names=FALSE)
head(df1)
数据处理
##利用循环处理具有重复的数据
data<-aggregate(E ~ Tax,data=df1,sum)
colnames(data)[2]<-"example"
for (i in colnames(df1)[2:length(colnames(df1))]){
#计算每列的和
data1<-aggregate(df1[,i]~Tax,data=df1,sum)
colnames(data1)[2]<-i
#合并列
data<-merge(data,data1,by="Tax")
}
df2<-data[,-2]
rownames(df2)=df2$Tax#修改行名
df3=df2[,-1]#删除多的列
#计算物种总丰度并降序排列
df3$rowsum <- apply(df3,1,sum)
df4 <- df3[order (df3$rowsum,decreasing=TRUE),]
df5 = df4[,-6]#删除求和列
#求物种相对丰度
df6 <- apply(df5,2,function(x) x/sum(x))
#取前10行
df7 <- df6[1:10,]
df8 <- 1-apply(df7, 2, sum) #计算剩下物种的总丰度
#合并数据
df9 <- rbind(df7,df8)
row.names(df9)[11]="Others"
#导出数据
write.table (df9, file ="genus_x.csv",sep =",", quote =FALSE)
绘图
# 行列注释信息
annotation_col = data.frame( g1 = c("a","a","b","c","c"))#行注释矩阵
rownames(annotation_col) = colnames(df9)
annotation_row = data.frame( g2 = factor(rep(c("c", "d", "e"), c(3, 4, 4))))#列注释矩阵
rownames(annotation_row) = rownames(df9)
#分组颜色
colors = list(g1 = c(a = "#1B9E77", b = "#D95F02",c = "red"),
g2 = c( c= "blue", d = "#E7298A", e = "#66A61E"))
#绘图
pheatmap(df9,
display_numbers = matrix(ifelse(df9 > 0.5, "*", ""), nrow(df9)), #设置条件并进行展示
annotation_col = annotation_col, annotation_row = annotation_row, #加入行列注释信息
annotation_colors = colors, #行列注释颜色设置
angle_col = "45", #行标签倾斜角度
cellwidth=25, cellheight=15, #格子大小
cluster_rows=F, treeheight_col = 15, #树的高度
legend_breaks=c(0,0.2,0.4,0.6,0.8) #图注标签展示
)