pythonbarcode碱基错配_扩增子分析解读2提取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

lengthcountexpectmax.errerror counts

3319959.903

444990.004

62311.902

8119.511

1110.311

1310.011

1590.029

17420.0242

182026260.02202626

19560.0256

2010.021

2110.021

3210.021

3810.021

3910.021

4110.021

30920.022

31010.021

31130.023

=== 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

lengthcountexpectmax.errerror counts

3219959.902

7178.010 1

8219.511 1

1061.214 2

1110.311

1230.112 1

1350.013 2

14240.0217 7

15510.0232 14 5

16710.0256 12 3

171340.0292 30 12

183270.02189 117 21

1910591750.021056863 2069 243

20138460.021817 10955 1074

217440.025 10 729

2210.021

2320.022

2410.021

2520.022

2750.025

2820.022

2920.022

3010.021

3120.022

32100.0210

4910.021

5110.021

16610.021

29160.026

40120.022

40910.021

44310.021

46020.022

46520.022

WARNING:

One or more of your adapter sequences may be incomplete.

Please see the detailed output above.