Alpha多样性稀释曲线rarefraction curve还不会画吗？快看此文
关于测试数据共享文件声明
百度云是一种非常方便的文件共享方式，但是有时会出现文件无法通过审核，导致大家访问失败？之前团队分享视频(百度管片最严，你懂的，上周六的纪录片将扩展名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次。
猜你喜欢
扩增子分析：1图表解读 2分析流程 3统计绘图 4美化进化树ggtree
16S预测功能 1KEEG通路PICRUSt 2元素循环FAPROTAX 3表型bugbase 4KO通路Tax4Fun
宏基因组分析：1拼接MEGAHIT 2基因注释Prokka 3Binning 4MetaPhlAn2
系列教程：1生信入门Biostar 2微生物组 3宏基因组
写在后面
宏基因组公众号经过半年发展，原创文章过百，关注人数12000+。为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外100+ PI，900+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”