【原】seqkit：序列梳理神器-统计、格式转换、长度筛选、质量值转换、翻译、反向互补、抽样、去重、滑窗、拆分等30项全能

写在前面

通过我几天的学习，我发现，seqkit十分好用，将序列的各种操作都囊括进去，加入多线程，我个人认为这将是非常好的胶水，在处理无论是基因组还是其他组学。定是一个必学神器。注意一下教程在0.15版本以后可用。低版本有些参数不可用。

支持Windows/Mac/Linux的32和64位系统。用户根据自己的系统自取。

最新版发布页面：https://github.com/shenwei356/seqkit/releases

• 软件优点

1. 体积小巧、安装方便、无依赖关系

2. 跨平台：Windows/Linux/Mac通用

3. 运行效率高

安装

本教程测试环境Ubuntu 16/18/20.04 LTS

目前可以安装比较最新的版本Version: 至少 >= 0.15.0。seqkit在github上的维护力度也比较大，功能比较完善，放心使用。

conda install -c bioconda seqkit

# 补充bedops(gtf2bed)和csvtk工具(可选)
conda install -c bioconda bedops -y
conda install -c bioconda csvtk -y

可选在 https://github.com/shenwei356/seqkit/releases 发布页直接下载适合的版本。

seqkit 一共有37个可用的命令，详细内容如下：

amplicon 通过引物检索扩增子(或其周围的特定区域)
bam 检查和在线绘制BAM记录文件的直方图
common 通过id/名称/序列查找多个文件的公共序列
concat 连接多个文件中具有相同ID的序列
convert 转换FASTQ质量编码格式：支持格式包括：桑格，Solexa和Illumina
duplicate 重复序列N次
faidx 创建FASTA索引文件并提取子序列
fish 使用局部比对在较大的序列中寻找短序列
fq2fa 转换FASTQ到FASTA
fx2tab 将FASTA/Q转换为表格格式(包含长度/GC含量/GC偏好)
genautocomplete 生成shell自动完成脚本
grep 通过ID/name/sequence/sequence motif搜索序列，允许错配
head 打印第一条序列
help 打印帮助信息
locate 定位序列，或者motifs，允许错配
mutate 编辑序列(点突变、插入、删除)
pair 匹配双端序列文件
range 打印一个范围内的序列
rename 重命名重复序列ID
replace 使用正则表达式修改名称或者序列
restart 重置环状基因组的起始位置
rmdup 通过id/名称/序列删除重复的序列
sample 按数量或比例对序列进行抽样
sana 清理损坏的单行fastq文件
scat real time recursive concatenation and streaming of fastx files
seq 转换序列(反向，补充，提取ID…)
shuffle 随机序列
sliding 序列滑窗提取，支持环形基因组
sort 按id/名称/序列/长度排序序列
split 按id/seq区域/大小/部件将序列拆分为文件(主要用于FASTA)
split2 按序列数量/文件数将序列拆分为多个文件(FASTA, PE/SE FASTQ)
stats FASTA/Q文件的简单统计
subseq 通过region/gtf/bed得到子序列，包括侧翼序列
tab2fx 转换表格格式为FASTA/Q格式
translate 翻译DNA/RNA到蛋白质序列(支持歧义碱基)
version 打印版本信息并检查是否更新
watch 序列特征的监测和在线直方图

参数

Flags:
--alphabet-guess-seq-length int seqkit根据第一个FASTA记录猜测序列类型的序列前缀的长度(0表示整个序列)(默认10000)
-h, --help 显示帮助
--id-ncbi FASTA头是ncbi风格的，例如>gi|110645304|ref|NC_002516.2
--id-regexp string 用于解析ID的正则表达式(default "^(\\S+)\\s?")，匹配空格前的部分为序列名
--infile-list string 输入文件列表中的文件 (one file per line), if given, they are appended to files from cli arguments
-w, --line-width int 输出FASTA格式时的行宽 (0 for no wrap) (default 60)
-o, --out-file string 输出 ("-" for stdout, suffix .gz for gzipped out) (default "-") -代表标准输出，加.gz可输出压缩文件
--quiet 保持安静，不要显示额外的信息
-t, --seq-type string 序列类型 (dna|rna|protein|unlimit|auto) (auto, 按第一个序列自动检测) (default "auto")
-j, --threads int CPU数量 (默认单核为1，多核为2) (default 2)

