科研星球

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.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,使用第一个空格之前的字符作为序列的名字


0 (7).png

需要注意的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/


QQ客服
电子邮箱
淘宝官店
没有账号?