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 […]

M	T	W	T	F	S	S
« Mar
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30	31