Archives

Categories

根据引物序列提取barcode: seqtk_barcode_inspect

seqtk_barcode_inspect 这个工具的使用场景比较少见,也许在这样的场景下会使用到:

1. 扩增子测序,barcode和引物在一起: barcode+primer
2. 你从NCBI下载了一些扩增子序列原始数据,比如16S测序结果,barcode没有被去除
3. 你的数据分析流程需要知道barcode是什么,然后将所有序列合并在一起,重新拆分和去除barcode

如果你有上面的使用场景,也许这个小程序可以帮助你,设计思路很简单(如果不嫌麻烦,可以肉眼去看):

1. 使用引物前几个非简并碱基进行定位,目标的前几个碱基就是我们要找的barcode序列
2. 默认在前10条序列搜索,如果搜索到非兼并碱基,终止,返回barcode序列,如果没有匹配,终止程序;

程序接口:

Usage: seqtk_subseq [options] <in.fa> <primer> <sample>
Options:
  -t INT     top n sequence to search, default: 10
  -v         print version number

使用用例:

$seqtk_barcode_inspect -t 10 raw.split.PCF.1.fq  GTGCC   F
F       AGCACAT

“GTGCC” 为 16S 引物对V3的前几个碱基, AGCACAT就是对应的barcode

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

Last Update: 2017-09-21 4:51 PM

ENA和SRA数据预处理: SRA Toolkit 和 seqtk_sra

前面一篇文章讨论了如何使用 aspera 上传/下载 神器 Aspera:

数据高速上传下载的利器:Aspera及其在生命科学中的应用

告诉我们如何高效(其实最重要的还是本地带宽, wget照样可以高效!)下载数据,遗留点下载数据后如何处理的问题,下面可以继续探讨。

SRA和ENA下载的数据,一般非标准的fastq文件格式,有时需要预处理才能进行下一步数据分析,比如:

1. 从sra文件提取fastq/sff 文件
2. 格式化fastq文件

一、数据下载

按照先前套路下载一个测试文件:http://www.ebi.ac.uk/ena/data/view/ERS980542

ENA

ascp_ena_download  ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz
ascp_ena_download  ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz

SRA

ascp_sra ERR1145336

二、数据格式化

2.1 安装 NCBI SRA Toolkit

下载地址: https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/

ascp_ncbi_download  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/sratoolkit.2.8.2-centos_linux64.tar.gz
cd sratoolkit.2.8.2-centos_linux64/bin
export PATH=$PWD:$PATH

2.2 提取fastq文件
ENA下载的就是fastq文件,所以不需要额外处理,现在针对 sra 文件, 使用 fastq-dump 程序

fastq-dump的命令行参数:

$fastq-dump

Usage:
  fastq-dump [options] <path> [<path>...]
  fastq-dump [options] <accession>

Use option --help for more information

fastq-dump : 2.8.2

zhangl@localhost:/project/zhangl/sra-utils/SRA/ERR1145336$ fastq-dump  --help

Usage:
  fastq-dump [options] <path> [<path>...]
  fastq-dump [options] <accession>

INPUT
  -A|--accession <accession>       Replaces accession derived from <path> in
                                   filename(s) and deflines (only for single
                                   table dump)
  --table <table-name>             Table name within cSRA object, default is
                                   "SEQUENCE
PROCESSING

Read Splitting                     Sequence data may be used in raw form or
                                     split into individual reads
  --split-spot                     Split spots into individual reads

...

参数很多,不过大部分情况只需要如下操作:

fastq-dump  --split-3  --defline-seq '@$sn/$ri' --defline-qual '+'   ERR1145336.sra  -O  .

这样生成的fastq文件兼容性比较好了,比如(为展示方便格式进行了调整,fastq文件一个记录只有4行):

@CI25.31_9/1
GTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTTAAGCTAGATGTGAAAGCCCC
GCGCTCAACGTGGGAGGGTCATTTAGAACTGGCAGACTAGAGTCTTGGAGAGGGGAGTGGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGG
AACATCAATGGCGAAGGCAACTCCCTGGCCAAAGACTGACGCTCATGTGCG
+
CFFGGGGGGGGGGGGGGGGHHHHGGGGGGGGGGGGGGGGGHHGGGGGHHHHHHHHGGGGHHHGGGGGGGGGGGGCGGHHHHHGGGHHHHHHHHHHHHHHG
GGGGGGGHHGHHGGHGGGGHHHHHHHHGGGGGGGGGGGGGGGGGGGGGFGFGGGGFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
FFAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBDDDFFFFFFFF

2.3 格式化序列

如果下载的fastq数据格式比较乱: 比如这种格式(为展示方便格式进行了调整,fastq文件一个记录只有4行):

@CI25.31_9 926R.7_206 HWI-M02034:5:000000000-A3P1T:1:1101:12893:1645 1:N:0: orig_bc=GATCTG new_bc=GATCTG \ 
      bc_diffs=0 orig_bc=GCTAC new_bc=GCTAC bc_diffs=0
GTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTTAAGCTAGATGTGAAAGCCCCGC
GCTCAACGTGGGAGGGTCATTTAGAACTGGCAGACTAGAGTCTTGGAGAGGGGAGTGGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAACA
TCAATGGCGAAGGCAACTCCCTGGCCAAAGACTGACGCTCATGTGCG
+
CFFGGGGGGGGGGGGGGGGHHHHGGGGGGGGGGGGGGGGGHHGGGGGHHHHHHHHGGGGHHHGGGGGGGGGGGGCGGHHHHHGGGHHHHHHHHHHHHHHGGG
GGGGGHHGHHGGHGGGGHHHHHHHHGGGGGGGGGGGGGGGGGGGGGFGFGGGGFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAF
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBDDDFFFFFFFF

虽然格式没有问题,但是兼容性可能有问题,我们还是希望序列的名称是这样的:

1. @CI25.31_9/1
2. HWI-M02034:5:000000000-A3P1T:1:1101:12893:1645/1

使用 seqtk_sra 可以进行格式化:

seqtk_sra命令行接口:

seqtk_sra
Usage: seqtk_subseq [options] <in.fa>
Options:
  -f INT     reads name option: 0 reads name, 1: the first elements in comment fields, 2: the second ..., default: 0
  -o INT     reads orientation option: 1: forward, 2: reverse default: 1
  -v         print version number

通过下面两种模式操作就可以获得期望序列:

seqtk_sra  -f 0  -o 1 CI25.31.fastq.gz 
seqtk_sra  -f 1  -o 1 CI25.31.fastq.gz 

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

Last Update: 2017-09-20 3:02 PM

数据高速上传下载的利器:Aspera及其在生命科学中的应用

一、Aspera 介绍

Aspera是IBM公司的一款高速传输软件,在 2013年被IBM收购 前属于Aspera, Inc 公司。

IBM® Aspera® 快速文件传输使企业能够利用现有的 WAN 基础架构,以可预测、可靠且安全的传送方式,快速传输大型文件和数据集(不论是结构化还是非结构化数据),而与文件大小、传输距离和网络状况无关。该文件传输软件包含客户端和服务器软件包,可以使用状态不佳的网络,跨多个位置(甚至是远程位置),传输对时间要求严格的文件。

本材料来源:http://www-03.ibm.com/software/products/zh/high-speed-file-transfer

一句话:就是传输速度快,基本可以占满客户的带宽。

二、Aspera 在生命科学领域内的应用

NCBI(数据下载接口)/ENA/GigaScience 等都采用Aspera技术提供数据下载服务。

三、 下载安装

这部分内容只包含Linux(测试平台: CentOS 7 操作系统)平台的安装以及常用数据的操作,NCBI提供了Web接口,通过浏览器插件也可以很方便的下载数据。

下载和安装:

wget  http://download.asperasoft.com/download/sw/connect/3.7.2/aspera-connect-3.7.2.141527-linux-64.sh
sh aspera-connect-3.7.2.141527-linux-64.sh

如需其它平台请移步: http://downloads.asperasoft.com/en/downloads/8?list

文件会被安装在:/home/biostack/.aspera 目录 目录结构如下:

.
└── connect
    ├── bin
    ├── etc
    ├── iso-swid
    ├── lib
    ├── localization
    ├── notices.txt
    ├── plugins
    ├── product-info.mf
    ├── res
    └── var

9 directories, 2 files

默认路径是浏览器可以找到的目录,如果不是有浏览器插件,只使用命令行工具 ascp , 那可以随便将 .aspera 拷贝任何目录: 比如:

mv connect aspera-connect-3.7.2
mv aspera-connect-3.7.2 /biostack/tools/utils
export PATH=/biostack/tools/utils/aspera-connect-3.7.2/bin:$PATH

这样就可以愉快的使用命令行工具了:

$ascp -h
Usage: ascp [OPTION] SRC... DEST
          SRC to DEST, or multiple SRC to DEST dir
          SRC, DEST format: [[user@]host:]PATH
  -h,--help                       Display usage
  -A,--version                    Display version.
  -T                              Disable encryption
  -d                              Create target directory, implied for file/file-pair lists
  -p                              Preserve file timestamp
  -q                              Disable progress display
  -v                              Verbose mode
  -6                              Use IPv6
  -D...                           Debug level
  -l MAX-RATE                     Max transfer rate
  -m MIN-RATE                     Min transfer rate
                                  RATE: G/g(gig),M/m(meg),K/k(kilo)
  -u USER-STRING                  User-specific string

这样还不够,如果你有管理权权限的话:

sudo cp etc/asperaweb_id_dsa.openssh  /etc/asperaweb_id_dsa.openssh

这样会比较方便,没有也没有关系:记住asperaweb_id_dsa.openssh 路径:/biostack/tools/utils/aspera-connect-3.7.2/etc/asperaweb_id_dsa.openssh 就可以了。

下面介绍下如何下载NCBI和ENA的数据。

四、NCBI和ENA 下载/上传

为什么一定要掌握Aspera这款工具呢? 上传原始数据到NCBI SRA 数据库, 如果数据库超过10G,需要采用”preload” 模式, 即先上传到申请的目录, 然后在从文件夹选择文件完成上传。

NCBI

下载NCBI家的数据, 比如 BLAST 程序: 先前我们需要这样下载:

wget  ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-x64-linux.tar.gz

现在我们可以这样:

ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-x64-linux.tar.gz .

上传数据时需要预先申请 ”preload 目录”, 比如:

ascp -i /etc/asperaweb_id_dsa.openssh  --overwrite=diff  -QT -l1000m -k1 -d  upload_R1.fastq.gz  subasp@upload.ncbi.nlm.nih.gov:uploads/申请目录/PRJNAXXXXXXX/

