SeqKit — 超快速的 FASTA和FASTQ 文件操作工具包

生信菜鸟团 · 公众号 · 生物 · 2024-09-24 22:00

正文

工欲善其事必先利其器

SeqKit

SeqKit 是沈伟（重庆医科大学附属第二医院）老师开发的一个跨平台且极速的生物信息分析工具，专为处理 FASTA/Q 文件而设计。它提供了一系列命令行工具和 API 接口，可以进行序列文件格式转换、序列处理和分析、序列统计等操作。其特性包括不限于：

易于安装和使用：提供一键式安装和静态链接的执行文件，无需任何配置，即可直接运行。
高性能：针对大规模数据进行了优化，执行速度快，处理能力强。
多功能：涵盖38个子命令，满足多种生物信息学需求。
兼容性强：支持各种操作系统和压缩格式，并能轻松融入工作流中。
跨平台：支持 Windows、Linux 和 macOS 系统。
无依赖关系：无需安装额外的依赖软件即可使用。
支持多种格式：支持读取和写入包括 GZIP、XZ、ZSTD 和 BZIP2 在内的多种压缩格式文件。

编程语言：GO
Github：https://github.com/shenwei356/seqkit
文档：https://bioinf.shenwei.me/seqkit/

发表文献

题目：SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation
期刊：PLoS ONE
日期：2016年10月5日
作者&单位：沈伟 & 陆军军医大学
DOI：https://doi.org/10.1371/journal.pone.0163962

如何安装

软件提供多种安装方式，安装非常简便。详见：https://bioinf.shenwei.me/seqkit/download/

二进制安装

下载对应系统的压缩文件，解压即可使用

#https://github.com/shenwei356/seqkit/releases/download/v2.8.2/seqkit_darwin_arm64.tar.gz

wget -c http://app.shenwei.me/data/seqkit/seqkit_linux_amd64.tar.gz
tar -xf seqkit_linux_amd64.tar.gz 

./seqkit

conda 安装

conda activate RNAseq
conda install seqkit

功能简介

SeqKit 提供了一系列强大的子命令，用于高效地处理和分析基因序列。例如提取或过滤序列、获取简单统计信息、定位子序列、序列搜索和位置查找、移除重复序列、拆分文件等。它的使用非常灵活，通过子命令的组合使用，用户可以根据需要实现复杂的数据处理任务。

一个简单示例

假设现在我们有一个样本的fastq1和fastq2，如果数据没保存好，发生了部分损坏，很不幸有没有备份数据。此时我们只能尽量恢复能用的数据。对于这个需求，我们该如何解决？

如果不使用seqkit

参考抢救你破碎的测序数据。首先我们定位到出错的行，然后截取正确的行。

##取出第三行
awk 'NR%4==3 {print NR " " $0}' Cko1_R2.fq > ./test/row3_Cko1_R2.fq.txt
##判断第三行是不是正常？
head row3_Cko1_R2.fq.txt |cut -d " " -f 2

##取出第三行不正常的序列
awk '$2 !="+" {print $0}' row3_Cko1_R2.fq.txt > err_row3_Cko2_R1.fq.txt

##查看开始出错的行
head -n 1 err_row3_Cko1_R2.fq.txt 

##保留正确的行号，出错的行减去3(理论上是，但是还需要查看一下原始数据，或者保险起见，可以多截取一部分)，就是需要保留的
head -n 1 err_row3_Cko1_R2.fq.txt |awk '{print $1-7}' 
87508348

##截取正确的reads
head -n 87508348 Cko1_R1.fq|pigz -p 4 >./test/Cko1_1.fastq.gz
head -n 87508348 Cko1_R2.fq|pigz -p 4 >./test/Cko1_2.fastq.gz

使用seqkit

如果了解seqkit的子命令，那么这个需求实现起来就很简单。只需要两个子命令。

sana ：恢复格式错误的fastq文件

#Cko1_R2.fq.gz 是损坏的fastq文件

seqkit sana -j 4 ../Cko1_R1.fq -o rescued_Cko1_R1.fq  2>err1.log
seqkit sana -j 4 ../Cko1_R2.fq -o rescued_Cko1_R2.fq  2>err.log
seqkit sana -j 4 ../Cko1_R2.fq.gz -o rescued_Cko1_R2.fq.gz  2>err2.log

ls -lha rescued_Cko1_R* |cut -d " " -f 5-
9.2G 9月  24 17:09 rescued_Cko1_R1.fq
7.6G 9月  24 16:25 rescued_Cko1_R2.fq
1.7G 9月  24 17:17 rescued_Cko1_R2.fq.gz

$wc -l err* rescued_Cko1_R*
          1 err1.log
   14404520 err.log
  106861528 rescued_Cko1_R1.fq
   87508348 rescued_Cko1_R2.fq
   
$zcat rescued_Cko1_R2.fq.gz|wc -l
87508348

又因为fq1，fq2文件reads数如果不一致，后续分析会报错

pair: 匹配两对fastq文件中的成对reads

seqkit pair -j 4  -1 rescued_Cko1_R1.fq -2 rescued_Cko1_R2.fq -O pair_out

查看结果

$cd pair_out/

$ls -lha *.fq |cut -d " " -f 5-
7.6G 9月  24 17:51 rescued_Cko1_R1.fq
7.6G 9月  24 17:51 rescued_Cko1_R2.fq

$wc -l *.fq
   87508348 rescued_Cko1_R1.fq
   87508348 rescued_Cko1_R2.fq

stats 对比结果

seqkit stats Cko1_2.fastq.gz 
seqkit stats Cko1_1.fastq.gz

seqkit stats rescued_Cko1_R2.fq 
seqkit stats rescued_Cko1_R1.fq

更多用法示例见：
https://bioinf.shenwei.me/seqkit/faq/
注：对于大多数命令来说，使用四个线程已经足够快了，其中 FASTA/Q 读写是性能瓶颈，使用更多线程不会提高速度。

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：