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

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

刘永鑫Adam

37571人浏览 · 2021-03-15 07:00:00

刘永鑫Adam · 2021-03-15 07:00:00 发布

写在前面

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

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

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

软件优点

体积小巧、安装方便、无依赖关系
跨平台：Windows/Linux/Mac通用
运行效率高

安装

本教程测试环境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.ensembl.org/pub/release-84/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# wget -c ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz

# 下载拟南芥基因组(120M，压缩包35M)  比较小 推荐 http://plants.ensembl.org/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.shenwei.me/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://doi.org/10.1371/journal.pone.0163962

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

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

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。