Alpha多样性稀释曲线rarefraction curve还不会画吗？快看此文

原创 metagenome 宏基因组

关于测试数据共享文件声明

百度云是一种非常方便的文件共享方式，但是有时会出现文件无法通过审核，导致大家访问失败？之前团队分享视频(百度管片最严，你懂的，上周六的纪录片将扩展名mkv修改为jpg才通过审核)时为避免被和谐，采用后台回复获取可更新下载链接的方式。而扩增子、宏基因组、网络教程的测序数据为方便大家使用，直接将链接贴在正文中，半年来累计下载转存超过5000次，为大家参考格式准备数据、自学分析提供了很大便利。最近百度更新了敏感词库，导致之前分享的测序数据全部被强制取消分享(连分享的AGCT纯文本数据都不小心组成了敏感词，防不胜防呀！)，给后续学习者带来不便见谅！

今后团队的教程中尽量展示输入文件的格式，即只读正文即可自己动手准备相关文件。对于必不可少的示例文件，一律采用后台回复方式获取下载链接(确认链接失效请文章下方留言通知我们更新)，本文相关数据暗号为“稀释曲线”。请读者仔细阅读原文，找到暗号，并准确回复获取网盘下载链接。请复制链接在电脑端正规浏览器中访问(如Chrome、IE、Firefox)。

经常有人用手机、微信访问百度链接无法打开，并误报链接失败，这是因为手机或微信的系统经常错误识别链接(如把链接前后的文字加入链接生成错误地址)，或是功能不完善无法访问一些百度云网址，去年下半年至少收到过20次以上误报链接失效，不是微信自动生成链接错误，就是复制了不准确的网址(地址前后包括文字或链接不完整)。再次声明，请在电脑端准确复制链接，并用Chrome等兼容性强的浏览器访问。如访问失败，请先核对地址链接是否正确，并尝试至少使用两种上面提到的正规浏览器访问获得的准确网址，确认无法访问后，采用文章留言或后台留言的方式联系团队成员更新链接，谢谢！

Alpha多样性之稀释曲线

比较菌群Alpha多样性，最常见的方法是采用箱线图组间比较，此外还存在另一种常见比较样品、组间多样性的方法就是稀释曲线(rarefraction curve)。

稀释曲线的使用实例，可以参考之前分享的一篇文章
《Microbiome: 简单套路发高分文章—杨树微生物组》，其中图1中4个子图全是稀释曲线。

今天我们将介绍最方便的方法是usearch + R，逐行实行绘制。

同时还提供QIIME命令行，和使用Docker QIIME实例网页版交互式稀释曲线的方法(前提是你能成功安装QIIME，或是有权限使用Docker)。

本实战，是基于本公众号之前发布文章《扩增子分析教程-3统计绘图-冲击高分文章》中提供测试数据的OTU表、实验设计文件上进行分析，保存过的直接查找相应文件，新手请后才回复“稀释曲线”获得本文所需OTU表和实验分组信息文件。

采用Usearch获得alpha rarefraction的数值

这里我们需要一个文件的OTU表，格式如下：

# Constructed from biom file
#OTU ID KO6     WT8     KO2     KO5     OE2  
OTU_52  73.0    15.0    29.0    62.0    1.0  
OTU_12  1109.0  182.0   575.0   886.0   110.0
OTU_26  87.0    53.0    66.0    117.0   12.0

也可以使用biom命令转换通常的biom为txt文件，此步不需要的请跳过：

# 获得文本OTU表
biom convert -i otu_table.biom -o otu_table.txt --table-type="OTU table" --to-tsv

采用usearch进行重抽样等量标准，并基于各样品等量数据生成稀释梯度数据

# 重采样为相同reads数量，默认为1万，测试数据比较小，这是抽3000条
usearch10 -otutab_norm otu_table.txt -output otu_table_norm.txt -sample_size 3000 # , normlize by subsample to 10000
# 取1%-100%的序列中OTUs数量
usearch10 -alpha_div_rare otu_table_norm.txt -output alpha_rare.txt -method without_replacement

最后我们获得了不同梯度下各样品中OTU的数据alpha_rare.txt文件，内容格式如下：

