扩增子分析解读3格式转换,去冗余,聚类

大家好，欢迎来到IT知识分享网。

本网对Markdown排版支持较差，请跳转“宏基因组”公众号阅读；

写在前面

之前发布的《扩增子图表解读》系列，相信关注过我的朋友大部分都看过了(链接直达7月文章目录)。这些内容的最初是写本实验室的学生们学习的材料，加速大家对同行文章的解读能力。

《扩增子分析解读》系列文章介绍

扩增子分析是目前宏基因组研究中最常用的技术，由于微生物组受环境影响大，实验间重复较差，更需要更多的实验重复和分析技术来保证结果的准确性、可重复性。

本系统文章叫分析解读，即有详细的扩增子分析流程代码，又有本人对使用参数、备选参数意义的解读，可以让大部分零基础的人，更好的理解数据分析过程，并可亲自实践在自己的课题上，获得更好、更合理的实验结果。

本文采用目前最主流的扩增子测序数据类型HiSeq2500 PE250类型数据为例，结合目前主流方法QIIME+USearch优点组合定制的分析流程。本课程中所需的测序数据、实验设计和课程分析生成的中间文件，均可以直去百度云下载。链接：http://pan.baidu.com/s/1hs1PXcw 密码：y33d。

本课程代码的运行，至少需要Linux平台+安装QIIME1.9.1，我之前发布过三种安装QIIME的方法详见文章目录，总有一款适合你。

第三节. 去冗余,聚类,生成OTU表

本节课程，需要完成扩增子分析解读1质控,实验设计,双端序列合并和2提取barcode,质控及样品拆分,切除扩增引物

分析前准备

# 进入工作目录 cd example_PE250

上一节回顾：我们提取barcode,质控及样品拆分,切除扩增引物，经历了两节课6步数据处理才拿到我们扩增的高质量目的片段(貌似基因组/RNA-Seq测序结果直接就是这个阶段了，可以直接mapping)。

接下来我们将这些序列去冗余、聚类为OTU、再去除嵌合体，这样就可以获得高质量的OTU(类似于参考基因组/转录组)，用于定量分析每个OTU的丰度。这一阶段我们使用著名的扩增子分析流程Usearch。

Usearch简介

软件作者不仅有Usearch一款软件，它的Muscle(多序列比对，引用18659+4212次), Uparse(OT聚类算法，引用1529次), Uchime(扩增子嵌合体检测，引用3558次)等众多流行工具，个人引用超4万次，而且发的软件大多由作者一人完成，佩服。

Usearch安装

# 下载的程序并重命名：下载链接来自邮件，请用户自行复制邮件中地址替换下面代码中的网址；或者在windows里面下载并重命名为usearch10 wget -O "usearch10" http://drive5.com/cgi-bin/upload3.py?license=XXXXXX # 添加可执行权限 chmod +x usearch10 # 运行程序测试，成功可显示程序版本、系统信息和用户授权信息 ./usearch10

7. 格式转换

做生信为什么要学Python/Perl/Shell这些语言，主要原因是各软件间要求的具体格式不同，需要进行格式转换，才能继续运行。因此想成为高手，不会语言基本寸步难行。

我们现在将QIIME拆分的结果类型，要转换成Usearch要求的格式。常见的解决思路是读Usearch帮助看它的格式要求，写个Python/Perl脚本转换格式。我这里使用了Shell脚本一行解决，优点是快，但缺点很多(人不容易看懂、不同Linux系统shell版本不同可能失效)。

我们要转换的序列文件其实一直是fasta格式，只是序列名称行格式不同

# 目前格式  >KO1_0 HISEQ:419:H55JGBCXY:1:1101:1931:2086 1:N:0:CACGAT orig_bc=TAGCTT new_bc=TAGCTT bc_diffs=0 # Usearch要求的格式  >KO1_0;barcodelabel=KO1; 

# 格式转换 sed 's/ .*/;/g;s/>.*/&&/g;s/;>/;barcodelabel=/g;s/_[0-9]*;$/;/g' temp/PE250_P5.fa > temp/seqs_usearch.fa

上面这条命令有点复杂。sed是linux的一条命令，又是一种语言，擅长文本替换。替换的思路分四步：首先s/ ./;/g将原文件空格后面的内容(全是无用信息)替换为分号；其次s/>./&&/g是将序列名重复一次；再次s/;>/;barcodelabel=/g将重复后的;>替换为;barcodelabel=；最后s/_[0-9]*;$/;/g替换序列编号为分号。这只是我的思路，分析数据如解答数学题，可以有多种解法，你够聪明还会想出更好的解法。
新人一定感觉这命令每句都不像人话，我告诉你Perl和Shell就是这样–难读但高效。改用易读的Python语言，肯定没有Shell简洁。

8. 去冗余

# 序列去冗余 ./usearch10 -fastx_uniques temp/seqs_usearch.fa -fastaout temp/seqs_unique.fa -minuniquesize 2 -sizeout

计算过程中出现如下信息：

00:06 607Mb 100.0% Reading temp/seqs_usearch.fa 00:06 574Mb CPU has 96 cores, defaulting to 10 threads 00:08 915Mb 100.0% DF 00:09 935Mb  seqs,  uniques,  singletons (90.9%) 00:09 935Mb Min size 1, median 1, max 18774, avg 1.85 62167 uniques written,  clusters size < 2 discarded (26.6%)

本条命令的详细使用，请阅读官方文档 http://www.drive5.com/usearch/manual/cmd_fastx_uniques.html

9. 聚类OTU

# 聚类OTU ./usearch10 -cluster_otus temp/seqs_unique.fa -otus temp/otus.fa -uparseout temp/uparse.txt -relabel Otu

程序运行过程会显示运行时间、进度，发现的OTU，以及嵌合体数据；结果如下：

03:39 84Mb 100.0% 5486 OTUs, 9187 chimeras

本条命令的详细使用，请阅读官方文档 http://www.drive5.com/usearch/manual/cmd_cluster_otus.html

小技巧：统计fasta文件中序列的数量

fasta文件每条序列以大于号(>)开始，其数量与序列数量相同，使用grep检索含有>的行，同时用-c参数对数量进行统计，即可快速获得fasta文件中序列数量。

# 查看OTU数量 grep '>' -c temp/otus.fa

写在后面

今天先到这里，本文已经讲了三个程序的使用，够大家学习一会的了。要想了解这些程序的更多功能，一定要阅读程序的帮助全文，才能有更深入的理解。

下节预告：扩增子分析解读4去嵌合体,生成OTU表和代表性序列

(宏基因组7月文章目录，更多精彩等你读)

Reference

1. http://www.drive5.com/usearch
2. http://www.drive5.com/usearch/manual/cmds_all.html
3. http://www.drive5.com/usearch/manual/cmd_fastx_uniques.html
4. http://www.drive5.com/usearch/manual/cmd_cluster_otus.html
5. http://drive5.com/usearch/manual/chimera_formation.html
6. http://www.cnblogs.com/xudongliang/p/6497465.html