实战

准备数据

# 创建练习目录并进入
mkdir -p seqkit
cd seqkit

# miRBase的RNA序列1.5M/785K/110K
wget -c ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.gz
wget -c ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz
wget -c ftp://mirbase.org/pub/mirbase/CURRENT/miRNA.diff.gz

# 下载人类基因组3G(压缩包840 MB)和基因注释gtf(44M) (可选)
# wget -c ftp://ftp./pub/release-84/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# wget -c ftp://ftp./pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz

# 下载拟南芥基因组(120M，压缩包35M) 比较小 推荐 http://plants./info/data/ftp/index.html
wget -c ftp://ftp.ensemblgenomes.org/pub/plants/release-49/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget -c ftp://ftp.ensemblgenomes.org/pub/plants/release-49/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.49.gtf.gz

# 下载fastq文件 用于测试
wget -c http://210.75.224.110/temp/meta/meta2010/seq/C1_1.fq.gz
wget -c http://210.75.224.110/temp/meta/meta2010/seq/C1_2.fq.gz

# 下载fa文件，宏基因组中prodigal预测结果
wget -c http://210.75.224.110/github/Note/Linux/data/gene.fa

根据gtf构建bed文件

下载的gtf文件似乎缺少一个transcript_id，这里补充一下。

zcat Arabidopsis_thaliana.TAIR10.49.gtf.gz|awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' | gtf2bed --do-not-sort | gzip -c > Arabidopsis_thaliana.TAIR10.49.bed.gz

提取1号染色体序列及注释作为示例

seqkit grep -p 1 Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -o chr1.fa.gz

# 注释信息按照染色体取子集：提取第一条染色体的基因组注释信息：'^1'
# 使用gtf文件提取
zcat Arabidopsis_thaliana.TAIR10.49.gtf.gz | grep -w '^1' | gzip -c > chr1.gtf.gz
# 提取第一条染色体的bed文件，用法相同
zcat Arabidopsis_thaliana.TAIR10.49.bed.gz | grep -w '^1' | gzip -c > chr1.bed.gz

stats FASTA/Q文件的简单统计

统计序列格式fasta(fa)/fastq(fq)、内容类型DNA/RNA/Protein，序列数量、总长度，最小、平均和最大长度

下面就四种格式序列构建和简单统计。

# FASTA DNA
echo -e ">seq\nacgtryswkmbdhvACGTRYSWKMBDHV" | seqkit stats
# RNA
echo -e ">seq\nACGUN\nACGUN" | seqkit stats
# Protein
echo -e ">seq\nabcdefghijklmnpqrstvwyz" | seqkit stats
# FASTQ DNA
echo -e "@read\nACTGCN\n+\n@IICCG" | seqkit stats

使用fq或者fa文件进行演示

# 一般模式
seqkit stats C1_1.fq.gz
#--输出结果tab分隔
seqkit stats C1_1.fq.gz -T
#--输出文件转化其他格式
seqkit stats C1_1.fq.gz -T| csvtk pretty -t
seqkit stats C1_1.fq.gz -T| csvtk csv2md -t
# 统计更多信息 -a
seqkit stats C1_1.fq.gz -a
# j多线程加速，尤其是对于具有多个序列文件会加速
# seqkit stats -j 2 *.fq.gz

seq 转换序列 （反向、互补/提取ID）

“-n”: 提取序列ID，包括“>”后面的全部内容

“-n -i”: 仅提取第一个空格前的ID

按长度过滤(常用)

扩增子分析时要筛选扩增长度相近的片段，过长或过短一般都要删除。宏基因组中比如组装的结果，经常要过滤<200/300bp的短片段，分箱时要筛选>1000/2000的长片段使用。本条命令非常多的应用场景。筛选后结果可用 > 写入文件

• -m 按照序列长度过滤，表示保留的最小长度，-M 此为保留的最大长度