richness        KO6     WT8     KO2     KO5     OE2     WT1     WT2   
1       39.0    27.0    39.0    37.0    33.0    27.0    34.0    31.0  
2       70.0    44.0    71.0    64.0    55.0    57.0    69.0    62.0  
3       92.0    60.0    91.0    91.0    70.0    76.0    85.0    90.0  
4       107.0   70.0    114.0   108.0   90.0    99.0    104.0   102.0 
5       128.0   84.0    119.0   127.0   101.0   106.0   125.0   116.0 
6       134.0   90.0    135.0   145.0   115.0   119.0   141.0   131.0

关于usearch10 -alpha_div_rare 的使用方法和介绍，请阅读其官方帮助及相关页面 http://www.drive5.com/usearch/manual/cmd_alpha_div_rare.html

其实这个表你在Excel中选中，就可以生成稀释曲线的图。但为我们还需要进一步的统计，以及为了实现可重复式计算，建议所有统计绘图全部在R语言中用代码记录，便于将来修改和检查错误，发表文章同时共享代码也是对成果的负责。

R ggplot2绘制稀释曲线

1. 按每个样品抽样的richness绘制曲线

实验设计的基本格式

SampleID        BarcodeSequence LinkerPrimerSequence    ReversePrimer   genotype        compartment     batch   species Description
KO1     TAGCTT  AACMGGATTAGATACCCKG     ACGTCATCCCCACCTTCC      KO      root    1       Aarabidopsis    KnockDown
KO2     GGCTAC  AACMGGATTAGATACCCKG     ACGTCATCCCCACCTTCC      KO      root    1       Aarabidopsis    KnockDown
KO3     CGCGCG  AACMGGATTAGATACCCKG     ACGTCATCCCCACCTTCC      KO      root    1       Aarabidopsis    KnockDown

在Rstudio中，首先设置工作目录为测试文件所在目录(Session - Set Working Directory - Choose Directory)

# 安装和载入相关包和数据文件
#source("https://bioconductor.org/biocLite.R")
#biocLite(c("ggplot2","reshape2")) # 没安装过此类包的请手动运行这两行代码安装包
library("ggplot2") # load related packages
library("reshape2")
design = read.table("design.txt", header=T, row.names= 1, sep="\t") 
rare = read.table("alpha_rare.txt", header=T, row.names= 1, sep="\t") 

# 提取样品组信息
sampFile = as.data.frame(design$genotype,row.names = row.names(design))
colnames(sampFile)[1] = "group"

# 转换宽表格为ggplot通用长表格格式
rare$x = rownames(rare) # 添加x轴列
rare_melt = melt(rare, id.vars=c("x")) # 转换为长表格
rare_melt$x = factor(rare_melt$x, levels=1:100) # 设置x轴顺序

# 添加分组信息
rare_melt3 = merge(sampFile,rare_melt, by.x="row.names", by.y="variable")
rare_melt3$variable=rare_melt3$Row.names

# 按样品分组，按组上色
p = ggplot(rare_melt3, aes(x = x, y = value, group = variable, color = group )) + 
  geom_line()+
  xlab("Rarefraction Percentage")+ylab("Richness (Observed OTUs)")+
  scale_x_discrete(breaks = c(1:10)*10, labels = c(1:10)*10) + theme_classic()
p
ggsave(paste("alpha_rare_samples.pdf", sep=""), p, width = 8, height = 5)

图1. 拆线图展示27个样品，3000条序列，在1% - 100% 抽样条件下的稀释曲线，并按组上色。

这里有一个问题，是曲线为锯齿。主要原因是测试数据量太少，数据量越大越平滑。另一个解决办法是增加抽样次数，如将多次计算的抽样在表取均值。或者也可以增加显示的步长，比如将1%的步长改为2%或5%，也可以消除据锯齿，读者可以自己试试。

2. 按组求均值和标准误展示组间差异

# 求各组均值
# 读取usearch rarefraction文件，上面己经修改，必须重新读入
rare = read.table("alpha_rare.txt", header=T, row.names= 1, sep="\t") 
mat_t = merge(sampFile, t(rare), by="row.names")[,-1]
# 按第一列合并求均值
mat_mean = aggregate(mat_t[,-1], by=mat_t[1], FUN=mean)
# 修正行名
mat_mean_final = do.call(rbind, mat_mean)[-1,]
geno = mat_mean$group
colnames(mat_mean_final) = geno