只需要将: 申请目录替换成提交SRA申请到的目录就可以了, PRJNAXXXXXXX 一般设置成申请到的项目名称。

现在应该会下载和上传了,但是,有时候记住那么多东西很麻烦,每次执行其实就是替换我们的下载ftp地址,使用脚本程序可以简化接口,比如:

ascp_ncbi_download ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-x64-linux.tar.gz

那如果 SRA/SRP/SRX 呢,其实熟悉NCBI FTP 目录的应该很容易想到直接告诉脚本SRA/SRP/SRX号就可以直接找到具体目录进行下载 比如:

ascp_sra  ERR000002  ERR000001

然后自动翻译成:

ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/ERR/ERR000/ERR000002  ./
ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/ERR/ERR000/ERR000001  ./

进行下载.

这样够简单吧。

ENA

有NCBI为什么还要使用ENA呢,ENA和NCBI数据存储方式不同,illumina数据ENA提供了fastq.gz存储, NCBI统一成了sra, 下一篇文章介绍下, NCBI的 SRA Toolkit 以及如何 hack 结果。

ENA的数据下载可以参考:http://www.internationalgenome.org/faq/how-download-ena-files-using-aspera/

还以ERR000002 ERR000001为例子很像

直接命令行模式:

ascp -i /etc/asperaweb_id_dsa.openssh –overwrite=diff -QTr -l6000m -P33001 era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR208/002/SRR2088062/SRR2088062_1.fastq.gz .

同样,我们可以简单化:

ascp_ena_download ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR208/002/SRR2088062/SRR2088062_1.fastq.gz

五、几点说明

1. 选用ENA的fastq还是NCBI的sra在Twitter上讨论过很多,很多人倾向于fastq格式,sra后续处理还存在各种问题,可参考: https://gist.github.com/lh3/acf621374581be03f0de654c45cf5ca8

2. 经过广泛测试 Windows 平台的客户端下载NCBI SRA 数据不光从稳定还是易用性来说,要好不少,续传问题也很灵活。

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

Last Update: 2017-09-19 1:52 PM

获取fasta/q文件的碱基组成和长度信息: seqtk_info 和 seqtk_fqchk

本部分内容主要介绍如何快速获取fasta/q文件的碱基数和长度信息,使用 klib 实现,klib介绍参考前文章:

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

一、工具介绍

程序接口:

秉承支持streaming的原则,接口很简单:

$seqtk_info
Usage: seqtk_info <fasta/q>
version 0.0.1

$seqtk_fqchk
Usage: seqtk_fqchk [options] <in.fq>
Options:
  -q INT     offset value: 33 for Sanger quality, 64 for Illumina quality, default: 33
  -v         print version number

接受(gz压缩)文件名或者标准输入流(使用 – 占位)

示例介绍:

time  pigz -dk -c -p 2   test.fq.gz | seqtk_info  -
#sequence       base    min_len max_len avg_len
1000000 177535684       100     230     177.54

real    0m1.749s
user    0m2.795s
sys     0m0.346s

time  pigz -dk -c -p 2   test.fq.gz | seqtk_fqchk   -  >/dev/null

real    0m3.735s
user    0m6.075s
sys     0m0.404s

seqtk_fqchk 结果和 seqtk fqchk 结果类似,包含一些基本统计信息,以及A/T/C/G/N 出现的频率和每个位置的Q值的频数,设计这个程序的主要目的是绘制FASTQCd的碱基组成和碱基质量分布图

#sequence: 1000000; base: 177535684; min_len: 100; max_len: 230; avg_len: 177.54; Q20: 0.9967; Q30: 0.9908;

碱基组成图:

base

质量分布图

qual

二、其它工具

seqkit

seqkit stats(pigz)

time  seqkit stats  test.fq.gz
file        format  type   num_seqs      sum_len  min_len  avg_len  max_len
test.fq.gz  FASTQ   DNA   1,000,000  177,535,684      100    177.5      230

real    0m1.782s
user    0m0.882s
sys     0m0.293s

seqkit stats(gzip)

time  seqkit stats  test.fq.gz
file        format  type   num_seqs      sum_len  min_len  avg_len  max_len
test.fq.gz  FASTQ   DNA   1,000,000  177,535,684      100    177.5      230

real 0m3.647s
user 0m0.870s
sys  0m0.288s

time  pigz -dk -c -p 2   test.fq.gz | seqtk fqchk   -  >/dev/null

real    0m1.877s
user    0m3.580s
sys     0m0.397s

三、总结

该小工具很小,但是很实用。

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

Last update:2017-09-18 5:20 PM

快速将 fasta/q 格式化成制表格式:seqkit fx2tab和seqtk_tab

该小工具是帮助seqkit的fx2tab测试小程序,使用klib实现, klib介绍参考前文章:

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

本测试的起因是因为 tseemann(Prokka/Snippy作者)测试 seqkit fx2tab 性能时发现速度比较慢,我协助沈伟进行简单测试,毕竟实现seqtk_tab这个程序初衷也是为测试fx2tab。

seqkit 使用 GO 语言实现,性能很好,而且功能很多,值得拥有,gz文件使用pigz显著提高文件读操作,所以使用seqkit时应检查系统路径是否包含pigz

一、测试

检查pigz

$which  pigz
/biostack/tools/utils/pigz-2.3.4/pigz

测试清单:

$seqkit stats EAOA2.fq.gz
file         format  type    num_seqs        sum_len        min_len  avg_len  max_len
EAOA2.fq.gz  FASTQ   DNA     16,017,817      2,402,672,550  150      150      150

