Archives

Categories

行计数 和 (组合)列 uniq 并统计计数工具: tabtk_nlines 和 tabtk_uniq

本文介绍两个小工具:tabtk_nlines 和 tabtk_uniq

1. tabtk_nlines

其实使用 wc -l 或者 grep -P '^#' | wc -l 可以很方便的解决行计数问题,但是处理数据时候需要做两个事情:

  • 添加样本信息
  • 过滤掉注释行

使用 tabtk_nlines 可以干净的解决这个问题:

命令行接口:

$ tabtk_nlines
Usage: tabtk_nlines [options] <text> <samples>
Options:
  -i  flag to ignore lines start with '#';
  -v  print version number

需要两个参数, 一个文本文件(支持压缩包)和样本信息

示例:

统计F1 样本 F1-ko.tsv.gz 注释到ko的基因个数,有表头信息,以#开头为注释信息

$ tabtk_nlines  -i  F1-ko.tsv.gz  F1
F1      7369
2. tabtk_uniq

tabtk_uniq 程序为了解决组合列进行Uniq操作的实现, 并可以通过 -c 可选项指示是否显示频数, -f 选项指定fields 信息

命令行接口:

$ tabtk_uniq
Usage: tabtk_uniq [options] <text>
Options:
  -f STR fields pattern: 1:2:3, [1];
  -c     display the counts;
  -i     flag to ignore lines start with '#';
  -v     print version number

示例:

测试文件 file.txt

#1      1       1
1       2       4
3       4       5
1       2       5
3       4       10
4       8       2
1       2       7
2       5       19

统计第一列和第二列组合Uniq

$tabtk_uniq  -i  -f1:2 file.txt
1       2
4       8
3       4
2       5

统计第一列和第二列组合Uniq并显示频数

$tabtk_uniq  -i -c   -f1:2 file.txt
1       2       3
4       8       1
3       4       2
2       5       1

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

Last Upate: 2017/10/16 21:20

Mask和过滤低复杂度序列:seqtk_sdust

处理 small rna 病毒组数据组装时, 经过multi-kmer, multi-coverage 调试最优参数,进行组装后,发现很多序列具有很低的复杂度, 使 在用 mdust 可以mask掉序列中低复杂度序列,但是还需要进一步过滤,不够优雅。

先前看到过Heng Li的一篇博客(A reimplementation of symmetric DUST )提到重新实了 SDUST 算法,并在 minimap2 找到了实现文件,sdust.c 和 sdust.h, 经过简单处理, 实现了支持根据复杂度区域覆盖度 -c 可选项进行序列过滤,默认设置 2.0 , 即不进行序列过滤。

程序在 Github seqtk_utils 仓库。

1 编译
gcc -O2   -D_NO_NT4_TBL   seqtk_sdust.c   sdust.c  -o  seqtk_sdust  -lz  kalloc.c
2 命令行接口
$ seqtk_sdust
    Usage: seqtk_sdust [-p 2.0] [-w 64] [-t 20] <in.fa>
    Options:
      -c FLOAT  coverage cutoff for different contigs/reads to filter, defalut [2.0]
      -w INT    window size W, defalut [64]
      -t INT    score threshold T, defalut [20]
      -v        print version number

保留了先前参数,添加了 -c 参数;

3 示例
seqtk_sdust  -c  0.9   -w 20  -t 20  test.fasta  >/dev/null

过滤掉低复杂区覆盖度超过90%的序列;

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

Last Upate: 2017-10-16 6:35 PM

使用 libxlsxwriter 实现将tsv文件转换成xlsx文件: tabtk_xlsx

这篇为tsv-utils系列的一个小工具,目的是将一堆制表符分隔的文本文件转换成Excel文件(xlsx格式,不是xls),转换成 xlsx 后体积会变小不少。

先前使用Python也写了一个版本, 性能不行,就直接使用 C库 libxlsxwriter实现了C版本, 程序见 Github tsv-utils 仓库 。

1 libxlsxwriter 介绍

libxlsxwriter 提供了很简洁的创建Excel文件的函数接口,100% 兼容xlsx, 官方文档 介绍的很详细,有很多实例, 接口都很简单,可以创建各种图形,是自动化生成Excel文档的不二选择。

libxlsxwriter

下面是 libxlsxwriter 的其它语言接口:

2 tabtk_xlsx 命令行接口

tabtk_xlsx 命令行接口如下:

