Alpha多样性稀释曲线rarefractioncurve还不会画吗？快看此文
关于测试数据共享文件声明
百度云是一种非常方便的文件共享方式，但是有时会出现文件无法通过审核，导致大家访问失败？之前团队分享视频(百度管片最严，你懂的，上周六的纪录片将扩展名mkv修改为jpg才通过审核)时为避免被和谐，采用后台回复获取可更新下载链接的方式。而扩增子、宏基因组、网络教程的测序数据为方便大家使用，直接将链接贴在正文中，半年来累计下载转存超过5000次，为大家参考格式准备数据、自学分析提供了很大便利。最近百度更新了敏感词库，导致之前分享的测序数据全部被强制取消分享(连分享的AGCT纯文本数据都不小心组成了敏感词，防不胜防呀！)，给后续学习者带来不便见谅！
今后团队的教程中尽量展示输入文件的格式，即只读正文即可自己动手准备相关文件。对于必不可少的示例文件，一律采用后台回复方式获取下载链接(确认链接失效请文章下方留言通知我们更新)，本文相关数据暗号为“稀释曲线”。请读者仔细阅读原文，找到暗号，并准确回复获取网盘下载链接。请复制链接在电脑端正规浏览器中访问(如Chrome、IE、Firefox)。
经常有人用手机、微信访问百度链接无法打开，并误报链接失败，这是因为手机或微信的系统经常错误识别链接(如把链接前后的文字加入链接生成错误地址)，或是功能不完善无法访问一些百度云网址，去年下半年至少收到过20次以上误报链接失效，不是微信自动生成链接错误，就是复制了不准确的网址(地址前后包括文字或链接不完整)。再次声明，lobbying请在电脑端准确复制链接，并用Chrome等兼容性强的浏览器访问。如访问失败，请先核对地址链接是否正确，并尝试至少使用两种上面提到的正规浏览器访问获得的准确网址，确认无法访问后，采用文章留言或后台留言的方式联系团队成员更新链接，谢谢！
Alpha多样性之稀释曲线
比较菌群Alpha多样性，最常见的方法是采用箱线图组间比较，此外还存在另一种常见比较样品、组间多样性的方法就是稀释曲线(rarefraction curve)。
稀释曲线的使用实例，可以参考之前分享的一篇文章
《Microbiome: 简单套路发高分文章—杨树微生物组》，其中图1中4个子图全是稀释曲线。
今天我们将介绍最方便的方法是uarch + R，逐行实行绘制。
同时还提供QIIME命令行，和使用Docker QIIME实例网页版交互式稀释曲线的方法(前提是你能成功安装QIIME，或是有权限使用Docker)。
本实战，是基于本公众号之前发布文章 《扩增子分析教程-3统计绘图-冲击高分文章》中提供测试数据的OTU表、实验设计文件上进行分析，保存过的直接查找相应文件，新手请后才回复“稀释曲线”获得本文所需OTU表和实验分组信息文件。
采用Uarch获得alpha rarefraction的数值
这里我们需要一个文件的OTU表，格式如下：
# Constructed from biom file
apec是什么意思#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 --table-type="OTU table" --to-tsv
采用uarch进行重抽样等量标准，并基于各样品等量数据生成稀释梯度数据
# 重采样为相同reads数量，默认为1万，测试数据比较小，这是抽3000条
uarch10 -otutab_norm -output otu_ -sample_size 3000 # , normlize by subsample to 10000
# 取1%-100%的序列中OTUs数量
uarch10 -alpha_div_rare otu_ -output -method without_replacement
overwhelmingly最后我们获得了不同梯度下各样品中OTU的数据碰撞检测文件，内容格式如下：
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
关于be good atuarch10 -alpha_div_rare 的使用方法和介绍，请阅读其官方帮助及相关页面 /uarch/manual/cmd_alpha_div_rare.html
其实这个表你在Excel中选中，就可以生成稀释曲线的图。但为我们还需要进一步的统计，
以及为了实现可重复式计算，建议所有统计绘图全部在R语言中用代码记录，便于将来修改和检查错误，发表文章同时共享代码也是对成果的负责。
初三英语听力
R ggplot2绘制稀释曲线
1. 按每个样品抽样的richness绘制曲线
实验设计的基本格式
SampleID BarcodeSequence LinkerPrimerSequence ReverPrimer 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 - Choo Directory)
# 安装和载入相关包和数据文件
#source("/biocLite.R")
#biocLite(c("ggplot2","reshape2")) # 没安装过此类包的请手动运行这两行代码安装包
library("ggplot2") # load related packages
library("reshape2")
design = read.table("", header=T, row.names= 1, p="\t")
rare = read.table("", header=T, row.names= 1, p="\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 (Obrved OTUs)")+
scale_x_discrete(breaks = c(1:10)*10, labels = c(1:10)*10) + theme_classic()
胎教英语p
ggsave(paste("alpha_rare_samples.pdf", p=""), p, width = 8, height = 5)callyoutonight
图1. 拆线图展示27个样品，3000条序列，在1% - 100% 抽样条件下的稀释曲线，并按组上色。