seqtk_tab

$time  seqtk_tab EAOA2.fq.gz  >/dev/null

real    0m29.091s
user    0m28.610s
sys     0m0.453s

$time  pigz -c -d  EAOA2.fq.gz | seqtk_tab -  >/dev/null

real    0m19.603s
user    0m36.531s
sys     0m4.786s

Shell实现

$time (gzip -c -d EAOA2.fq.gz | paste - - - - | cut -f 1,2,4 > /dev/null)

real    0m38.195s
user    1m7.577s
sys     0m5.141s

$time (pigz -c -d EAOA2.fq.gz | paste - - - - | cut -f 1,2,4 > /dev/null)

real    0m20.407s
user    1m2.547s
sys     0m8.182s

seqkit fx2tab

#with pigz
$time seqkit fx2tab EAOA2.fq.gz > /dev/null

real    0m21.934s
user    0m17.370s
sys     0m3.735s

#with pigz
$time seqkit fx2tab -j 8  EAOA2.fq.gz > /dev/null

real    0m23.544s
user    0m21.857s
sys     0m5.047s

#without pigz
$time seqkit fx2tab EAOA2.fq.gz > /dev/null

real 0m40.135s
user 0m19.697s
sys  0m4.412s

#without pigz
$time seqkit fx2tab -j 8  EAOA2.fq.gz > /dev/null

real 0m40.259s
user 0m22.970s
sys  0m5.013s

二、结论

首先这个场景下Shell的版本性能很高,这个在先前做 双端序列 interleave / deinterleave 测试场景已经得到过结论,不是主要讨论点, 序列interleave/deinterleave都基于klib实现了, seqtk mergepe 实现了interleave功能。

参考本测试以及沈伟做的各种不同机器下的测试:

update: 更新了最新测评统计图:

Github Issues请移步:https://github.com/shenwei356/seqkit/issues/25

结论显而易见: seqkit fx2tab 性能很好!

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

Last Upate: 2017/9/17 11:07

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

本部分内容主要介绍如何根据列表快速提取 fasta/q 序列子集/补集,主要解决下面两个问题场景:

 1. 交集 寻找出现在列表的序列子集
 2. 补集 寻找未出现在列表的序列子集

并集的问题我们也经常会遇到,比如我们手头上有很多序列集合,但是序列有重叠(序列名称一致、序列碱基/氨基酸一致),这个问题可以通过程序 seqtk_uniq 解决, 下次再讲.

一、 seqtk_subseq 介绍