#--提取序列长度大于60的并统计长度信息
zcat hairpin.fa.gz | seqkit seq -m 60 | seqkit stats
# 设置最小序列长度和最大序列长度，用于过滤序列，并统计
zcat hairpin.fa.gz | seqkit seq -m 100 -M 1000 | seqkit stats
# 保存>100且<1000长度的序列
seqkit seq -m 100 -M 1000 hairpin.fa.gz > hairpin100-1000.fa
seqkit stat hairpin100-1000.fa

提取ID

head gene.fa
## 名称全行
seqkit seq gene.fa -n | head
# 仅仅打印ID
seqkit seq gene.fa -n -i | head

# 使用正则表达式提取名字中的信息
zcat hairpin.fa.gz | head
# 提取ID中第二个字段作为ID
seqkit seq hairpin.fa.gz -n -i --id-regexp "^[^\s]+\s([^\s]+)\s" | head

单行/多行转换

• -s提取并展示序列

• -w 代表每行的碱基数量，0代表不换行

#仅仅提取序列 -s
seqkit seq gene.fa -s -w 0|head
#--将多行序列转化为标准4行FASTQ
seqkit seq C1_1.fq.gz -w 0|head

反向/互补

• -r 序列反向

• -p序列互补

# 序列反向互补,-r反向，-p互补
seqkit seq hairpin.fa.gz -r -p|head

删除gap/大小写转换

• -g 去除序列中的间隔,将中间的横杠去掉

• -u转化序列为大写字母展示

echo -e ">seq\nACGT-ACTGC-acc" | seqkit seq -g -u

RNA转为DNA

• —rna2dna 将RNA序列转化为DNA序列

echo -e ">seq\nUCAUAUGCUUGUCUCAAAGAUUA" | seqkit seq --rna2dna

subseq通过指定区域

• -r 通过区域来截取序列

如1:12提取前12个碱基，-12:-1提取序列结尾12个碱基；
for last 12 bases, 13:-1 for cutting first 12 bases. type “seqkit subseq -h” for more examples

#-提取序列前1：12个碱基
zcat C1_1.fq.gz | seqkit subseq -r 1:12 |head
#-提取序列最后1：12个碱基
zcat C1_1.fq.gz | seqkit subseq -r -12:-1 |head

#取第12至倒数第12个碱基，即前11和后11个碱基去掉
zcat C1_1.fq.gz | seqkit subseq -r 12:-12| head

基于gtf/bed信息挑选子序列。

• —gtf  根据gtf文件挑选基因，这部分功能用于根据基因注释快速提取基因序列，在宏基因组、转录组、重测序中常用。—chr 选择染色体，—feature cds选择序列类型

以拟南芥基因组的序列和注释数据演示：提取第一条染色体上的CDS基因信息，并统计基本信息

seqkit subseq --gtf Arabidopsis_thaliana.TAIR10.49.gtf.gz --chr 1 --feature cds Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz > chr1.gtf.cds.fa
seqkit stats chr1.gtf.cds.fa

• -u 可以提取目标基因上游的序列

• -f 目标区域不展示

#--挑选序列并多加上上游的3个碱基
seqkit subseq --gtf Arabidopsis_thaliana.TAIR10.49.gtf.gz Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -u 3 |head

# 仅提取上游序列，如提取启动子区2k：-f仅定位不输出位置序列，-u输出上游序列，此处示例3bp
seqkit subseq --gtf Arabidopsis_thaliana.TAIR10.49.gtf.gz Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -u 3 -f |head

sliding 滑窗提取序列，支持环状基因组

#-s 步长为3，-W 序列长度为6个碱基
echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6
# -g 贪婪模式，后面不足6个那也取
echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6 -g
# 环状DNA模式-C，首尾不算中断，环状
echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6 -C

步长为5，取30个碱基序列，然后统计GC含量

• fx2tab：统计fasta/fastq序列的信息为表格

• -n仅输出ID，不输出序列

• -g为GC含量

zcat hairpin.fa.gz | seqkit sliding -s 5 -W 30 | seqkit fx2tab -n -g |head

faidx 创建FASTA索引文件并提取子序列

# 解压
zcat hairpin.fa.gz > hairpin.fa

# 建索引*.fai文件
seqkit faidx hairpin.fa

