Archives

fasta/q序列抽取和排序: seqtk_reoder

一、使用场景

给定fasta/q文件和一个序列列表,按照给定的序列顺序抽取序列,如不存在则忽略,先前的 seqtk_subseq( 如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq ) 抽取出序列按照标识符在序列文件中出现的顺序输出,有时我们需要按照给定的标识符输出序列。

seqtk_reoder使用reoder作为名字输出就是指序列输出顺序按照指定顺序输出。

二、工具介绍

通过seqtk_info 查看测试序列信息:

$seqtk_info shear.fastq #sequence base min_len max_len avg_len 47525 23905075 503 503 503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_reorder Usage: seqtk_reorder <fasta/q> <list> version: 0.0.1

seqtk_reorder 接受两个参数, 序列文件和列表文件, 支持stream流模式;

实例

$cat shear.idx | seqtk_reorder shear.fasta – $cat shear.idx | seqtk_reorder shear.fastq.gz – $cat shear.idx | seqtk_reorder […]

fasta/q序列拆分三剑客: seqtk_split、seqtk_partition 和 seqtk_shred

一、使用场景

本部分内容主要解决以下几个问题:

对指定fasta/q序列文件拆分成每份N条序列问题,比如执行InterProScan。 对指定fasta/q序列文件拆分成指定N份问题,比如我们有40核心,程序不支持多线程(或者支持的不理想),那拆分40份多进程执行, 比如HMMER 对指定fasta/q序列文件序列拆分,每条序列保留最大l长度,并保证分割非序列存在N个碱基交叠。使用的场景多,比如 454 Newbler拼装软件,不支持很长度序列,如果我们对不同拼装版本的组装结果混合组装,将序列撕碎(保留比较大的overlap)再进行组装,就没有问题了。 二、工具介绍

以下三个程序分别解决以上三个问题。

通过seqtk_info 查看测试序列信息:

$seqtk_info shear.fastq #sequence base min_len max_len avg_len 47525 23905075 503 503 503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_split Usage: seqtk_split [options] <fasta/q> <size> <prefix> Options: -t fastx file suffix: default: [fasta]; -v print version number

seqtk_split 接受三个参数, 序列文件、 序列大小(条数)和输出文件前缀。 可以通过 -t 指定输入文件类型, 比如 […]

NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

本文介绍NGS数据(这里我们针对FASTQ序列数据)质量控制, 质量控制常见的内容有:

过滤低复杂度比例高的序列; 包含有’N’s的序列; 低质量区域,一般指5’端和3’端; phix序列; Expected errors 期望误差(这里按照这个字面解读);

其它几种比较常见, 这篇文章主要介绍期望误差值的方式。

期望误差可以参考:Expected errors predicted by Phred (Q) scores,简而言之就是将Phred值使用错误概率值(P_error)描述,对所有碱基的的P_error求和,该值我们可以称之为”E”值, 可以使用该值对序列进行过滤,文章: Error filtering, pair assembly and error correction for next-generation sequencing reads 描述的和详细, 并推荐使用E值阈值为 1.0。

P_error 、Q值 和 ASCII 表映射关系如下

Robert Edgar 在其网站上也描述了为什么使用平均Q值过滤不好, 参考:Average Q score

下面介绍几个小工具,如何使用E值模式进行序列质量控制。

一、fastq-filter

fastq-filter为 Usearch 实现的序列质量控制工具,命令行参数介绍:cmd_fastq_filter

支持的命令行选项:

-fastqout filename FASTQ 输出文件; […]

安全的合并序列文件:seqtk_cat

该小工具只解决一个小问题:安全的合并一堆 FASTA/Q 文件。 当一个文件末尾没有’\n’换行符的时候,合并序列文件会导致一条序列的’@’ 或者’>’ 会在另外一条序列的尾巴上,会造成很多程序出现故障。

seqtk_cat 命令行接口:

seqtk_cat Usage: seqtk_cat <fasta/q…> version: 0.0.1

实例:

seqtk_cat OS2_R1.fastq OS2_R2.fastq >/dev/null seqtk_cat *.fastq >/dev/null

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

Last Update: 2017-10-25 5:07 PM

双端序列链接: seqtk_join 和 fastq_join

本文介绍一个小工具: seqtk_join, 直接连接PE序列双端,序列中使用’N’和质量值通过’I’连接, 默认8个。 该程序主要使用场景为 metagenome PE序列合并,中间通过N并不影响 Kmer空间, 遇到N kmer会终止。

1. seqtk_join

命令行接口:

$ seqtk_join Usage: seqtk_join [options] <in.fa/q> Options: -n INT inserting padding number of ‘N’s, default:[8] -v print version number

实例:

seqtk_join -n 0 OS2_X_R1.fastq OS2_X_R2.fastq | head seqtk_join -n 8 OS2_X_R1.fastq OS2_X_R2.fastq | head 2. fastq_join

USEARCH 系统也有一个类似子命令:fastq_join

命令行模式:

usearch -fastq_join OS2_1.R1.fq […]

FASTQ/A 反向互补操作: seqtk_revcomp、seqtk seq -r 和 fastx_revcomp