下面讲如何使用seqtk_namesseqtk_subseq 这两个小工具解决这些问题, 程序实现主要是使用 klib (高效的轻量级C库,大部分是使用宏实现,学习如何使用klib, 请 移步 Wiki

下面是当前主要功能模块:

Common components

khash.h: generic hash table based on double hashing.
kbtree.h: generic search tree based on B-tree.
ksort.h: generic sort, including introsort, merge sort, heap sort, comb sort, Knuth shuffle and the k-small algorithm.
kseq.h: generic stream buffer and a FASTA/FASTQ format parser.
kvec.h: generic dynamic array.
kdq.h: generic double-ended queue (de-queue).
klist.h: generic single-linked list and memory pool.
kstring.{h,c}: basic string library.
kmath.{h,c}: numerical routines including MT19937-64 pseudorandom number generator, basic nonlinear programming and a few special math functions.
kson.{h,c}: simple JSON parser (no streaming)
kthread.c: simple multi-threading models.

Components for more specific use cases

kexpr.{h,c}: parsing and evaluating mathematical expressions
ksa.c: constructing suffix arrays for strings with multiple sentinels, based on a revised SAIS algorithm.
knetfile.{h,c}: random access to remote files on HTTP or FTP.
kurl.{h,c}: FILE-like interfaces to libcurl.
kopen.c: smart stream opening.
khmm.{h,c}: basic HMM library.
ksw.{h,c}: Striped Smith-Waterman algorithm.
knhx.{h,c}: Newick tree format parser.

有效使用这些函数,还用什么 Perl/Python 等动态语言完成计算/IO密集度比较高的数据处理任务呢?

因为我们经常遇到的是很大的数据集,为了测试下面小程序, 这里使用的测试数据都比较大, 并且从一个大文件中随机抽取 200000 条序列和 200 条序列以减小序列分布上的偏差。

材料清单:

序列集EAOA2.fq.gz合基本信息

seqtk_info EAOA2.fq.gz

#sequence       base            min_len    max_len    avg_len
16017817        2402672550      150        150        150.00

序列随机抽取:

seqtk sample   EAOA2.fq.gz  200 | gzip >200.fq.gz
seqtk sample   EAOA2.fq.gz  200000 | gzip >200000.fq.gz

seqtk_names: 快速获取序列的名字列表:

seqtk_names  200.fq.gz >200.txt
seqtk_names  200000.fq.gz >200000.txt

seqtk_subseq: 交集/补集

seqtk_subseq 命令行接口:

$ seqtk_subseq
Usage: seqtk_subseq [options] <in.fa> <name.list>
Options:
  -s INT     subset option: 1 intersection, 0 supplementary set, default: 1
  -v         print version number

使用-s 可选参数, 默认是交集模式:

$seqtk_subseq  EAOA2.fq.gz  200.txt
$seqtk_subseq -s 1 EAOA2.fq.gz  200.txt

前面两个命令是等价的, 补集模式:

seqtk_subseq -s 0 EAOA2.fq.gz  200.txt

另外: 命令行模式程序支持 streaming 模式可以提供很大便利, 比如我获取 文件 200.fq.gz 和 EAOA2.fq.gz文件的交集, 那么:

cat  200.fq.gz | seqtk_names - | seqtk_subseq  -s 1  EAOA2.fq.gz  -  | gzip  > intersection.fq.gz

就可以很愉快的解决问题;

二、其它相关工具介绍

前面两个小工具是本着一个程序解决一个问题的原则来实现的, 下面要介绍两个瑞士军刀类工具: seqtkseqkit

seqtk工具

seqtk 是 klib 作者 liheng 实现的处理序列文件的工具集合,下面是根据各种接口:

Usage:   seqtk <command> <arguments>
Version: 1.2-r94

Command: seq       common transformation of FASTA/Q
         comp      get the nucleotide composition of FASTA/Q
         sample    subsample sequences
         subseq    extract subsequences from FASTA/Q
         fqchk     fastq QC (base/quality summary)
         mergepe   interleave two PE FASTA/Q files
         trimfq    trim FASTQ using the Phred algorithm

         hety      regional heterozygosity
         gc        identify high- or low-GC regions
         mutfa     point mutate FASTA at specified positions
         mergefa   merge two FASTA/Q files
         famask    apply a X-coded FASTA to a source FASTA
         dropse    drop unpaired from interleaved PE FASTA/Q
         rename    rename sequence names
         randbase  choose a random base from hets
         cutN      cut sequence at long N
         listhet   extract the position of each het

seqtk的subseq可以解决这个问题:

seqtk subseq EAOA2.fq.gz  200.txt >/dev/null

但是只能求交集模式, 很多时候我们需要补集模式

seqkit 工具

下面是seqkit的命令接口, 相当丰富,可以解决很多很多问题

SeqKit -- a cross-platform and ultrafast toolkit for FASTA/Q file manipulation

Version: 0.7.1-dev

Author: Wei Shen <shenwei356@gmail.com>

Documents  : http://bioinf.shenwei.me/seqkit
Source code: https://github.com/shenwei356/seqkit
Please cite: https://doi.org/10.1371/journal.pone.0163962

Suggestion : Install pigz to gain better parsing performance for gzipped data

Usage:
  seqkit [command]

Available Commands:
  common          find common sequences of multiple files by id/name/sequence
  concat          concatenate sequences with same ID from multiple files
  convert         convert FASTQ quality encoding between Sanger, Solexa and Illumina
  dup             duplicate sequences N times
  faidx           create FASTA index file
  fq2fa           convert FASTQ to FASTA
  fx2tab          convert FASTA/Q to tabular format (with length/GC content/GC skew)
  genautocomplete generate shell autocompletion script
  grep            search sequences by pattern(s) of name or sequence motifs
  head            print first N FASTA/Q records
  help            Help about any command
  locate          locate subsequences/motifs
  range           print FASTA/Q records in a range (start:end)
  rename          rename duplicated IDs
  replace         replace name/sequence by regular expression
  restart         reset start position for circular genome
  rmdup           remove duplicated sequences by id/name/sequence
  sample          sample sequences by number or proportion
  seq             transform sequences (revserse, complement, extract ID...)
  shuffle         shuffle sequences
  sliding         sliding sequences, circular genome supported
  sort            sort sequences by id/name/sequence/length
  split           split sequences into files by id/seq region/size/parts
  stats           simple statistics of FASTA/Q files
  subseq          get subsequences by region/gtf/bed, including flanking sequences
  tab2fx          convert tabular format to FASTA/Q format
  version         print version information and check for update

Flags:
      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)
  -h, --help                            help for seqkit
      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...
      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")
  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)
  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")
      --quiet                           be quiet and do not show extra information
  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")
  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

Use "seqkit [command] --help" for more information about a command.

seqkit 执行提取序列

seqkit  common -j 8  EAOA2.fq.gz    200000.fq.gz  > /dev/null
seqkit  grep -f <(seqkit seq -n -i 200000.fq.gz) EAOA2.fq.gz

seqkit在使用多线程情况可以很大程序提高执行效率。

三、比较和总结

seqtk_subseq 使用klib 实现, 采用了很简单的Hash装载序列集合,然后遍历序列文件,所以对查询序列多少不是很敏感。

下面简单测试下几种情况下执行效率:

seqtk_subseq 版本

time seqtk_names 200000.fq.gz |  seqtk_subseq EAOA2.fq.gz  -  >/dev/null

real    0m30.265s
user    0m30.076s
sys     0m0.245s

time seqtk_names EAOA2.fq.gz  | seqtk_subseq  200000.fq.gz -  >/dev/null

real    0m37.530s
user    0m42.271s
sys     0m1.330s

time seqtk_names <(pigz -dk -c -p 2 200000.fq.gz) | \
seqtk_subseq  <(pigz -dk -c -p 2 EAOA2.fq.gz) -  >/dev/null

real    0m20.460s
user    0m10.840s
sys     0m1.440s    

seqtk subseq 版本

time seqtk_names 200000.fq.gz |  seqtk subseq EAOA2.fq.gz  -  >/dev/null

real    0m31.372s
user    0m31.201s
sys     0m0.242s

