Archives

快速对扩增子序列进行数据拆分: seqtk_demultiplex

扩增子测序常用的模式是在双端引物两边各加上几个碱基的条形码(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

Comments are closed.