【R语言编程】---利用三代测序绘制菌群聚类热图与物种丰度图

前言：仍然是三代测序数据的分析，宏基因组的文章中经常出现聚类热图和物种丰度图，用来直观地识别与某些疾病或者表型相关的菌群构成。

1.读取数据
一共有11个样本，每一个样本的测序reads都经过Nanopore官方的Epi2Me程序鉴定了物种，下表中第一列是被鉴定的菌种，第二列是该样本中每个物种产生的reads数目。

pecies	Sample1
Prevotella bivia	12159
Gemella asaccharolytica	9192
Mycoplasma hominis	1165
Veillonella montpellierensis	935
Parvimonas micra	931
Dialister micraerophilus	761
Anaerococcus tetradius	607
Aerococcus christensenii	552
Lactobacillus iners	263
Mycoplasma equirhinis	176
Gemella bergeri	152
Prevotella amnii	141
Peptoniphilus harei	114
Sneathia sanguinegens	114

首先导入到R语言中，合并所有样本到一个数据框：

#创建一个空列表
L <- list()
i <- 1
#样本数为11，共11个输入文件
spnum <- 11
#读取文件作为元素循环添加至列表中
for (i in c(1:spnum)) {
  sp <- dir('F://研究生学习/Work Job/Job_6/RawData2/Rd/')[i]
  setwd('F://研究生学习/Work Job/Job_6/RawData2/Rd/')
  spr <- read.csv(sp,header = TRUE,sep=',',stringsAsFactors = FALSE)
  L[[i]] <- spr
  i <- i+1
}
#根据物种名合并列表中的数据框
rawdf <- L[[1]]
for(i in c(1:(spnum-1))){
  rawdf <- merge(rawdf,L[i+1],all = TRUE,by = "Species")
}
rownames(rawdf)<-rawdf$Species
newdf <- rawdf[,-1]
#将NA转换成0
newdf[is.na(newdf)]<-0
#替换数字中的千分位符号“，”
for (x in c(1:spnum)) {
  newdf[,x]<-gsub(',',"",newdf[,x])
}
#将数据框转换为数值型，方便后续归一化处理
newdf <- as.data.frame(lapply(newdf,as.numeric))

2.绘制热图
经过上一步，我们得到了列名为样本，行名为菌种的reads数据框，然后就可以绘制热图，进行聚类分析了：

#测序Reads数据归一化
newdf_norm = newdf
rownum = dim(newdf)[1]
for(i in 1:rownum){
  sample_sum=apply(newdf, 1, sum)
  for(j in 1:spnum){
    newdf_norm[i,j]=newdf[i,j]/sample_sum[i]
  }
}
#修改行名为物种名
rownames(newdf_norm)<-rownames(rawdf)
#绘制聚类热图
pheatmap(newdf_norm,fontsize = 15,border_color="black")

绘制结果：

菌群聚类热图

3.绘制物种丰度图
丰度图，其实就是堆积图，把每个样本的reads数目转换为百分数，然后作图就可以了：

#绘制丰度图
#复制新的数据框
newAb <- rawdf[,-1]
#将NA值转换为0
newAb[is.na(newAb)]<-0
#替换数字中的千分位符号“，”
for (x in c(1:spnum)) {
  newAb[,x]<-gsub(',',"",newAb[,x])
}
#将数据框转换为数值型
newAb_num <-  as.data.frame(lapply(newAb,as.numeric))
#将每一列的reads数目转换为丰度百分数
for (y in c(1:spnum)) {
  newAb_num[,y] <- newAb_num[,y] / sum(newAb_num[,y])
}
rownames(newAb_num)<-rownames(newAb)
Taxonomy <- rownames(newAb_num)
txdf <- data.frame(newAb_num,Taxonomy)
txdf <- melt(txdf, id='Taxonomy')
names(txdf)[2] <- "sample_id"
#绘图
ggplot(txdf, aes(sample_id, fill=Taxonomy, value*100))+
  geom_col(position='stack')+
  labs(x='Samples', y='Relative Abundance (%)')+
  scale_y_continuous(expand=c(0, 0))+
  theme(axis.text.x=element_text(angle=45, hjust=1,size = 20))

绘制结果：

物种丰度图