# ID信息：hsa-let-7a-1 hsa-let-7a-2
seqkit faidx hairpin.fa hsa-let-7a-1 hsa-let-7a-2
# -f 标题行全部显示
seqkit faidx hairpin.fa hsa-let-7a-1 hsa-let-7a-2 -f
# 提取序列并选择区域显示
seqkit faidx hairpin.fa hsa-let-7a-1:1-10
seqkit faidx hairpin.fa hsa-let-7a-1:-10--1
seqkit faidx hairpin.fa hsa-let-7a-1:1
# 检查hsa开头的序列并统计
seqkit faidx hairpin.fa hsa -r | seqkit stats

watch 序列质量的监测和在线直方图

#-取对数展示直方图
seqkit watch -L -f ReadLen hairpin.fa

# 每五千个做一个图保存在pdf文件中
seqkit watch -p 500000 -O qhist.pdf -f MeanQual C1_1.fq.gz

• 从有错误记录的fastq文件中挽救可用的读取

sana：清理损坏fastq文件

这里我专门将C1_1.fq的第一个序列进行了错位，进行测试。这个操作往往在进行数据整合的时候可以有很大作用。

zcat C1_1.fq.gz|sed '2 s/^/A/' > C1_1_bad.fq
seqkit sana C1_1_bad.fq -o rescued.fq.gz

fq2fa 将fq转为fa格式

seqkit fq2fa C1_1.fq.gz -o C1_1.fa

fx2tab & tab2fx 序列转化表格格式

这一转化很有用，往往用于表格/矩阵处理的时候。

seqkit fx2tab hairpin.fa.gz | head -n 2

通过矩阵格式的序列文件统计序列长度和质量值

• -l 统计序列长度

• -g 统计平均GC含量

• -i 只打印名称(不打印序列)

• -H 打印标题行

# 打印序列长度、GC含量
seqkit fx2tab hairpin.fa.gz -l -g -n -i -H | head

• tab2fx 和表格格式转化为序列格式

# seqkit tab2fx 表格形式转化为序列形式
zcat hairpin.fa.gz | seqkit fx2tab | seqkit tab2fx | head

# 转化为表格，然后排序，然后转化回去
zcat hairpin.fa.gz \
| seqkit fx2tab -l \
| sort -t "`echo -e '\t'`" -n -k4,4 \
| seqkit tab2fx | head
# 等同于下面的命令
seqkit sort -l hairpin.fa.gz | head

# 通过这个转化可以将很多在表格中实现的数据处理方法用于序列
例如下面的提取1000个序列：seqkit fx2tab hairpin.fa.gz | head -n 1000 | seqkit tab2fx | head

translate 翻译DNA/RNA为蛋白质序列

#--转化为蛋白序列
seqkit translate gene.fa|head

# 去除'X' 和 '*'
seqkit translate hairpin.fa
seqkit translate hairpin.fa --trim | head

grep  序列匹配

# 生成一个ID列表
grep '>' C1_1.fa|cut -c2-|head -n 10 > id.txt

# 使用序列id列表进行搜索(不包含空格)
seqkit grep -f id.txt C1_1.fq.gz -o result.fq.gz
# 使用序列名称列表进行搜索(它们可能包含空格)
seqkit grep -n -f id.txt C1_1.fq.gz -o result.fa.gz
zcat result.fa.gz
# 提取hsa开头的序列
zcat hairpin.fa.gz | seqkit grep -r -p ^hsa |head

# -v参数v用于移除序列
zcat hairpin.fa.gz | seqkit grep -r -p ^hsa -p ^mmu -v | head

# 提取ID
zcat miRNA.diff.gz | grep ^# -v | grep NEW | cut -f 2 > list
head list
# 根据ID提取文件
zcat hairpin.fa.gz | seqkit grep -f list > new.fa
head new.fa

# 提取包含特定碱基组合的序列
cat hairpin.fa.gz | seqkit grep -s -i -p aggcg |head
# 统计
cat hairpin.fa.gz | seqkit grep -s -i -p aggcg | seqkit stats

# 去除 包含特定组合碱基的序列
zcat hairpin.fa.gz | seqkit grep -s -r -i -p ^aggcg |head

locate 输出匹配位置