time seqtk_names EAOA2.fq.gz  | seqtk subseq  200000.fq.gz -  >/dev/null

real    0m38.925s
user    0m46.647s
sys     0m1.575s

time seqtk_names <(pigz -dk -c -p 2 200000.fq.gz) | \
seqtk subseq <(pigz -dk -c -p 2 EAOA2.fq.gz)  -  >/dev/null

real    0m20.082s
user    0m11.580s
sys     0m1.344s

seqkit common 版本

单线程模式(pigz)

time  seqkit  common   EAOA2.fq.gz    200000.fq.gz  > /dev/null
[INFO] read file: EAOA2.fq.gz
[INFO] read file: 200000.fq.gz
[INFO] find common seqs ...
[INFO] 200000 unique sequence IDs found in 2 files, which belong to 200000 records in the first file: EAOA2.fq.gz
[INFO] extract seqs from the first file: EAOA2.fq.gz

real    2m2.577s
user    2m42.880s
sys     0m16.566s

time  seqkit  common  200000.fq.gz EAOA2.fq.gz  > /dev/null
[INFO] read file: 200000.fq.gz
[INFO] read file: EAOA2.fq.gz
[INFO] find common seqs ...
[INFO] 200000 unique sequence IDs found in 2 files, which belong to 200000 records in the first file: 200000.fq.gz
[INFO] extract seqs from the first file: 200000.fq.gz

real    1m41.730s
user    2m43.950s
sys     0m12.143s

多线程模式(pigz)

time  seqkit  common -j 20  EAOA2.fq.gz    200000.fq.gz  > /dev/null
[INFO] read file: EAOA2.fq.gz
[INFO] read file: 200000.fq.gz
[INFO] find common seqs ...
[INFO] 200000 unique sequence IDs found in 2 files, which belong to 200000 records in the first file: EAOA2.fq.gz
[INFO] extract seqs from the first file: EAOA2.fq.gz

real    1m35.691s
user    3m17.335s
sys     0m16.924s

time  seqkit  common -j 20  200000.fq.gz  EAOA2.fq.gz  > /dev/null
[INFO] read file: 200000.fq.gz
[INFO] read file: EAOA2.fq.gz
[INFO] find common seqs ...
[INFO] 200000 unique sequence IDs found in 2 files, which belong to 200000 records in the first file: 200000.fq.gz
[INFO] extract seqs from the first file: 200000.fq.gz

real    1m1.913s
user    2m34.953s
sys     0m13.042s

seqkit grep 模式(pigz)

time  seqkit  grep -f <(seqkit seq -n -i 200000.fq.gz) EAOA2.fq.gz  > /dev/null

real    0m22.403s
user    0m16.333s
sys     0m3.403s

time  seqkit  grep -f <(seqkit seq -n -i EAOA2.fq.gz) 200000.fq.gz  > /dev/null

real    0m28.713s
user    0m42.685s
sys     0m6.438s

seqkit grep 模式(gzip)

time  seqkit  grep -f <(seqkit seq -n -i 200000.fq.gz) EAOA2.fq.gz  > /dev/null

real 0m41.321s
user 0m14.727s
sys  0m3.402s

time  seqkit  grep -f <(seqkit seq -n -i EAOA2.fq.gz) 200000.fq.gz  > /dev/null

real 0m49.428s
user 0m47.589s
sys  0m7.809s

综上: seqkit grep 模式 为做交集最佳实践 (系统需要安装pigz), 使用小文件进行查找性能更佳, 但是综上结果显示各个软件最佳状态差距没有太大,那就随意选吧,选最方便/最喜欢的那个。

update:

综上: 从seqkit在gzip和pigz模式执行效率看,性能还是比较依赖pigz, 当然C实现seqtk subseq 或者 seqtk_subseq 都在pigz解压(线程设置为2)时,seqkit性能还是没有超过seqtk subseq 或者 seqtk_subseq, 不过接口更干净。

不过今天要提到的重点是同时支持交集和补集功能的 seqtk_subseq, pgzip 还是可以提供很大的优势。

写本文档过程中和seqkit作者沈伟(Github)进行了有效沟通,再次表示感谢。

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

Last Upate: 2017/9/16 22:00

如何从命令行模式获取病毒序列

本文介绍如何使用命令行模式获取病毒序列,面对我们的只有Linux终端和命令行工具

NT库中的病毒序列

一、从 NCBI NT/NR 库中抽取序列

1.1 下载数据

从 NCBI NT 库中获得病毒序列无疑是最全面的,数据资源可在 BLAST 的 db库 下载。

ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nt.gz .

这里如何安装和使用 ascp 下载工具下一个小文章再写。

自从2016年9月份 NCBI放弃GI号 后, 下载的 NCBI FASTA 文件就只有Accession号,如果需要获得序列的 NCBI Taxonomy 信息, 我们需要做的就是获得 Accession号 到 Taxonomy Taxon号的映射文件, NCBI Taxonomy 数据库提供了次信息, 只需要下载 accession2taxid 文件:

accession2taxid 包含了以下映射文件:

nucl_est.accession2taxid.gz
TaxID mapping for live nucleotide sequence records of type EST.

nucl_gss.accession2taxid.gz
TaxID mapping for live nucleotide sequence records of type GSS.

nucl_wgs.accession2taxid.gz
TaxID mapping for live nucleotide sequence records of type WGS or TSA.

nucl_gb.accession2taxid.gz
TaxID mapping for live nucleotide sequence records that are not EST, GSS, WGS or
TSA.