$ tabtk_xlsx
Usage: tabtk_xlsx  <xlsx> <sheet:file_path> [<sheet:file_path>]
version: 0.0.1

接口很简单,变量第一个位置为输入xlsx文件名称, 后面为需要转换的文件列表,格式为 sheet名称:文件路径名称

3 tabtk_xlsx 测试和实例

首先测试下Python版本和C版本性能差异:

对NCBI Taxonomy文件进行转换, 1560796行。

$time  tabtk_xlsx.py  ncbi.map.xlsx  ncbi:ncbi.map
real    2m47.532s
user    2m46.537s
sys     0m0.829s

$time  tabtk_xlsx  ncbi.map.xlsx  ncbi:ncbi.map
real    0m12.577s
user    0m11.980s
sys     0m0.582s

善于使用工具可以节省时间,当然也可以节省金钱。

对16S的不同水平的物种分类进行合并生成 xlsx 文件多个表文件示例:

tabtk_xlsx  taxonomy.xlsx  phylum:phylum.txt order:order.txt  genus:genus.txt species:species.txt

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

Last Upate: 2017-10-12 2:52 PM

CentOS 7 手动设置静态IP的方法

这里介绍两种方法设置CentOS 7 静态IP地址的方法,力求图文并茂。

方法1: 使用 GNOME 桌面网络管理窗口

第一步: 从管理界面进入设置界面

第一步

第二步: 选择网络管理

第二步

第三步: 选择设备, 一般选择em1, 根据插入的网口决定

第三步

第四步: 设置自动连接选项

第四步

第五步: 设置IP地址、网关、DNS地址

第五步

设置点击应用后应该就生效了。

方法2: 修改文件方法:

进入网络脚本目录

# cd /etc/sysconfig/network-scripts
# vim ifcfg-em1

修改内容如下:

TYPE=Ethernet
BOOTPROTO=none
DEFROUTE=yes
IPV4_FAILURE_FATAL=no
IPV6INIT=yes
IPV6_AUTOCONF=yes
IPV6_DEFROUTE=yes
IPV6_FAILURE_FATAL=no
NAME=em1
UUID=d2325ba4-d9f5-4363-bb25-d9a7ddbbd991
ONBOOT=yes
DNS1=192.168.0.1
IPADDR=192.168.0.2
PREFIX=24
GATEWAY=192.168.0.1
IPV6_PEERDNS=yes
IPV6_PEERROUTES=yes

根据实际情况,修改ONBOOT、DNS1、IPADDR、PREFIX和GATEWAY 参数值, 然后重启网络:

# systemctl restart network

设置IP就生效了

注: UUID冲突的解决方案

由于各种原因, /etc/sysconfig/network-scripts 目录的UUID可能会变化,如果出现:

# service network restart
Restarting network (via systemctl):  Job for network.service failed because the control process exited with error code. \
See "systemctl status network.service" and "journalctl -xe" for details.

错误信息, 根据:

# systemctl status network.service
# journalctl -xe

出现这样的信息:

ifcfg-rh: cannot load /etc/sysconfig/network-scripts/ifcfg-em1 due to \

conflicting UUID for /etc/sysconfig/network-scripts/ifcfg

需要使用:

nmcli con

查看网卡的UUDI,并进行修改。

如何并行下载NCBI拼装数据库的Genbank和Refseq数据

NCBI 拼装数据库目前(时间:2017-10-09)已经超过了 10万个基因组数据, Genbank拼装版本已经超过了12万, 下面介绍如何下载这些数据,很多 metagenome reads 序列分类软件都需要这些参考基因组数据, 下载流程主要通过NCBI ASSEMBLY Report 获取基因组数据在FTP的位置.

前面介绍过 Aspera: 数据高速上传下载的利器:Aspera及其在生命科学中的应用,为什么还要使用wget下载呢, Aspera下载小文件其实不是很划算,建立连接很花时间,还没有wget快, 当然可以直接下载整个ASSEMBLY目录, 但是里面有很多个头比较大的东西我们分根本不需要, 根据实践经验, 使用NCBI ASSEMBLY Report文件作为索引下载需要的文件是比较好的选择.

1 获取NCBI的Assembly数据库目录结构

grep -P -v  ^"#" assembly_summary_refseq.txt | cut -f20  | grep -w -v  "na"   >gcf.txt
grep -P -v  ^"#" assembly_summary_genbank.txt | cut -f20  | grep -w -v  "na"  >gca.txt