＃　其他两种输出格式
zcat hairpin.fa.gz | seqkit locate -i -d -p AUGGACUN --bed
zcat hairpin.fa.gz | seqkit locate -i -d -p AUGGACUN --gtf

fish 使用局部比对在较大的序列中寻找短序列

echo -e '>seq\nACGACGACGA' \
| seqkit locate -p ACGA -G | csvtk -t pretty

echo -e '>seq\nACGACGACGA' \
| seqkit fish -F ACGA -a 2>&1 | csvtk -t pretty

amplicon 通过引物检索扩增子(或其周围的特定区域)

echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt

# 设置输出格式为bed，匹配的位点等信息
echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt --bed

#- 使用引物文件，这用于刘老师组的高通量分菌的序列拆分速度应该很客观
# cat seqs4amplicon.fa | seqkit amplicon -p primers.tsv --bed# 输出除去引物之外的部分，-r输出几个碱基
echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt -r 4:7
# 输出格式为bed
echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt -r 4:7 --bed

duplicate 对序列重复N次

# 重复序列1次，但是名字没有修改
cat hairpin.fa | seqkit head -n 1 \
| seqkit duplicate -n 2
# 对重复序列改名，使其独一无二
cat hairpin.fa | seqkit head -n 1 \
| seqkit duplicate -n 2 | seqkit rename

rmdup 通过id/名称/序列删除重复的序列

# 去除重复的序列
zcat hairpin.fa.gz | seqkit rmdup -s -o clean.fa.gz

# 保存重复序列得到文件 -D duplicated.detail.txt
zcat hairpin.fa.gz \
| seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

common ：通过id/名称/序列查找多个文件的公共序列

这里同时支持fa和fq文件

# 通过ID匹配，文件夹下全部的fa序列公共部分输出来
seqkit common *.fq.gz

# 通过-n实现全部名字匹配，-o输出结果
seqkit common *.fq.gz -n -o common.fq
# 通过-s序列匹配
seqkit common *.fq.gz -s -i -o common.fq

split 拆分序列为子文件

按名称ID、给定区域的子序列、文件大小或序列数量将序列拆分为文件

可用于将大文件拆分后，并行处理，加速分析。如从contig中预测基因。

#按照10000个序列为一个文件拆分，结果为hairpin.fa.gz.split/目录 ，文件名为hairpin.part_00x.fasta，-s
seqkit split hairpin.fa.gz -s 10000 -2

# 将序列拆分为四个部分(常用，等分然后并行)
seqkit split hairpin.fa.gz -p 5 -2

# 复杂一点的就是按照ID区分
seqkit split hairpin.fa.gz -i --id-regexp "^([\w]+)\-" -2

# 按照前三个序列碱基来区分
seqkit split hairpin.fa.gz -r 1:3 -2

split2 拆分文件 升级版本

Flags:
-l, --by-length string split sequences into chunks of >=N bases, supports K/M/G suffix
-p, --by-part int 按照拆分出来的数量，比如：拆分成两个子文件2。-s, --by-size int 按照序列数量拆分
-f, --force 强制覆盖文件
-h, --help 查看帮助文件
-O, --out-dir string 输出文件夹 (default value is $infile.split)
-1, --read1 string (gzipped) 双端序列第一个
-2, --read2 string (gzipped) 双端序列第二个

同时支持fa和fq文件。单端和双端序列拆分实例

-f强制覆盖结果，适合重复计算时使用

seqkit split2 hairpin.fa.gz -s 10000 -f

# 双端序列拆分(重点)，p指定拆分数量，-O指定输出目录，-f覆盖结果，默认为压缩
seqkit split2 -1 C1_1.fq.gz -2 C1_2.fq.gz -p 2 -O out -f

pair 拼接两个fastq文件

留下匹配的，去除不匹配的，这里我们使用扩增子的双端序列做一个演示：

注意：双端序列在两个文件中的顺序最好是一样的，否则会消耗大量内存去匹配。

seqkit pair -1 C1_1.fq.gz -2 C1_2.fq.gz -O result
# -u 输出未匹配上的文件
seqkit pair -1 C1_1.fq.gz -2 C1_2.fq.gz -O result -u -f