prot.accession2taxid.gz
TaxID mapping for live protein sequence records.

为了获得每个Taxon Node 是否是属于病毒序列,需要回溯序列整个 “祖源” Lineages,以判断是否是病毒序列,所以需要下载 taxdump.tar.gz 文件

1.2 提取数据

现在可以开始提取NT库中病毒序列:

cat  nucl_wgs.accession2taxid.gz  nucl_gss.accession2taxid.gz  nucl_gb.accession2taxid.gz  nucl_est.accession2taxid.gz | \
./taxon-clade  ncbi.map  -   10239 | \
cut -f1 | seqtk_subseq  nt.gz -  1   \
>10239.accession.fasta

从此例可以看到管道可以减少多中间结果文件,善于使用命令行工具可以有效提高数据分析效率。

如果需要提取NR库中的病毒序列,可以如法炮制

二、从 NCBI Assembly/Refseq 数据库下载数据

2.1 NCBI Refseq Release 下载:

ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/refseq/release/viral .

直接解压就可以使用了, 序列标识符也是accession号。

2.2 NCBI Assembly 数据库下载:

ascp -i /etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m anonftp@ftp.ncbi.nlm.nih.gov:/genomes/Viruses .

需要合并各个基因组的信息才好用, 如果使用 Kraken 之类的工具, 需要对基因组序列进行加上标签,符合分类器的文件格式。

Last update: 2017-09-16 1:42 AM

使用 Vagrant 在 CentOS 7 平台上使用命令行模式访问VirtualBox虚拟机镜像

前面介绍过在 Windows 平台安装虚拟机解决一些生物信息问题,现在要讲的是在 Linux 平台安装和使用虚拟机, 当然虚拟机平台还是 选用 VirtualBox,测试场景还是 CentOS 7 操作系统。

在安装VirtualBox和vagrant之前,需要准备好sudo用户或者root用户权限

这里不讨论为什么不使用 Docker 而使用 VirtualBox。

1 安装 VirtualBox

直接rpm安装,也可以下载下来安装:

sudo rpm -ivh  http://download.virtualbox.org/virtualbox/5.1.26/VirtualBox-5.1-5.1.26_117224_el7-1.x86_64.rpm

这篇文档 VirtualBox 不是主角,所以写的不多。

2 安装 vagrant

直接rpm安装,也可以下载下来安装:

sudo rpm -ivh  https://releases.hashicorp.com/vagrant/2.0.0/vagrant_2.0.0_x86_64.rpm

vagrant的镜像叫 box, 这里有一个 vagrant box 仓库 vagrantbox 可以使用, 包含了各种Linux的发行版本, 提及有大有小, 看阉割的程度了。

vagrant 有丰度的命令行接口:

Usage: vagrant [options] <command> [<args>]

    -v, --version                    Print the version and exit.
    -h, --help                       Print this help.

Common commands:
     box             manages boxes: installation, removal, etc.
     connect         connect to a remotely shared Vagrant environment
     destroy         stops and deletes all traces of the vagrant machine
     global-status   outputs status Vagrant environments for this user
     halt            stops the vagrant machine
     help            shows the help for a subcommand
     init            initializes a new Vagrant environment by creating a Vagrantfile
     login           log in to HashiCorp's Vagrant Cloud
     package         packages a running vagrant environment into a box
     plugin          manages plugins: install, uninstall, update, etc.
     port            displays information about guest port mappings
     powershell      connects to machine via powershell remoting
     provision       provisions the vagrant machine
     push            deploys code in this environment to a configured destination
     rdp             connects to machine via RDP
     reload          restarts vagrant machine, loads new Vagrantfile configuration
     resume          resume a suspended vagrant machine
     share           share your Vagrant environment with anyone in the world
     snapshot        manages snapshots: saving, restoring, etc.
     ssh             connects to machine via SSH
     ssh-config      outputs OpenSSH valid configuration to connect to the machine
     status          outputs status of the vagrant machine
     suspend         suspends the machine
     up              starts and provisions the vagrant environment
     validate        validates the Vagrantfile
     version         prints current and latest Vagrant version

For help on any individual command run `vagrant COMMAND -h`

Additional subcommands are available, but are either more advanced
or not commonly used. To see all subcommands, run the command
`vagrant list-commands`.

下面几个命令,第一次使用会用到:

vagrant init 
vagrant up
vagrant ssh
vagrant global-status

3 安装 测试

现在 CentOS 5 以及不再维护 (维护截至日期 2017-03-31), 所以这里测试用例不再使用CentoS 5, 从vagrantbox可以选一个 CentOS 7, 最好包含 VirtualBox Guest Additions, 这样文件共享也比较方便, 不过这里演示 ssh 登录

mkdir ~/centos-7.0 && cd ~/centos-7.0
wget https://github.com/tommy-muehle/puppet-vagrant-boxes/releases/download/1.1.0/centos-7.0-x86_64.box
vagrant box add centos-7.0  centos-7.0-x86_64.box

通过 vagrant box list, 可以查看到已经加入box列表;

[biostack@localhost centos-7.0]$ vagrant  box list
centos-7.0 (virtualbox, 0)

下面可以开始挂载:

vagrant init  centos-7.0
vagrant up

初始化完后会生成一个 Vagrantfile 配置文件, 经过一系列的初始化工作,通过 global-status 可以查看已挂在的box