2 下载数据

genome-download  gcf.txt gff  GCF  2>log.txt&
genome-download  gca.txt gff  GCF  2>log.txt&

通过使用 log.txt 可以查看失败的下载条目. 可以重新调用genome-download下载, 通过这两步可以轻松下载NCBI的基因组拼装数据.

genome-download 可以通过Github 仓库: ncbi-utils 下载.

后面会介绍如何用好这些下载的数据, 比如: 菌株鉴定,序列分类等.

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

Last Upate: 2017-10-09 11:54 AM

CentOS 7 “The system network services are not compatible with this version” 解决方案

CentOS7 GNOME 桌面系统网络设置如果遇到:

The system network services are not compatible with this version

这个错误提示,应该是:NetworkManager 服务没有启动, 最简单解决就是启动NetworkManager服务,并设置开机自动启动, 命令执行:

# service NetworkManager start
# chkconfig NetworkManager on

Linux 平台轻松解决Python3环境: conda

一、 引子(需求场景)

因为 Python2 和 Python3 之间的语法差异, 现在Python开发者主要有两大阵营,Python 2.7.x 版本阵营 和 Python 3 阵营, 为了使用不同阵营的包, 导致我们要维护两套Python系统。 当前Linux发行版本主要还是Python2环境,所以需要额外一套方案使得Python2和Python3共存, 最简单的方案应该就是切换环境变量,这样做的好处是我们可以包括任意多的版本,这里使用conda来切换环境。

这篇可以作为介绍Bioconda的引子,只介绍如何安装 miniconda 解决 Python 环境的问题,如何使用bioconda 安装生物信息软件,后面再讲。

bioconda

二、使用Conda进行Python包管理

下面就讲几点:

  • miniconda 介绍
  • miniconda 安装
  • 设置环境变量
  • 安装python包
  • 设定和管理特定Python环境

1 miniconda 介绍:

先摆出几个概念,Anaconda, conda, Miniconda, pip,可以参考官方文档:
What’s the difference between Anaconda, conda, and Miniconda?
也可以参考一篇博客:Anaconda、Miniconda、Conda、pip的相互关系 , 讲的很清楚, 为了完整性,这里简单介绍:

  • conda, 包管理器,安装和管理包的命令接口;
  • Miniconda, conda的安装程序,会安装Python、conda及其各种依赖程序;
  • Anaconda, Python科学计算生态系统中的完整分布,包含了Conda及其各种库。
  • pip, Python的包管理器.

Windows平台请参考知乎:[CONDA] – 上手环境包管理系统

2 miniconda 安装:

wget  https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
sh  Miniconda3-latest-Linux-x86_64.sh

 >>> ENTER
 >>> yes
 >>> /biostack/bioconda
 >>> no

这里没有将conda 环境变量直接添加到 .bashrc,是为了更方便隔离系统自带的Python及其pip包管理。

3 设置环境变量:

配置和设置 activate 文件

cd /biostack/bioconda/envs
echo 'export PATH=/biostack/bioconda/bin:$PATH' >activate.sh
vim ~/.bashrc

在.bashrc设置alias

alias bioconda='source /biostack/bioconda/envs/activate.sh'
alias pip='pip --trusted-host pypi.douban.com '

立即生效

$source ~/.bashrc

设置切换 Python3 环境模式:

$bioconda
$which pip
    alias pip='pip --trusted-host pypi.douban.com '
    /biostack/bioconda/bin/pip

通过这样设置,我们可以在需要的时候才切换conda环境,执行 ·bioconda· 之前不会对系统产生任何影响;

4 安装包:

现在我们有两种模式安装Python包,pip和conda。
安装模式都很简单, 而且二者兼容,可以二选一。

pip 模式(可以参考:Pythoner的福利,豆瓣的PyPI源 设置douban的python包镜像):

pip install scipy

conda 安装模式

conda install numpy

conda 直接访问国外网站速度比较慢,需要设置国内镜像:
执行:

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --set show_channel_urls yes

也可以手动添加到 ~/.condarc, 内容如下:

channels:
  - bioconda
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - defaults
show_channel_urls: true

这样我们就可以使用命令:

conda install matplotlib

利用 tuna 的镜像安装 python 包;

安装路径,其实安装的Python包都在python发行版本根目录:

$ python -c "import scipy; print (scipy.__file__)"
/biostack/bioconda/lib/python3.6/site-packages/scipy/__init__.py

使用conda list 和 pip list 都可以获取已安装的包;

conda安装模式可以指定仓库(channel), 比如安装bioconda的软件:

conda config --add channels bioconda
conda install -c  bioconda  agg

5 设定特定Python环境:

需求:比如我们需要安装2.7的Python环境:

通过conda create进行安装py-default环境,设置python版本 2.7 (实际安装版本: 2.7.13)

$conda create --name py-default   python=2.7

激活py-default环境

$source activate py-default

查看python版本

$which python
/biostack/bioconda/envs/py-default/bin/python

安装包,并安装到py-default环境

$conda  install -n  py-default docopt

查看安装路径:

$python -c "import docopt; print (docopt.__file__)"
/biostack/bioconda/envs/py-default/lib/python2.7/site-packages/docopt.pyc

在py-default环境生效前提下可以不需要指定-n py-default

$conda  install gffutils
/biostack/bioconda/envs/py-default/lib/python2.7/site-packages/gffutils/__init__.pyc

如果不在py-default环境下指定-n py-default docopt 也可以正确安装。

退出环境变量:

source deactivate py-default

列出所有的环境:

$conda-env list
# conda environments:
#
py-default               /biostack/bioconda/envs/py-default
root                  *  /biostack/bioconda

删除Python环境:

conda-env  remove -n py-default

这样在/biostack/bioconda/envs/py-default 构建的py-default变量就会被删除了。

三、后续

基于conda的bioconda软件分发已经占领了R和Python包管理和安装,很多其它生物信息软件也陆续被加入bioconda阵营,值得深入研究,后续会研究:

  • 如何构建channel镜像;
  • 通过 anaconda 构建 conda源。

对扩增子序列进行引物匹配: usearch 和 primersearch

primer
经典问题场景:
数据:16S 扩增子序列,已完成双端合并,amplicon.fasta 文件, 一对引物序列:primers.txt

我们想做:

  1. 需要知道引物是否匹配,过滤掉一些噪音序列;
  2. 多对引物在同一个文库,需要根据引物把目标序列提取出来;

下面介绍两个比较简单的引物识别程序(双端引物序列模式)

一、USEARCH

强大的数据处理工具集合,但是: 只能免费获取32位版本,最大内存使用4G, 大部分的数据分析任务还是可以应付。

下载地址:32Bit , 选择需要下载的版本,地址会通过邮件进行发送。

如果想购买64bit 版本可以移步:buy64bit , 也可通过代理商购买,现在的销售模式是大版本免费升级,当前大版本为10.0

这里只涉及到引物匹配子命令:search_pcr

命令行接口:

-db, 引物序列文件, fasta格式;
-strand, 支持 plus 和 both 两种模式
-maxdiffs, 单端引物去最大差异
-minamp,匹配区域的最小长度;
-maxamp,匹配区域的最大长度;
-pcrout, 制表符分隔的文件,记录比对的详细信息;
-threads,  线程数;
-ampout,    匹配区域的序列文件;
-log, 输出log 信息;

-pcr_strip_primers 可以切除引物区域,如果引物中存在简并碱基,可能程序会报错,使用需要谨慎。

cat primers.txt 

>fw
GTGCCAGCMGCCGCGG
>rw
CCGTCAATTCMTTTRAGTTT

现在执行:

usearch10.0.240_i86linux32  -search_pcr amplicon.fasta   -db  primers.txt -strand both -threads 5 \
-maxdiffs 2  -minamp 390  -maxamp 430  --pcrout amplicon-pcrout.txt -log  amplicon.log            \
-ampout  amplicon-pcrout.fasta;

现在:amplicon-pcrout.fasta就是我们需要的文件,对于16S序列,可能存在正向和方向序列共存的情况,以及正确切除引物的需求, 这个问题留到 seqtk_primer_strip 介绍。

pcrout 输出文件格式见: pcrout, 后面对序列的匹配情况进行过滤需要使用这个文件,因为同一个序列可能在不同位置都匹配,这样都会report结果。

二、EMBOSS

EMBOSS 是老牌的数据处理工具箱,内容很丰度,小工具相当多, 达到300个之多,今天需要介绍的是: primersearch 做引物匹配:

primersearch 命令行接口如下:

$ primersearch -h
Search DNA sequences for matches with primer pairs
Version: EMBOSS:6.6.0.0

   Standard (Mandatory) qualifiers:
  [-seqall]            seqall     Nucleotide sequence(s) filename and optional
                                  format, or reference (input USA)
  [-infile]            infile     Primer pairs file
  [-mismatchpercent]   integer    [0] Allowed percent mismatch (Any integer
                                  value)
  [-outfile]           outfile    [*.primersearch] Whitehead primer3_core
                                  program output file

   Additional (Optional) qualifiers: (none)
   Advanced (Unprompted) qualifiers: (none)
   General qualifiers:
   -help               boolean    Report command line options and exit. More
                                  information on associated and general
                                  qualifiers can be found with -help -verbose

主要参数:

-seqall 扩增子序列文件;
-infile 引物文件;
-mismatchpercent 错配百分比; 引物长度为20的话,设置0.5,就是一个错配;
-outfile 输出文件;

引物文件序列格式,制表符分隔:

V4-V5   GTGCCAGCMGCCGCGG        CCGTCAATTCMTTTRAGTTT

程序执行:

primersearch  -infile  primers.txt   -seqall  amplicon.fasta   -outfile  amplicon.txt -mismatchpercent  5

结果有点小麻烦:

Primer name V4-V5
Amplimer 1
        Sequence: 2226

        GTGCCAGC[CAM]GCCGCGG hits forward strand at 1 with 0 mismatches
        CCGTCAATTC[CAM]TTT[GAR]AGTTT hits reverse strand at [1] with 0 mismatches
        Amplimer length: 412 bp

序列的标识符丢掉了,对后续的解析程序会造成很大麻烦,使用小工具primersearch_hack, 可以将序列标识符放在序列注释部分,这样就可以解决问题了:

primersearch_hack   amplicon.fasta  >amplicon_hack.fasta

再重新执行,并查看结果,是不是很愉悦了。

Primer name V4-V5
Amplimer 1
        Sequence: 2226
        HISEQ:593:HW27CBCXY:2:1101:18543:2226
        GTGCCAGC[CAM]GCCGCGG hits forward strand at 1 with 0 mismatches
        CCGTCAATTC[CAM]TTT[GAR]AGTTT hits reverse strand at [1] with 0 mismatches
        Amplimer length: 412 bp

后续还是需要根据匹配的位置,错配等进行过滤,以后再讲。

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

Last Upate: 2017-09-29 4:17 PM

快速对扩增子序列进行数据拆分: 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

对TSV文件进行行筛选:tabtk_subset

继续针对制表符分隔的文件操作进行解读,今天介绍一个工具: abtk_subset ,与先前介绍的一个小工具:seqtk_subseq 功能模式差不多,

如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq

实现方式也差不多, 使用一个HASH SET 存储列表标识符, 然后遍历库文件。

使用 grep 也可以解决TSV文件筛选问题, 但是:

1. grep 无法指定列,求集合操作;
2. 查询列表比较大的时候,效率很低

一、tabtk_subset 介绍

tabtk_subset 命令行接口比较简单, 两个必须参数,一个待搜索库,一个标识符列表;

$tabtk_subset
Usage: tabtk_subset [options] <tsv> <list>
Options:
  -c INT     target col for select, default: 1
  -s INT     subset option: 1 intersection, 0 supplementary set, default: 1
  -v         print vesion number

该程序接收两个可选参数:

-c, 指定列, 默认为1;
-s, 子集类型, 1 为交集, 0 为补集;

二、使用场景实例及其用法

使用场景设置:

我们手头有一个序列文件,并进行KEGG注释, 获得ko.tsv文件,但是我们只关注几个功能类型的基因,比如 K01190(beta-galactosidase), K00865( glycerate kinase), 将K01190、K00865列入文件:list.txt, 我们可以进行如下操作:

tabtk_subset -c 2  -s 1  ko.tsv list.txt
cat  list.txt | tabtk_subset -c 2  -s 1  ko.tsv  -
cat  ko.tsv | tabtk_subset -c 2  -s 1 - list.txt
tabtk_subset -c 2  -s 1  ko.tsv.gz  list.txt.gz

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

tabtk_subset -c 2  -s 1  <(pigz -dk -c -p 2  ko.tsv.gz)  <(pigz -dk -c -p 2  list.txt.gz)

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

Last Upate: 2017-09-27 2:06 PM