rare=as.data.frame(round(mat_mean_final))
rare$x = rownames(rare)
rare_melt = melt(rare, id.vars=c("x"))

# 求各组标准误
# 转置rare表格与实验设计合并，并去除第一列样品名
se = function(x) sd(x)/sqrt(length(x)) # function for Standard Error
mat_se = aggregate(mat_t[,-1], by=mat_t[1], FUN=se) # se 为什么全是NA
mat_se_final = do.call(rbind, mat_se)[-1,]
colnames(mat_se_final) = geno

rare_se=as.data.frame(round(mat_se_final))
rare_se$x = rownames(rare_se)
rare_se_melt = melt(rare_se, id.vars=c("x"))

# 添加标准误到均值中se列
rare_melt$se=rare_se_melt$value

rare_melt$x = factor(rare_melt$x, levels=c(1:100))

p = ggplot(rare_melt, aes(x = x, y = value, group = variable, color = variable )) + 
  geom_line()+
  geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.5) +
  xlab("Percentage")+ylab("Richness (Observed OTUs)")+
  theme(axis.text.x=element_text(angle=90,vjust=1, hjust=1))+
  scale_x_discrete(breaks = c(1:10)*10, labels = c(1:10)*10)+
  theme_classic()
p
ggsave(paste("alpha_rare_groups.pdf", sep=""), p, width = 8, height = 5)

图2. 拆线图展示3组均值和标准误在1% - 100% 抽样条件下的稀释曲线，并按组上色。

对于组间无差异也是很正常的，不是所有的实验都在物种多样性差异。

大部分结果按组+标准误是不是差别更明显了，谁多谁少一目了解。对于其中的OE为直线是测试数据量太少造成的，实际公司返回的数据都有3-10万。一般取3-5万重抽样，过低的重复可以直接扔掉。

QIIME绘制稀释曲线

前提是你能装好QIIME，最近QIIME更新了Conda方式安装，是很方便，可是我在Ubuntu16.04上快速安装成功，但就是不好使。

我还是采用了之前的手动安装，或Docker安装的方式都成功了，有需要的可以参考之前发布过QIIME1.9.1安装的三种方法如下：

虚拟机安装：适合在Windows上学习，但分析效率低。
Docker安装：Linux上最简单的安装方法，需要管理员帮忙并给你开通部分权限。
管理员直接安装：直接安装QIIME1.9.1相关的上百个程序和包，不同环境依赖关系不同，需要极丰富经验，建议管理员安装。

QIIME进行OTU表步长重抽样，计算Alpha指数

# 调置重抽样的起、始和步长，默认每步还抽10次增加曲线平滑
multiple_rarefactions.py -i otu_table.biom -m 100 -x 3000 -s 100 -n 30 -o a_rare/ # min, max, step
# 批量计算shannon和observed_otus指数
alpha_diversity.py -i a_rare/ -m shannon,observed_otus -o a_div/ # -t ${result}/rep_seqs.tree 
# 合并结果
collate_alpha.py -i a_div/ -o a_collated/
# 实验设计前添加#符合QIIME mappingfile要求
cat <(echo -n "#") design.txt > design_rare.txt # make mapping file
# 绘制释释曲线
make_rarefaction_plots.py -i a_collated/ -m design_rare.txt -o a_rare_plots/

如果你QIIME安装不完整，出现报错Segmentation fault (core dumped)，可以在VirtualBox中运行，代码不变。推荐使用Docker更方便。

Docker虚拟机中运行生成稀释曲线的方法

docker run --rm -v `pwd`:/home --name=qiime yoshikiv/basespace-qiime-191-dev make_rarefaction_plots.py -i home/a_collated/ -m home/design_rare.txt -o home/

结果可打开rarefaction_plots.html网页方部，图片在average_plots目录中。

网页非常好用，可以选择不同的Alpha多样性指数，和不同的分组类别，有多种结果选择。
缺点是结果不是矢量图。

图3. QIIME的梯度稀释曲线，步长只设置了30个，曲线更平滑，而usearch为100次。

写在后面

宏基因组公众号经过半年发展，原创文章过百，关注人数12000+。为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外100+ PI，900+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

继续滑动看下一个