这篇文章介绍的是最常见的序列操作: 反向互补,介绍三个Linux命令行实现, Web版本请移步: http://reverse-complement.com/

1. seqtk_revcomp

seqtk_revcomp 设计的目的当你需要查看一个引物的反向互补序列的时候,直接在命令行输入字符串,即可获得反向互补序列

比如:

$ seqtk_revcomp CCGTCAATTCmTTTRAGTTT AAACTYAAAkGAATTGACGG 2. seqtk seq -r

seqtk的序列处理工具箱,开箱即用。

seqtk seq -r seq.fna 3. fastx_revcomp

USEARCH 的 fastx_revcomp 子命令,支持fastq和fasta输入:

命令接口:

usearch -fastx_revcomp reads.fastq -label_suffix _RC -fastaout reads_rc.fa

可以通过-label_suffix修饰序列标识符。

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

对FASTA/Q序列添加注释信息: seqtk_comment

前面有一个小文章介绍如何去除序列文件的注释信息,今天遇到一个新的场景,需要将注释信息添加到序列文件里, 所以就有了 seqtk_comment。

如何快速去除fasta/q的注释信息:seqtk seq

命令接口

seqtk_comment Usage: seqtk_comment [options] <db> <in.fa/q> Options: -v print version number

需要两个参数,第一个参数为序列的注释文件, 制表符分隔; 第二个参数为序列文件;

注释信息是追加模式,可以连续添加,如果需要替换先前的注释信息,可以偶联seqtk seq -C去除所有的注释信息再添加就可以了。

实例:

cat seq.fna | seqtk_comment size.txt – | seqtk_comment annotation.txt – |head

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

Last […]

双端序列合并:fastq_mergepairs

当前16S扩增子测序实验中,比较常规的测序平台为HiSeq2500和Miseq,测序模式分别为 2×250 和 2×300 双端测序, 但目标片段长度一般大于单端序列(比如 V4区,V4-V5区), 双端序列通过交叠,然后合成一条序列可以完成整个扩增子序列测序,这里的合成就涉及双端序列合并问题(示例如下图)。

该类工具比较多,比如:Usearch fastq_mergepairs、FLASH、PANDASEQ、 COPE、BBmerge、PEAR、fastq-join、read-stitcher等。

后面会逐步使用 tag:mergepairs 汇总这类工具的使用。

fastq_mergepairs 为 USEARCH 程序包中的双端序列合并子命令,具体算法参考文章:Error filtering, pair assembly and error correction for next-generation sequencing reads 以及 cmd_fastq_mergepairs命令行介绍页面。

命令行接口:

usearch -fastq_mergepairs *R1*.fastq -relabel @ -fastq_maxdiffs 10 -fastq_pctid 80 \ -fastqout merged.fq

主要命令行参数:

-fastq_mergepairs R1端序列; -reverse 反向序列,如不提供,默认将R1替换成R2 -interleaved 支持interleaved reads 格式; -fastqout FASTQ […]

切除引物序列获取扩增目标序列:seqtk_pcr_strip

前面一篇文章: 对扩增子序列进行引物匹配: usearch 和 primersearch 遗留了一个问题,就是识别引物后如何处理? 保留引物序列还是切除引物序列的问题。

后面会讲到基于100%相似度构建OTU(或者feature)表的问题,引物区域带来的噪音就必须要考虑进去,所以这里介绍一个工具: seqtk_pcr_strip 切除引物序列获取扩增目标序列。

seqtk_pcr_strip 命令行接口:

$seqtk_pcr_strip Usage: seqtk_pcr_strip <amplicon> <fasta/q> Options: -l INT max overhung length, default [1] -v print version number

这里有一个参数: -l, 限制引物在扩增子序列位置, 设置为0, 意味着引物必须在序列的最左边和最右边,

实例

seqtk_pcr_strip amplicon-pcrout.txt amplicon-pcrout.fasta >amplicon-pcr_strip.fasta seqtk_pcr_strip amplicon-pcrout.txt amplicon-pcrout.fastq >amplicon-pcr_strip.fastq

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line […]

序列转换之RNA转DNA:seqtk_rna_to_dna 和 seqtk_base_mask

SILVA 为(当前版本 128)老牌的RNA序列数据库,提供的序列为RNA序列,顾名思义,序列里是 ‘U’ 而不是 ‘T’,但很多时候我们需要的是DNA序列,这样我们需要将序列里所有的’U’替换成’T’。

替换操作可以使用seqtk_rna_to_dna 专用程序或者 seqtk_base_mask 通用程序解决。

1. seqtk_rna_to_dna

seqtk_rna_to_dna 命令接口

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

实例:

seqtk_rna_to_dna SILVA_128_LSURef_tax_silva.fasta >SILVA_128_LSURef_tax_silva_dna.fasta 2. seqtk_base_mask

seqtk_base_mask 命令接口

$ seqtk_base_mask Usage: seqtk_base_mask <fasta/q> <base:from> <base:to> version: 0.0.1

需要显示指定从碱基’X’变换成碱基’Y’, 可以识别大写和小写。

实例:

seqtk_base_mask SILVA_128_LSURef_tax_silva.fasta U T >SILVA_128_LSURef_tax_silva_dna.fasta

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux […]