扩增子分析解读2提取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://qiime.org/scripts/extract_barcodes.html

 

5. 质控及样品拆分

上步对序列方向进行了调整所有为正向，并切开了barcode与扩增序列。下面使用split_libraries_fastq.py对混池根据barcode拆分样品，同时筛选质量大于20(即准确度99%)的序列进行下游分析。参数解释以下：

-i 输入序列文件，上步产生；

-b 输入barcode文件，上步产生；

-m 实验设计，依赖样品barcode列表将序列按样品从新命名

-q 测序质量阈值，20为99%准确率，30为99.9%准确

--max_bad_run_length 容许连续的低质量碱基调用的最大值，默认值为3

--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://qiime.org/scripts/split_libraries_fastq.html

 

6. 切除引物序列

这里咱们介绍一款高通量引物切除软件，cutadapt，2017-6-16最新版1.14；

https://pypi.python.org/pypi/cutadapt

此软件2011年发布至今一直在更新，Google Scholar截止17年8月8日引用2263次。

 

Cutadapt 1.14下载和安装

# 下载，请尽可能检查主页下载最新版源码
wget https://pypi.python.org/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 3.6.1
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 73.83 s (58 us/read; 1.04 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
8 1 19.5 1 1
11 1 0.3 1 1
13 1 0.0 1 1
15 9 0.0 2 9
17 42 0.0 2 42
18 202626 0.0 2 202626
19 56 0.0 2 56
20 1 0.0 2 1
21 1 0.0 2 1
32 1 0.0 2 1
38 1 0.0 2 1
39 1 0.0 2 1
41 1 0.0 2 1
309 2 0.0 2 2
310 1 0.0 2 1
311 3 0.0 2 3

=== 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
3 2 19959.9 0 2
7 1 78.0 1 0 1
8 2 19.5 1 1 1
10 6 1.2 1 4 2
11 1 0.3 1 1
12 3 0.1 1 2 1
13 5 0.0 1 3 2
14 24 0.0 2 17 7
15 51 0.0 2 32 14 5
16 71 0.0 2 56 12 3
17 134 0.0 2 92 30 12
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
22 1 0.0 2 1
23 2 0.0 2 2
24 1 0.0 2 1
25 2 0.0 2 2
27 5 0.0 2 5
28 2 0.0 2 2
29 2 0.0 2 2
30 1 0.0 2 1
31 2 0.0 2 2
32 10 0.0 2 10
49 1 0.0 2 1
51 1 0.0 2 1
166 1 0.0 2 1
291 6 0.0 2 6
401 2 0.0 2 2
409 1 0.0 2 1
443 1 0.0 2 1
460 2 0.0 2 2
465 2 0.0 2 2

WARNING: One or more of your adapter sequences may be incomplete. Please see the detailed output above.