[biostack@localhost centos-7.0]$ vagrant global-status
id       name    provider   state   directory
------------------------------------------------------------------------
9fe8b66  default virtualbox running /home/biostack/centos-7.0

The above shows information about all known Vagrant environments
on this machine. This data is cached and may not be completely
up-to-date. To interact with any of the machines, you can go to
that directory and run Vagrant, or you can use the ID directly
with Vagrant commands from any directory. For example:
"vagrant destroy 1a2b3c4d"

登录box

[biostack@localhost centos-7.0]$ vagrant ssh
Last login: Thu Jul 16 08:48:31 2015 from 10.0.2.2
Welcome to your Vagrant-built virtual machine.
[vagrant@localhost ~]$

用户名和密码都是vagrant, 并且有sudo用户权限,然后通过一些工具可以和宿主进行交互文件, 比如 rsync, 如果没安装:

sudo  yum  install rsync

轻松搞定:

从宿主机器拷贝文件:

[vagrant@localhost ~]$ rsync  -avP  biostack@192.168.0.101:/biostack/distribution/binaries/zstd-1.3.1-x86_64.tar.gz  .
The authenticity of host '192.168.0.101 (192.168.0.101)' can't be established.
ECDSA key fingerprint is 03:23:c0:ea:5e:54:81:90:3b:f7:aa:a9:5e:08:c1:87.
Are you sure you want to continue connecting (yes/no)? yes
Warning: Permanently added '192.168.0.101' (ECDSA) to the list of known hosts.
biostack@192.168.0.101's password:
receiving incremental file list
zstd-1.3.1-x86_64.tar.gz
     2936323 100%   56.01MB/s    0:00:00 (xfer#1, to-check=0/1)

sent 30 bytes  received 2936797 bytes  345509.06 bytes/sec
total size is 2936323  speedup is 1.00

然后你可以做你想做的事情了。

虚拟机文件存放在:

/home/biostack/VirtualBox VMs/centos-70_default_1505385849643_38205

虚拟机可以通过 VBoxManage list vms 查看到名字和UUID

不用了清除掉虚拟机

vagrant destroy  9fe8b66

记得 下次想再使用从 vagrant up 开始,但是先前的/home/biostack/VirtualBox VMs/centos-70_default_1505385849643_38205 执行完 vagrant destroy 就没有了, 下次加载是全新的虚拟镜像。

4 结语

了解到 Vargant和 VirtualBox是通过Liheng的 Github 仓库 centos5-vm, 使用centos5-vm初衷是编译Linux平台兼容的生物信息程序,上面介绍的也比较详细了, 这里只留下 sourceforge devtools 的地址。

Vargant 还有很多丰度的功能,后面慢慢介绍。

Last update: 2017-09-14 7:04 PM

CentOS 7 平台安装R 环境

本文档在CentOS 7 环境测试

1 下载R

wget https://cran.r-project.org/src/base/R-3/R-3.4.1.tar.gz

2 编译R

tar xzvf  R-3.4.1.tar.gz
cd R-3.4.1
./configure  --prefix=$PWD   --enable-R-shlib  --with-tcltk -with-tcl-config=/usr/lib64/tcl8.5/tclConfig.sh --with-tk-config=/usr/lib64/tkConfig.sh
make -j
make install

几点注意:

安装 tcl/tk 开发环境 sudo yum install tcl tcl-devel tk tk-devel
添加编译选项   --enable-R-shlib 可以支持其它语言和R进行交互,比如rpy2

3 添加环境变量及其镜像

设置环境变量:

export PATH=/biostack/tools/devtools/R-3.4.1/bin:$PATH

设置镜像:

永久设置CRAN镜像, 镜像汇总 http://cran.r-project.org/mirrors.html , 方法参考 Stackoverflow: Set default CRAN mirror permanent in R

vim ~/.Rprofile

options(repos=structure(c(CRAN="http://mirrors.ustc.edu.cn/CRAN/")))

安装及设置Bioconductor镜像:

设置中国镜像:http://mirrors.ustc.edu.cn/bioc/ , 参考:http://www.tengfei.name/cn/2013/03/ustc-bioc/

source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite()

4 安装一些常用的包

biocLite('goseq')
biocLite('qvalue')
biocLite('topGO')
biocLite('GOstats')
biocLite('DEseq2')
biocLite('DEseq')
biocLite('edgeR')
install.packagese('data.table')
install.packages('hash')
install.packages('stringr')
install.packages('rlist')
install.packages('devtools')
install.packages('docopt')
install.packages('ggplot2')
install.packages('reshape')
install.packages('dplyr')
install.packages('vegan')
install.packages('gtools')
install.packages('pheatmap')
install.packages('optparse')
install.packages('data.table')
install.packages('data.table')
install.packages('plot2groups')

Last update: 2017-09-13 5:55 PM

CentOS 7 系统三步手动挂载硬盘 – GUI篇

第一步:选中Disk管理

方式: Applications -> Utilities -> Disk

Disk

第二步:格式化磁盘

通过配置按钮,选择格式化,进行格式化硬盘,默认选择 ext4 格式

2-1

2-2

第三步:挂在磁盘

选择挂在磁盘目录, 比如 Biostack 系列服务器, 选择挂在在根目录

创建根目录挂载点:

cd /
mkdir biostack

在挂在选项挂载点设置: /biostack, 并关闭自动挂载选项

此处输入图片的描述

来自Biostack团队: 2017-08-03 版本