seqkit：FASTA/FASTQ序列处理

seqkit 是用 Go 写的 FASTA/FASTQ 处理工具，处理 50GB 序列文件比手写 awk 快近 10 倍。本文覆盖 seqkit 最高频的 8 个场景：统计、过滤、抽样、格式转换、序列提取。全部在 Debian 12 上实测。

1. 安装与验证#

1
conda install -c bioconda seqkit -y
2
seqkit version
3
# seqkit v2.13.0

2. stats——拿到FASTQ第一件事#

这是我最常用的 seqkit 命令，比 FastQC 快但给出最核心的信息：

1
seqkit stats *.fastq.gz

输出：

1
file          format  type  num_seqs    sum_len  min_len  avg_len  max_len
2
sample_R1.gz  FASTQ   DNA     30,000  4,500,000    150    150.0      150
3
sample_R2.gz  FASTQ   DNA     30,000  4,500,000    150    150.0      150

列	含义	关注点
num_seqs	总序列数	双端文件应该相等
sum_len	总碱基数	换算测序量
min_len/max_len	长度范围	发现异常短/长 reads
avg_len	平均长度	确认 read 类型

如果 min_len 远小于 avg_len——说明有严重降解或接头污染，需要先跑 fastp。

3. 序列过滤——按长度、质量、名称#

1
# 过滤≥100bp的序列
2
seqkit seq -m 100 input.fastq.gz -o long.fastq.gz
3

4
# 过滤≤1000bp的序列（保留短的）
5
seqkit seq -M 1000 input.fastq.gz -o short.fastq.gz
6

7
# 只保留ATCG组成的序列（去掉含N的）
8
seqkit grep -s -r -p "^[ATCG]+$" input.fasta -o clean.fasta
9

10
# 按名称过滤（正则匹配）
11
seqkit grep -n -r -p "chr1" genome.fa -o chr1.fa

4. 抽样——大文件测试时的救星#

分析 500GB 的全基因组数据前，先抽取 1% 跑一遍验证流程：

1
# 按比例抽样（1%）
2
seqkit sample -p 0.01 input.fastq.gz -o sample.fastq.gz
3

4
# 按数量抽样（取10000条）
5
seqkit sample -n 10000 input.fastq.gz -o sample_10k.fastq.gz
6

7
# 设置随机种子（保证可重复）
8
seqkit sample -n 10000 -s 42 input.fastq.gz -o sample_10k.fastq.gz

-s 设置随机种子的重要性： 审稿人可能问你”怎么抽样的”，有种子就能精确复现。

5. 格式转换——FASTA↔FASTQ#

1
# FASTQ → FASTA（丢掉质量信息）
2
seqkit fq2fa input.fastq.gz -o output.fasta
3

4
# FASTA → FASTQ（质量全设为最高，仅测试用）
5
seqkit fa2fq input.fasta -o output.fastq.gz

6. 序列提取——按ID列表批量取#

1
# 从id_list.txt提取对应的序列
2
seqkit grep -n -f id_list.txt input.fasta -o subset.fasta
3

4
# 排除这些ID（反向选择）
5
seqkit grep -n -v -f exclude.txt input.fasta -o filtered.fasta

id_list.txt 格式一行一个 ID，不需要 > 前缀。

7. 序列去重与排序#

1
# 按序列去重（保留第一条）
2
seqkit rmdup -s input.fasta -o dedup.fasta
3

4
# 按名称排序
5
seqkit sort -n input.fasta -o sorted.fasta
6

7
# 按序列长度排序（从长到短）
8
seqkit sort -l -r input.fasta -o by_length.fasta

8. 碱基组成统计#

1
seqkit fx2tab -n -g -c input.fasta
2
# 输出：序列ID  序列  GC含量
3

4
# 只看GC含量
5
seqkit fx2tab -n -g -c input.fasta | awk '{print $1, $NF}'

9. 踩坑#

坑1：gzip文件需要额外处理吗？ 不需要。seqkit 自动识别 .gz 后缀并解压，直接喂就行。

坑2：输出文件格式 seqkit 默认根据 -o 参数的后缀自动选格式（.gz=压缩，.fa=FASTA，.fq=FASTQ）。

坑3：大文件 grep 很慢 seqkit grep -n -f big_list.txt 在大文件上比 seqtk subseq 慢。如果列表很大（>10万条），用 seqtk 更快。

10. 小结#

场景	命令	比awk快多少
快速统计	`seqkit stats`	5-10x
按长度过滤	`seqkit seq -m 100`	3-5x
抽样	`seqkit sample -p 0.1`	10x+
格式转换	`seqkit fq2fa`	2-3x
序列提取	`seqkit grep -n -f list`	3-5x

本文于 2025-04-01 在 Debian 12 上实测，seqkit v2.13.0。