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

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

首先导入到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))

绘制结果：

物种丰度图