Biostack.org

扩增子测序常用的模式是在双端引物两边各加上几个碱基的条形码（barcode）进行区分样本，这样三十到四十个样本可以混在一起进行建库，节省实验成本。

一、seqtk_demultiplex 介绍

程序见 Github seqtk_utils 仓库。

命令行接口

seqtk_demultiplex 因为参数比较多，所以全部采用了选项模式传递命令行参数，接口如下：

$ seqtk_demultiplex
Usage: barcode_demutiplex [options]
Options:
  -1 STR  forward reads file;
  -2 STR  forward reads file;
  -b STR  barcode mapping file;
  -d STR  output filename directory;
  -l INT  minimum barcode length, default: 5；
  -v      print version number.
  -h      display usage.

Example:
barcode_demutiplex  -b barcodes.txt -1 I365.R1.fastq.gz  -2  I365.R2.fastq.gz -l 5  -d  I365

该程序比如接收如下几个参数

-1， 测序结果正向序列， 标准fastq文件，支持gz压缩文件；
-2， 测序结果反向序列， 标准fastq文件，支持gz压缩文件；
-b， barcode的文件；
-d,  输入文件目录；
-l,  barcode 序列长度（如长度大小不一致，填写最短的序列长度）， 默认是5个；

barcode文件格式

采用制表符分隔：三列， 第一列为样本名称，第二列为F端barcode， 第三列是R端barcode

F_1     CAATT   GCTGTGA
F_2     GCAGGA  CGGAAT
F_3     TGAAGG  ATGTCA
F_4     ATGTGCT TAAGA
F_5     CAGAGAT GCTGG
F_6     GCGTT   CATGATA

二、实例及其用法

原始测序数据文件：I365.R1.fastq.gz 和 I365.R2.fastq.gz
barcode 文件： barcodes.txt

seqtk_demultiplex -b barcodes.txt -1 I365.R1.fastq.gz  -2  I365.R2.fastq.gz  -l 5  -d  I365

上面为几种常见执行模式, 支持streaming 操作，也直接支持gz压缩文件, 针对大文件也可以使用pigz进行加速读操作， 比如：

time  seqtk_demultiplex -b barcodes.txt -1 <(pigz -dk -c -p 2 I365.R1.fastq.gz) -2 <(pigz -dk -c -p 2 I365.R2.fastq.gz)  -l 5  -d  I365

real    0m1.546s
user    0m3.368s
sys     0m0.848s

速度很快。

三、其它工具

先介绍fastq-multx这个工具

下载安装

wget  https://github.com/brwnj/fastq-multx/archive/v1.3.1.tar.gz
tar xzvf fastq-multx-1.3.1.tar.gz
cd fastq-multx-1.3.1
make

命令行接口:

fastq-multx
Usage: fastq-multx [-g|-l|-B] <barcodes.fil> <read1.fq> -o r1.%.fq [mate.fq -o r2.%.fq] ...
Version: 1.3.1

Output files must contain a '%' sign which is replaced with the barcode id in the barcodes file.
Output file can be n/a to discard the corresponding data (use this for the barcode read)

Barcodes file (-B) looks like this:

<id1> <sequence1>
<id2> <sequence2> ...

Default is to guess the -bol or -eol based on clear stats.

If -g is used, then it's parameter is an index lane, and frequently occuring sequences are used.

If -l is used then all barcodes in the file are tried, and the *group* with the *most* matches is chosen.

Grouped barcodes file (-l or -L) looks like this:

<id1> <sequence1> <group1>
<id1> <sequence1> <group1>
<id2> <sequence2> <group2>...

Mated reads, if supplied, are kept in-sync

Options:

-o FIL1     Output files (one per input, required)
-g SEQFIL   Determine barcodes from the indexed read SEQFIL
-l BCFIL    Determine barcodes from any read, using BCFIL as a master list
-L BCFIL    Determine barcodes from <read1.fq>, using BCFIL as a master list
-B BCFIL    Use barcodes from BCFIL, no determination step, codes in <read1.fq>
-H          Use barcodes from illumina's header, instead of a read
-b          Force beginning of line (5') for barcode matching
-e          Force end of line (3') for barcode matching
-t NUM      Divide threshold for auto-determine by factor NUM (1), > 1 = more sensitive
-G NAME     Use group(s) matching NAME only
-x          Don't trim barcodes off before writing out destination
-n          Don't execute, just print likely barcode list
-v C        Verify that mated id's match up to character C (Use ' ' for illumina)
-m N        Allow up to N mismatches, as long as they are unique (1)
-d N        Require a minimum distance of N between the best and next best (2)
-q N        Require a minimum phred quality of N to accept a barcode base (0)

命令行接口比较复杂，这里只解决主要使用方法：
主要参数：

-o， 输出文件，一个输入文件一个输出文件流，格式： %.R1.fq.gz， %为barcode对应的样本字符串，
     如有多个文件，需要多个 -o 参数;
-m， barcod允许的主要错配个数，一般设置为0， 默认设置了1；
-B,  barcode文件，允许单端和双端barcode;
-n,  打印barcode序列；
-b,  从序列的5'端碱基开始匹配barcode;
-e， 从序列3'端开始匹配序列；

其它参数主要是控制匹配参数， 比如

-q, 控制barcode碱基的最小phred quality值，默认为0， 不控制；
-d 控制匹配的最佳barcode位置和此佳barcode位置的位置， 两个匹配距离不能超过2个碱基；

barcode文件格式：

单端barcode 只需要提供两列数据，双端barcode需要中间加上’-‘, 示例如下：

F_1 CAATT-GCTGTGA
F_2 GCAGGA-CGGAAT
F_3 TGAAGG-ATGTCA
F_4 ATGTGCT-TAAGA
F_5 CAGAGAT-GCTGG
F_6 GCGTT-CATGATA

使用情况， 参考：

$time fastq-multx   -B mapping_file.txt  -m 0  -b   I365.R1.fastq.gz  I365.R2.fastq.gz    -o   I365/%.R1.fq.gz  -o I365/%.R2.fq.gz  >/dev/null
Using Barcode File: mapping_file.txt
End used: start

real    0m7.107s
user    0m14.406s
sys     0m0.778s

性能很不错， 但是！！！， R端序列的barcode没有被切除， 留待Hack.

四、注意事项

1. 本工具常用使用场景为对扩增子数据进行拆分， 设计模式很简单，就是使用Hash Map去匹配，所以速度绝对没有问题。
2. 数据拆分不允许barcode出现错配；
3. barcode识别从第一个碱基开始；
4. 支持barcode序列长度不一致问题，使用最短的序列长度设置 -l 参数。
5. 不支持gz压缩文件输出格式

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料， 版权归 上海逻捷信息科技有限公司 所有

Last Upate: 2017-09-28 7:05 PM