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