sample 按数量或比例对序列进行抽样。

按照百分比例和序列数量进行抽样

# 抽样百分之十
zcat C1_1.fq.gz | seqkit sample -p 0.1 -o sample.fq.gz
# 抽样1000条
zcat C1_1.fq.gz | seqkit sample -n 1000 -o sample.fq.gz

注意：1000条并不是很准确，可能是900多条，为什么呢？看这里了解问题。https://bioinf./seqkit/note/#effect-of-random-seed-on-results-of-seqkit-sample

这里为大家展示一下减少内存的序列抽样方法

# 抽样 seqkit sample
zcat hairpin.fa.gz \
| seqkit sample -p 0.1 \
| seqkit head -n 1000 -o sample.fa.gz

# 设置随机种子，方便重复结果: -s 11
zcat hairpin.fa.gz \
| seqkit sample -p 0.1 -s 11 |head

# 抽样后打乱序列 :seqkit shuffle
zcat hairpin.fa.gz \
| seqkit sample -p 0.1 \
| seqkit shuffle -o sample.fa.gz

range 打印序列 按照一个范围

# 打印一个范围内的序列
cat hairpin.fa | seqkit range -r 1:10
# 打印最后几行的序列
cat hairpin.fa | seqkit range -r -100:-1 | seqkit stats
# 打印中间范围的序列
cat hairpin.fa | seqkit range -r 101:150 | seqkit stats

repeat 使用正则表达式替换名称/序列。

# 修改序列名称：删除空格后内存
echo -e ">seq1 abc-123\nACGT-ACGT" \
| seqkit replace -p "\s.+"

# 修改序列名：替换
echo -e ">seq1 abc-123\nACGT-ACGT" \
| seqkit replace -p "\-" -r '='

# 修改序列：去除序列间隔
echo -e ">seq1 abc-123\nACGT-ACGT" \
| seqkit replace -p " |-" -s

# 修改序列：给每一个碱基加上空格
echo -e ">seq1 abc-123\nACGT-ACGT" \
| seqkit replace -p "(.)" -r '$1 ' -s

# 使用字符加数据重命名序列-用于扩增子代表序列改名非常优秀
echo -e ">abc\nACTG\n>123\nATTT" \
| seqkit replace -p .+ -r "ASV_{nr}"

echo -e ">abc\nACTG\n>123\nATTT" \
| seqkit replace -p .+ -r "SAV_{nr}" --nr-width 5

rename 重命名重复的ID

# 重命名：相同序列会在后面加上_2 来处理
echo -e ">a comment\nacgt\n>b comment of b\nACTG\n>a comment\naaaa" \
| seqkit rename

concat 连接多个文件中具有相同ID的序列

#这里演示组合前面两个碱基和最后两个碱基的用法
seqkit concat <(seqkit subseq -r 1:2 C1_1.fq.gz) <(seqkit subseq -r -2:-1 C1_2.fq.gz)|head

shuffle 随机打乱序列 默认全部读入内存

seqkit shuffle hairpin.fa.gz -2 > shuffled.fa

sort 按id/名称/序列/长度排序序列

• —quiet 屏幕不输出过程

• 
  

• -i 排序忽略大小写

• 
  

• -l 按照序列长度排序

# ID排序
echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" \
| seqkit sort --quiet

# 按照ID排序，忽略大小写
echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" \
| seqkit sort --quiet -i

# 按照序列长度排序，由小到大
echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAAnnn\n>seq3\nacgt" \
| seqkit sort --quiet -l

其他功能

mutate 编辑序列(点突变、插入、删除)

echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n"

# 修改第一个碱基
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -p 1:x

# 修改第五个位置的碱基，输出信息隐藏
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -p 5:x --quiet

# 可以同时修改多个碱基
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -p 1:x -p -1:x --quiet

# 删除碱基
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -d 1:1 --quiet

# 删除倒数三个碱基
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -d -3:-1 --quiet

# 插入碱基
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -i 0:xx --quiet

head 展示开头N行的序列

seqkit head -n 1 hairpin.fa.gz

locate 定位子序列或者保守序列位置

