微生物组16S rRNA数据分析小结: vsearch

本次笔记内容：

• vsearch超简介绍
• install vsearch (linux)
• 常用command

vsearch超简介绍

usearch的可靠替代，代码开源，没有数据量大小限制，速度更快，免费。
其文档vsearch_manual.pdf见https://github.com/torognes/vsearch

install vsearch (linux)

git clone https://github.com/torognes/vsearch.git
cd vsearch
./autogen.sh
./configure
make
sudo make install   # as root or sudo make install

常用command

--fastx_filter

trim或者filter fasta/fastq格式的序列文件。下列命令中的seqs.fastq为一个已经割库后，包含了所有samples的大fastq文件。可以由qiime1的multiple_split_libraries.py得到。当然也可以用vsearch来jion双端序列再split，这里不详细描述。
--fastq_maxee: 根据质量来filter, 例如expected error大于1.0则去除
--fastq_trunclen: 240，将序列修剪到240的长度，短于240则去除
--fasta_width: 0, 默认设置为80，设置output的fasta中，每行有多少个核苷酸。设置为0则将核苷酸排列为一行。

vsearch --fastx_filter seqs.fastq \
        --fastq_maxee 1.0 --fastq_trunclen 240 \
        --fastaout filtered.fa \
        --fasta_width 0

--derep_fulllength

查找重复序列，并将其归拢到一起。output显示unique的序列及这些unique序列出现的次数。与usearch的 fastx_uniques类似。
--minuniquesize: 在去重之后统计各unique序列的丰度，小于该值则去除
--sizeout: add abundance annotation to ouput fasta file
--uc: 生成一个tab分隔的表格，是每条序列的具体重复情况。第一列是S(cluster centroid,即重复序列中的一个中心序列)，H(hit assigned to the cluster)，C(cluster records)构成的分类变量列。具体每一列代表什么，见vsearch_manual.pdf.

vsearch --derep_fulllength filtered.fa \
        --minuniquesize 2 \
        --sizeout --fasta_width 0 \
        --output derep.fa \
        --uc derep.uc

precluster(预聚类) --cluster_size

后续会基于这一步产生的all.precluster.fa去嵌合体

vsearch --cluster_size derep.fa \
        --id 0.97 \
        --sizeout --fasta_width 0 \
        --uc all.precluster.uc \
        --centroids all.precluster.fa

去嵌合体： 通过denovo和参考数据库两种方式去除

自行下载rdp_gold.fa, 嵌合体相关的数据库。经过这两步得到all.ref.nonchimeras.fa为已经去除了嵌合体，基于预聚类的序列。

vsearch --uchime_deno all.precluster.fa \
        --sizein --sizeout --fasta_width 0 \
        --nonchimeras all.denovo.nonchimeras.fa

vsearch --uchime_ref all.denovo.nonchimeras.fa \
         --db rdp_gold.fa \
         --sizein --sizeout --fasta_width 0 \
         --nonchimeras all.ref.nonchimeras.fa

将去除嵌合体的信息mapping回原始数据上

由于上一步去除嵌合体是基于预聚类的序列，这一步将去除嵌合体的信息mapping回去，然后再正式cluster. vsearch这样做应该是为了缩短时间，避免文件过大等问题。（？有待确定）

perl map.pl derep.fa all.precluster.uc all.ref.nonchimeras.fa > all.nonchimeras.derep.fa
# perl map.pl filtered_uparsed.fa all.derep.fa all.nonchimeras.derep.fa > all.nonchimeras.fa
perl map.pl derep.fa all.derep.uc all.nonchimeras.derep.fa > all.nonchimeras.fasta

以97%相似程度聚类：得到代表序列和OTU table

注意这里input的all.nonchimeras.fa是在上一步中，将嵌合体信息mapping回去，得到的处理后序列。all.otu.fa为代表序列。otu_raw.tab为原始OTU table.

vsearch --cluster_size all.nonchimeras.fa \
        --id 0.97 \
        --sizein --sizeout --fasta_width 0 \
        --uc all.cluster.uc \
        --relabel OTU \
        --centroids all.otus.fa

vsearch --usearch_global filtered_uparsed.fa \
         --db all.otus.fa \
         --id 0.97 \
         --uc map.uc \
         --otutabout otu_raw.tab

基于OTU table和代表序列的后续分析

可以参考微生物组16S rRNA数据分析小结：从fastq测序数据到OTU table的usearch部分，得到物种注释结果，alpha及beta的初步分析。

附：
嵌合体数据库rdp_gold.fa wget http://drive5.com/uchime/rdp_gold.fa
map.pl脚本：https://gist.github.com/GeoMicroSoares/29ca37f321d03b0e2325d4ed80c81fd5
...因为时间仓促，没有整理的很周全。以后如果工作中需要会慢慢补上orz