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

关于测试数据共享文件声明 百度云是一种非常方便的文件共享方式,但是有时会出现文件无法通过审核,导致大家访问失败?之前团队分享视频(百度管片最严,你懂的,上周六的纪录片将扩展名mkv修改...

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

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

image

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

attachments-2018-01-parNdWJ85a5df17108ae6.png

图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)

attachments-2018-01-rmnRxgkr5a5df186eb2a2.png

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

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

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

QIIME绘制稀释曲线

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

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

  1. 虚拟机安装:适合在Windows上学习,但分析效率低。
  2. Docker安装:Linux上最简单的安装方法,需要管理员帮忙并给你开通部分权限。
  3. 管理员直接安装:直接安装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
# 批量计算shannonobserved_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多样性指数,和不同的分组类别,有多种结果选择。 缺点是结果不是矢量图。

attachments-2018-01-WvvGC6Ug5a5df19a7e827.png

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

猜你喜欢

  • 发表于 2018-01-16 00:05
  • 阅读 ( 5396 )
  • 分类:其他组学

0 条评论

请先 登录 后评论
不写代码的码农
刘永鑫

工程师

64 篇文章

作家榜 »

  1. 祝让飞 118 文章
  2. 柚子 91 文章
  3. 刘永鑫 64 文章
  4. admin 57 文章
  5. 生信分析流 55 文章
  6. SXR 44 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章