【原】扩增子分析解读2提取barcode,质控及样品拆分,切除扩增引物

宏基因组 2020-10-09

展开全文

写在前面

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

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

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

16S扩增子测序数据主要来自HiSeq2500产出的双端各250 bp (PE250)数据，因为读长长且价格便宜(性价比高)。HiSeqX PE150和MiSeq PE300也比较常见，但PE150过短分辨率低，而PE300价格高且末端序列质量过低；此外454在之前研究较多，但现在设备已经停产；PacBio读长长可直接测序16S全长1.5kb代表未来的趋势。

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

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

第二节. 提取barcode,样品拆分及质控,切除扩增引物

本节课程，需要完成扩增子分析解读1质控,实验设计,双端序列合并

先看一下扩增子分析的整体流程，从下向上逐层分析。

分析前准备

# 进入工作目录
cd example_PE250

上一节回顾：我们拿到了双端数据，进行了质控、并对实验设计进行了填写和检查无误、最后将双端数据合并为单个文件进行下游分析。

接下来我们将序列末端的barcode标签切下来，因为它们是人为添加的，不属于实验对象；再根据标签序列与实验设计文件比对，对每条序列属于哪个样品进行分类；最后我们切除掉扩增使用的引物，因为它们是人工合成的相似序列，并不是物种的真实序列。这样我们就获得了高质量的扩增区域数据，并且序列名称中包括了样品信息。

4. 提取barcode

为什么扩增子有barcode？
我以前做过sRNA-Seq, RNA-Seq, ChIP-Seq等等，都是一个文库对应一个样品，为什么没有用过barcode拆分。
原因是扩增子目前研究对象细菌、真菌多样性没有表达基因数量大，一般是几百到千的水平，对数据量要求最多10万条序列即可饱合。而Illumina测序仪的通量很高，采用Index来区分每个文库，仍然每个文库的数据量达到千万的级别，加上建库测序的成本也不会低于千元。对于扩增子动辄成百上千的样品即太贵，又浪费。因此将扩增子样本添加上barcode(标签)，通常将48/60个样品混合在一起，构建一个测序文库，达到高通量测序大量样品同时降低实验成本的目的。

通常的测序仪下机数据，只经过Index比对，拆分成来自不同文库的数据文件，分发给用户。而扩增子的一个文库包括几十个样品，还区要通过每个样品上标记的特异Barcode进一步区分，再进行下游分析。

注：如果你是自己构建测序文库，返回数据可以按下文拆分样品；如果是公司建库并测序，他们有样品的barcode信息，通常会返给用户样品拆分后的数据，无需本节中的操作。但原理还是要理解，对整体分析的把握非常重要。

Barcode在扩增子的位置和类型？

Barcode位于引物的外侧，比较典型的有三种，上图展示的为最常用的barcode位于左端(正向引物上游)，此外还有右端和双端两类也比较常用。
本分析采用的数据类型为右端barcode。

extract_barcodes.py是QIIME中用于切除barcode的脚本，支持你想到的所有类型。
-f参数为输入文件，即上文中合并双端数据后的文件；
-m为实验设计文件；
-o为输出切下的barcode的数据目录；
-c为barcode类型选择，目前的barcode_paired_stitched选项支持右端和双端类型，如果是左端类型，请改为barcode_single_end；
—bc1_len设置左端barcode的长度，按实验设计添写，本数据只有右端则为零；
—bc2_len设置右端barcode的长度，按实验设计添写，本数据右端长度为6bp，常用的还有8，10；
-a是使用实验设计中的引物来调整序列的方向，本实验的测序无方向，必须开此选项调整，有些公司的建库类型有方向，但开了总没错,顶多多花点计算时间；
—rev_comp_bc2是将右端barcode取反向互补再与左端相连，建议打开，这样你的实验设计中barcode一栏直接将添写原始barcode即可，不用再考虑方向了；如果是双端则将两个barcode直接连起来填在barcode列，不用考虑方向。

# 提取barcode
extract_barcodes.py -f temp/PE250_join/fastqjoin.join.fastq     -m mappingfile.txt     -o temp/PE250_barcode     -c barcode_paired_stitched --bc1_len 0 --bc2_len 6 -a --rev_comp_bc2

这步任务会在输出目录temp/PE250_barcode中生成5个文件

barcodes.fastq # 切下来的barcode，用于后续拆分样品
barcodes_not_oriented.fastq # 方向不确定序列的barcode。连引物都不匹配，质量太差，建议不再使用
reads1_not_oriented.fastq # 方向不确定序列的序列，可能barcode切错方向。连引物都不匹配，质量太差，不建议使用
reads2_not_oriented.fastq # 空文件
reads.fastq # 序列文件，与barcode对应，用于下游分析

更多说明建议阅读帮助 http:///scripts/extract_barcodes.html

5. 质控及样品拆分