cat gene.fa | seqkit locate -p ACT | csvtk pretty -t
# 调整错配 最大错配为1
cat gene.fa \
| seqkit locate -p ACTG -m 1 \
| csvtk pretty -t

# 简并碱基
zcat hairpin.fa.gz \
| seqkit locate -i -d -p AUGGACUN \
| head -n 4

restart 重置循环基因组的起始位置

这个使用的比较少

echo -e ">seq\nacgtnACGTN"

echo -e ">seq\nacgtnACGTN" | seqkit restart -i 2
echo -e ">seq\nacgtnACGTN" | seqkit restart -i -2

convet 二代测序质量值的转化为 Sanger

# Illumina-1.8+ -> Sanger
# seqkit convert tests/Illimina1.8.fq.gz | seqkit head -n 1
# 40分以上的都认为40 Illumina-1.8+ -> Sanger
# seqkit convert tests/Illimina1.8.fq.gz -f | seqkit head -n 1

#Illumina-1.8+ -> Illumina-1.5+
# seqkit convert tests/Illimina1.8.fq.gz --to Illumina-1.5+ | seqkit head -n 1

# Illumina-1.8+ -> Illumina-1.5+ -> Sanger.
# seqkit convert tests/Illimina1.8.fq.gz --to Illumina-1.5+ | seqkit convert | seqkit head -n 1

# Solexa -> Sanger
# seqkit convert tests/Illimina1.8.fq.gz --from Solexa

# Illumina-1.5+ -> Sanger
# seqkit convert tests/Illimina1.5.fq | seqkit head -n 1

技术细节和使用

seqkit处理压缩文件

pigz 或者 gzip 在部分操作中不能加速，所以在v.0.8.1版本以便关注了，然而还是可以使用命令：

pigz -d -c seqs.fq.gz | seqkit xxx

因为seqkit使用了pgzip去写gzip。这比gzip和pigz更快。（10X of gzip, 4X of pigz），而且gzip压缩文件比较大。

从数据处理格式来讲

seqkit无缝支持fa和fq格式数据，并且可以自动识别。除了 faidx之外，全部命令都可以处理这两种格式的数据。

只有fa格式支持命令(subseq, split, sort和shuffle)利用FASTA索引(by flag -two-pass)下提高大文件的性能。

序列类型的检测DNA/RNA/Protein，会使用子序列进行，默认检测第一条子序列，通过—alphabet-guess-seq-length参数默认为10000，如果长度小于10000，则检查整条序列。

序列名字

所有的软件，包括seqkit，使用第一个空格之前的字符作为序列的名字：

需要注意的NCBI等一些序列的格式并不是如此，例如：

>gi|110645304|ref|NC_002516.2| Pseudomona

想要在seqkit中识别出来的序列ID为：NC_002516.2。

此时使用参数--id-regexp "\|([^\|]+)\| "，或者添加参数--id-ncbi，但如果是只要前面的gi数字作为ID的话，添加参数：--id-regexp "^gi\|([^\|]+)\|"。

注意：.seqkit.fai不同于samtools产生的.fai格式文件，seqkit使用整个序列开头而不是ID作为索引。

并行运算

单核CPU默认线程：—threads 1，多个CPU，线程默认2.

内存占用

seqkit许多的命令都不需要将整个序列读入到内存中。包括：stat, fq2fa, fx2tab, tab2fx, grep, locate, replace, seq, sliding, subseq。

注意：如果使用subseq —gtf | —bed时，如果GTF或者BED文件太大，内存使用量会暴增，可以通过指定染色体：—chr，或者—feature去限制特征。

有一些命令需要将文件读入内存，但是可以用过rmdup 和 common减少内存使用。

随机—抽样

抽样命令sample和shuffle使用了随机功能，为了保证重现性，可以使用-s设置随机种子。这可以保证在不同的环境中可以有相同的抽样结果。

编辑：文涛

责编：刘永鑫

参考文献

Wei Shen, Shuai Le, Yan Li & Fuquan Hu. (2016). SeqKit: a cross-platform and ultrafast toolkit for fasta/q file manipulation. PloS One 11, e0163962, doi: https:///10.1371/journal.pone.0163962

软件发布页：https://github.com/shenwei356/seqkit/releases

帮助文档：https://bioinf./seqkit/