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

Python017

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

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

1.读取数据

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

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

2.绘制热图

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

绘制结果:

3.绘制物种丰度图

丰度图,其实就是堆积图,把每个样本的reads数目转换为百分数,然后作图就可以了:

绘制结果:

生态位宽度(Bi)采用Colwell等(1971)加权修正的Levins指数

式中,Bi为物种i的生态位宽度值,其范围在[0-1]之间,值越大说明该物种生态位宽度越宽;Pij为物种i在第j资源状态下的个体数占该种所有个体数的比例,r为资源总数。

生态位重叠指数  生态位重叠指数采用Pianka指数(Pianka et al ., 1973):

式中, Oik 为物种i与物种k的生态位重叠指数,其值越大表示生态位重叠程度越高,取值范围为[0,1];Pij为物种i利用资源状态 j 的个体数占该种个体总数的比例,r为资源(季节或样点)的总数;Pkj为物种k利用资源状态j的个体数占该种个体总数的比例。 R 表示生态响应速率。

下面我们进行R语言操作:

library(spaa) ##加载包

df1=read.csv("物种丰度.csv") #读取丰度数据

这里的丰度数据为:行为样本,列为物种;

(1)计算生态位宽度:levins指数

b=as.data.frame(niche.width(df1, method ="levins")) #计算生态位宽度

colnames(b)="levins" #重新命名

由于spaa包中的niche.width()函数是没有/r的,如图一,公式中,因此继续用"levins"/r,(r为样点总个数)

b$levins1=b$levins/nrow(df1)   #计算加权的levins生态位指数,既公式一(图一)的生态位指数

参考文献:

Colwell RK, Futuyma DJ. 1971. On the measurement of Niche Breadth and Overlap. Ecology, 52(4): 567-576.

前言:接上一篇,很多文献中为了更直观的展示一个微环境中的菌群分布,常常将样本聚类与物种丰度同时展示。

1.数据结构

首先需要准备丰度数据表Abundance.txt和分组信息group.txt

丰度数据以样本为列名,以菌种为行名:

分组信息以列出了每个样本属于病例/对照组:

3.再画丰度图

丰度图其实就是堆叠图:

4.结果展示