上步对序列方向进行了调整全部为正向，并切开了barcode与扩增序列。下面使用split_libraries_fastq.py对混池根据barcode拆分样品，同时筛选质量大于20(即准确度99%)的序列进行下游分析。参数解释如下：
-i 输入序列文件，上步产生；
-b 输入barcode文件，上步产生；
-m 实验设计，依赖样品barcode列表将序列按样品重新命名
-q 测序质量阈值，20为99%准确率，30为99.9%准确
—min_per_read_length_fraction 最小高质量序列比例，0.75即最少有连序75%的碱基质量高于20
—max_barcode_errors barcode允许的碱基错配个数，建议设置0为不允许，即使修改为1,2通常也不允许错配，不信你试试
barcode_type调置barcode类型，默认为golay_12，并支持错配；我们通常设置为整数，对应barcode的长度总和，本实验中添0+6=6。

# 质控及样品拆分
split_libraries_fastq.py -i temp/PE250_barcode/reads.fastq     -b temp/PE250_barcode/barcodes.fastq     -m mappingfile.txt     -o temp/PE250_split/     -q 20 --max_bad_run_length 3 --min_per_read_length_fraction 0.75 --max_barcode_errors 0 --barcode_type 6

运行结果会输出三个文件

histograms.txt # 所有序列长度分布数据，可知长分范围308-488,峰值为408
seqs.fna # 质控并拆分后的数据，序列按样品编号为SampleID_0/1/2/3
split_library_log.txt # 日志文件，有基本统计信息和每个样品的数据量；查看可知每个样品最大数据量为110454，最小值为10189

更多说明建议阅读帮助 http:///scripts/split_libraries_fastq.html

6. 切除引物序列

这里我们介绍一款高通量引物切除软件，cutadapt，2017-6-16最新版1.14；
https://pypi./pypi/cutadapt
此软件发2011年发布至今一直在更新，Google Scholar截止17年8月8日引用2263次。

Cutadapt 1.14下载和安装

# 下载，请尽量检查主页下载最新版源码
wget https://pypi./packages/16/e3/06b45eea35359833e7c6fac824b604f1551c2fc7ba0f2bd318d8dd883eb9/cutadapt-1.14.tar.gz
# 解压
tar xvzf cutadapt-1.14.tar.gz
# 进入程序目录
cd cutadapt-1.14/
# 安装在当前用户目录，不需管理员权限
python setup.py install --user

cutadapt切除双端引物及长度控制，参数详解：
-g 5’端引物
-a 3’端引物，需要将引物取反向互补
-e 引物匹配允许错误率，我调置0.15，一般引物20bp长允许3个错配，为了尽量把引物切干净
-m 最小序列长度，根据情况设置，本实验扩增V5-V7区，长度主要位于350-400，故去除长度小于300bp的序列，无经验可不填此参数
—discard-untrimmed 如果引物末切掉的序列直接丢弃，防止了现假OTU
seqs.fna 为输入文件
PE250_P5.fa 为输出文件

cutadapt -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -e 0.15 -m 300 --discard-untrimmed temp/PE250_split/seqs.fna -o temp/PE250_P5.fa

程序运行结果会在屏幕上输出结果统计报告，如去接头比例、去掉过短序列比例等；还有去除引物的长度分布，看看有没有异常。如3’端序列没有取反向互补会无法去除19bp序列，而是几bp的错误序列。

Cutadapt结果报告

This is cutadapt 1.14 with Python 2.7.12
Command line parameters: -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -e 0.15 -m 300 --discard-untrimmed temp/PE250_split/seqs.fna -o temp/PE250_P5.fa
Trimming 2 adapters with at most 15.0% errors in single-end mode ...
Finished in 41.03 s (32 us/read; 1.87 M reads/minute).

=== Summary ===

Total reads processed:               1,277,436
Reads with adapters:                 1,277,194 (100.0%)
Reads that were too short:               8,849 (0.7%)
Reads written (passing filters):     1,268,345 (99.3%)

Total basepairs processed:   522,379,897 bp
Total written (filtered):    495,607,409 bp (94.9%)

=== Adapter 1 ===

Sequence: GGAAGGTGGGGATGACGT; Type: regular 3'; Length: 18; Trimmed: 202757 times.

No. of allowed errors:
0-5 bp: 0; 6-12 bp: 1; 13-18 bp: 2

Bases preceding removed adapters:
  A: 96.3%
  C: 1.5%
  G: 0.8%
  T: 1.3%
  none/other: 0.0%
WARNING:
    The adapter is preceded by "A" extremely often.
    The provided adapter sequence may be incomplete.
    To fix the problem, add "A" to the beginning of the adapter sequence.

Overview of removed sequences
length    count    expect    max.err    error counts
3    3    19959.9    0    3
4    4    4990.0    0    4
6    2    311.9    0    2
18    202626    0.0    2    202626
19    56    0.0    2    56
=== Adapter 2 ===

Sequence: AACMGGATTAGATACCCKG; Type: regular 5'; Length: 19; Trimmed: 1074437 times.

No. of allowed errors:
0-5 bp: 0; 6-12 bp: 1; 13-19 bp: 2

Overview of removed sequences
length    count    expect    max.err    error counts
18    327    0.0    2    189 117 21
19    1059175    0.0    2    1056863 2069 243
20    13846    0.0    2    1817 10955 1074
21    744    0.0    2    5 10 729
WARNING:
    One or more of your adapter sequences may be incomplete.
    Please see the detailed output above.