SeqKit -- a cross-platform and ultrafast toolkit for FASTA/Q file manipulation
01 fasta文件与文件处理工具
1.1 fastq格式(常简写为.fq)
Fastq格式是NGS测序出来的原始数据格式。Fastq格式的序列一般有4行。第一行是元数据, 由‘@’开始,后面跟着序列的描述信息, 通常被认为是标志行。第二行是序列,由{ACGTN}组成, 其中N代表不确定,表示该位置可能是ACG或T。第三行由‘+’开始, 后面也可以跟着序列的描述信息, 通常和@后面的内容一样。第四行是第二行序列的质量评价,因此字符数和第二行一样,表示对应碱基的错误概率的相关信息,越大代表错误概率越低,通常用ASCII表示。该格式已成为生物信息学领域的一项标准。尤其在NGS和三代测序领域使用最频繁。
核酸数据
@A00402:166:HFTWGDSXY:4:1101:4526:1000 1:N:0:CCGAAG
GGCAGAAGCAATCTTAAGATCATCGGATTCGGGGGCGGCGGCGCACGGGATTTTTCGCGAATCCTCCTTCACCATCAAGCGATCAAAATCGCTGATCTCGTCTTCGTCATCCGTGATATCTTCGTGCGCAGCGCGGGGTATGATGAGCGAC
+
FFFFFFFFFF:FFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFF
1.2 fasta格式(常简写为.fa)
fasta格式,又称Pearson格式,主要发明人是威廉·皮尔森(William Raymond Pearson)和戴维德.李普曼(David J. Lipman),威廉·雷蒙德·皮尔森是美国弗吉尼亚大学的生物化学与分子遗传学教授,戴维德.李普曼在1989年至2017年期间担任NCBI主任,他也是BLAST算法的发明人之一。fasta格式是一种基于文本用于表示核酸序列或多肽序列的格式。其中核酸或氨基酸均以单个字母来表示,且允许在序列前添加序列名及注释。该格式和fastq格式已成为生物信息学领域的一项标准。
Fasta格式通常有两行。第一行由大于号“>”(较常用)开始, 后面跟着序列的描述信息, 通常被认为是标志行,可能包括仪器名, 测序平台, 短读的标识等信息。第二行是序列,序列信息通常由ACGT四个碱基组成,有时会有N代表不确定该位置是ACG还是T。仅从数据内容的格式上而言,相比于Fastq格式数据, 它缺少了Fastq格式中的第三行和第四行内容。
核酸序列
>TRINITY_DN0_c0_g1_i1.p1 TRINITY_DN0_c0_g1~~TRINITY_DN0_c0_g1_i1.p1 ORF type:3prime_partial (-),score=52.42 len:1017 TRINITY_DN0_c0_g1_i1:3-1019(-)
ATGCGTGTGTTTTGTGTTTGTTACATTTTCAGGTGTTTCGAATCATTAATTGCGGATCTG
TTGCTCTTGTTTTCACCAACAAAAACCCTAGCTATCGGTCGTTGTTCGTCGCAGTTTTTG
TTTCACATATTGCATGTGCTTAGCTTAGCTATGGATCGGTACCAGAAGGTGGAGAAGCCG
AAGGTTGAGACGCCAATCGCTGAGAATGAAATTCGGATTACAAGCCAGGGCAGGATGCGA
AGCTACATCACTTACGCCATGACTCTGCTTCAGGAAAAGGGCTCAGATGAAGTAGTATTC
AAGGCAATGGGCAGAGCCATCAATAAAACTGTTACAATTGTGGAGTTAATCAAGAGAAGA
ATTGTTGGTCTTCACCAAATTACAGCAATTCAATCCACTGATATAACTGACACATGGGAA
蛋白序列
>TRINITY_DN0_c0_g1_i1.p1 TRINITY_DN0_c0_g1~~TRINITY_DN0_c0_g1_i1.p1 ORF type:3prime_partial (-),score=52.42 len:339 TRINITY_DN0_c0_g1_i1:3-1019(-)
MRVFCVCYIFRCFESLIADLLLLFSPTKTLAIGRCSSQFLFHILHVLSLAMDRYQKVEKP
KVETPIAENEIRITSQGRMRSYITYAMTLLQEKGSDEVVFKAMGRAINKTVTIVELIKRR
IVGLHQITAIQSTDITDTWEPLEEGLQILETTRKVSMVTITLSKKDLDMNNVGYQPPIPA
DQVKVSTELEYDGEGSPIGRGRGRGGRGRGRPRGGPAGNGYAPAEFDDGGYDRSRGYPRG
RGRGRGRNFRGRGRGGYYQSDAQNDAGGRGRGGYYQSDAQNDAGGRGRGGYYQSDAQNDA
GGRGRGGYYQSDAQNDAGGRGRGGYYQSDAQNDAPRYNQ
1.3 序列编辑工具
如何对.fa文件进行编辑,大家有用记事本(win),或者其他win版本txt编辑工具诸如typora等等,同时伴随着markdown语法的兴起,笔者观察到越来越多的研究生在使用win电脑或者mac电脑时候开始用txt做记录,用word等的频率开始下降。
SeqKit - a cross-platform and ultrafast toolkit for FASTA/Q file manipulation横空出世,在github上已获得1k以上star,极大地帮助了linux用户批量编辑处理.fa/.fq文件。
- Ultrafast (see technical-details and benchmark)
- Seamlessly parsing both FASTA and FASTQ formats
- Supporting (
gzip
/xz
/zstd
/bzip2
compressed) STDIN/STDOUT and input/output file, easily integrated in pipe - Reproducible results (configurable rand seed in
sample
andshuffle
) - Supporting custom sequence ID via regular expression
- Supporting Bash/Zsh autocompletion
- Practical functions supported by 38 subcommands
02 参考
https://github.com/shenwei356/seqkit #官网
Author: Wei Shen <shenwei356@gmail.com>
Documents : http://bioinf.shenwei.me/seqkit
Source code: https://github.com/shenwei356/seqkit
Please cite: https://doi.org/10.1371/journal.pone.0163962
03 安装
wget -c https://github.com/shenwei356/seqkit/releases/download/v2.5.1/seqkit_linux_amd64.tar.gz
#通过这个链接安装最新版,then使用即可
#或者conda安装亦可
Light weight and out-of-the-box, no dependencies, no compilation, no configuration
conda install -c bioconda seqkit
04 使用
##38个子命令支持的实用功能
使用说明:
- `seqkit [命令]`
可用命令:
- amplicon 通过引物提取扩增子(或其周围特定区域)
- bam 监控BAM记录特征的在线直方图
- common 通过id/name/sequence找出多个文件的共有序列
- concat 从多个文件中连接具有相同ID的序列
- convert 在Sanger、Solexa和Illumina之间转换FASTQ质量编码
- duplicate 将序列复制N次
- fa2fq 通过FASTA文件检索相应的FASTQ记录
- faidx 创建FASTA索引文件并提取子序列
- fish 使用局部对齐在较大序列中查找短序列
- fq2fa 将FASTQ转换为FASTA
- fx2tab 将FASTA/Q转换为表格格式(及长度、GC含量、平均质量等)
- genautocomplete 生成shell自动补全脚本(bash|zsh|fish|powershell)
- grep 通过ID/名称/序列/序列基序搜索序列,允许不匹配
- head 打印前N个FASTA/Q记录
- head-genome 打印名称中有共同前缀的第一个基因组的序列
- locate 定位子序列/基序,允许不匹配
- merge-slides 合并由seqkit sliding生成的滑动窗口
- mutate 编辑序列(点突变、插入、删除)
- pair 从两个fastq文件中匹配配对末端读取
- range 在范围内打印FASTA/Q记录(start:end)
- rename 重命名重复的ID
- replace 通过正则表达式替换名称/序列
- restart 为圆形基因组重置起始位置
- rmdup 通过ID/名称/序列移除重复序列
- sample 按数量或比例抽取序列
- sana 清理损坏的单行FASTQ文件
- scat 实时递归连接和流式传输fastx文件
- seq 转换序列(提取ID、按长度过滤、移除间隙、反向互补等)
- shuffle 洗牌序列
- sliding 在滑动窗口中提取子序列
- sort 按id/name/sequence/length排序序列
- split 按id/序列区域/大小/部分将序列分割成文件(主要针对FASTA)
- split2 按大小/部分将序列分割成文件(FASTA,PE/SE FASTQ)
- stats FASTA/Q文件的简单统计
- subseq 通过区域/gtf/bed获取子序列,包括侧翼序列
- sum 为FASTA/Q文件中的所有序列计算消息摘要
- tab2fx 将表格格式转换为FASTA/Q格式
- translate 将DNA/RNA翻译成蛋白质序列(支持模糊基础)
- version 打印版本信息并检查更新
- watch 监控序列特征的在线直方图
标志:
--alphabet-guess-seq-length int 根据第一个FASTA记录的序列前缀长度猜测序列类型(0表示整个序列)(默认10000)
--compress-level int gzip、zstd、xz和bzip2的压缩级别。键入 "seqkit -h" 查看每种格式的范围和默认值(默认-1)
-h, --help seqkit的帮助
--id-ncbi FASTA头部是NCBI风格的,例如:>gi|110645304|ref|NC_002516.2|Pseud...
--id-regexp string 解析ID的正则表达式(默认 "^(\\S+)\\s?")
--infile-list string 输入文件列表的文件(每行一个文件),如果给定,它们将被附加到cli参数的文件中
-w, --line-width int 输出FASTA格式时的行宽(0表示不换行)(默认60)
-o, --out-file string 输出文件("-"表示标准输出,后缀.gz表示压缩输出)(默认"-")
--quiet` 保持安静,不显示额外信息
-t, --seq-type string 序列类型(dna|rna|protein|unlimit|auto)(对于auto,它会自动通过第一个序列检测)(默认"auto")
-j, --threads int CPU数量。也可以通过环境变量SEQKIT_THREADS设置)(默认4)
使用 "seqkit [命令] --help" 获取有关命令的更多信息。
SeqKit 使用 pgzip([https://github.com/klauspost/pgzip](https://github.com/klauspost/pgzip))包来读写gzip文件,输出的gzip文件会略大于GNU gzip生成的文件。
SeqKit 写gzip文件非常快,远快于多线程的pigz,因此无需将结果导入到gzip/pigz中。
自v2.2.0起,SeqKit 也支持读写xz(.xz)和zstd(.zst)格式。自v2.4.0起,支持bzip2格式。
压缩级别:
- 格式 范围 默认值 评论
- gzip 1-9 5 [https://github.com/klauspost/pgzip](https://github.com/klauspost/pgzip) 将5设为默认值。
- xz NA NA [https://github.com/ulikunitz/xz](https://github.com/ulikunitz/xz) 不支持。
- zstd 1-4 2 大致等同于zstd 1、3、7、11。
- bzip 1-9 6 [https://github.com/dsnet/compress](https://github.com/dsnet/compress)
4.1 seq
用法
转换序列(提取ID,按长度过滤,移除间隙,反向互补...)
用法:
seqkit seq [标志]
标志:
-k, --color 序列着色 - 需要通过 "less -R" 来查看
-p, --complement 补全序列,建议开启 '-v' 标志
--dna2rna DNA转RNA
-G, --gap-letters string 用于-g/--remove-gaps的间隙字母(默认 "- \t.")
-h, --help seq的帮助
-l, --lower-case 以小写字母打印序列
-M, --max-len int 只打印长度小于或等于最大长度的序列(-1表示无限制)(默认-1)
-R, --max-qual float 只打印平均质量小于此限制的序列(-1表示无限制)(默认-1)
-m, --min-len int 只打印长度大于或等于最小长度的序列(-1表示无限制)(默认-1)
-Q, --min-qual float 只打印平均质量大于或等于此限制的序列(-1表示无限制)(默认-1)
-n, --name 只打印名称/序列头部
-i, --only-id 打印ID而非完整头部
-q, --qual 只打印质量值
-b, --qual-ascii-base int ASCII基数,33表示Phred+33(默认33)
-g, --remove-gaps 移除由-G/--gap-letters设置的间隙字母,例如空格、制表符和破折号(对齐序列中的间隙"-")
-r, --reverse 反向序列
--rna2dna RNA转DNA
-s, --seq 只打印序列
-u, --upper-case 以大写字母打印序列
-v, --validate-seq 根据字母表验证碱基
-V, --validate-seq-length int 验证序列的长度(0表示整个序列)(默认10000)
4.2 subseq
用法
通过区域/gtf/bed获取子序列,包括侧翼序列。
注意事项:
1. 使用 "seqkit grep" 来提取序列的子集。
"seqtk subseq seqs.fasta id.txt" 等同于
"seqkit grep -f id.txt seqs.fasta"
推荐:
1. 使用纯FASTA文件,这样seqkit可以利用FASTA索引。
2. 推荐使用 -U/--update-faidx 标志,以确保 .fai 文件与FASTA文件匹配。
区域的定义是基于1的,并有一些自定义设计。
示例:
1-基索引 1 2 3 4 5 6 7 8 9 10
负数索引 0-9-8-7-6-5-4-3-2-1
序列 A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
用法:
seqkit subseq [标志]
标志:
--bed string 通过制表符分隔的BED文件
--chr strings 当使用 --gtf 或 --bed 时,选择有限的序列ID(支持多个值,忽略大小写)
-d, --down-stream int 下游长度
--feature strings 选择有限的特征类型(支持多个值,忽略大小写,仅与GTF一起使用)
--gtf string 通过GTF(版本2.2)文件
--gtf-tag string 输出此标签作为序列注释(默认 "gene_id")
-h, --help subseq的帮助
-f, --only-flank 只返回上/下游序列
-r, --region string 通过区域。例如1:12表示前12个碱基,-12:-1表示最后12个碱基,13:-1表示切除前12个碱基。键入 "seqkit subseq -h" 了解更多示例
-R, --region-coord 为 -r/--region 追加序列ID的坐标
-u, --up-stream int 上游长度
-U, --update-faidx 如果存在,则更新fasta索引文件。如果您不确定fasta文件是否更改,请使用此选项
4.3 sliding
用法
在滑动窗口中提取子序列
用法:
seqkit sliding [标志]
标志:
-c, --circular 循环基因组(与 -C/--circular-genome 相同)
-C, --circular-genome 循环基因组(与 -c/--circular 相同)
-g, --greedy 贪婪模式,即,即使窗口大小小于最后的子序列也导出
-h, --help sliding的帮助
-s, --step int 步长大小
-S, --sufix string 添加到序列ID的后缀
-W, --window int 窗口大小
4.4 stats
用法
FASTA/Q文件的简单统计
列:
1. file 输入文件,"-" 表示标准输入(STDIN)
2. format FASTA或FASTQ
3. type DNA, RNA, 蛋白质或无限制
4. num_seqs 序列数量
5. sum_len 基础或残基数量,计算间隙或空格
6. min_len 最小序列长度,计算间隙或空格
7. avg_len 平均序列长度,计算间隙或空格
8. max_len 最大序列长度,计算间隙或空格
9. Q1 序列长度的第一四分位数,计算间隙或空格
10. Q2 序列长度的中位数,计算间隙或空格
11. Q3 序列长度的第三四分位数,计算间隙或空格
12. sum_gap 间隙数量
13. N50 N50. https://en.wikipedia.org/wiki/N50,_L50,_and_related_statistics#N50
14. Q20(%) 质量得分大于20的基础百分比
15. Q30(%) 质量得分大于30的基础百分比
16. AvgQual 平均质量
17. GC(%) GC含量的百分比
注意事项:
1. 序列长度指标(sum_len, min_len, avg_len, max_len, Q1, Q2, Q3)
计算间隙或空格数量。你可以用 "seqkit seq -g" 移除它们:
seqkit seq -g input.fasta | seqkit stats
提示:
1. 对于许多小文件(特别是在SSD上),使用大的 '-j' 值来并行计数。
2. 使用csvtk(https://github.com/shenwei356/csvtk)提取一个指标:
seqkit stats -Ta input.fastq.gz | csvtk cut -t -f "Q30(%)" | csvtk del-header
用法:
seqkit stats [标志]
别名:
stats, stat
标志:
-N, --N strings 其他类似N50的统计。值范围[0, 100],支持多个值,
例如,-N 50,90 或 -N 50 -N 90
-a, --all 所有统计,包括序列长度的四分位数、sum_gap、N50
-b, --basename 只输出文件的基本名
-E, --fq-encoding string fastq质量编码。可用值: 'sanger', 'solexa',
'illumina-1.3+', 'illumina-1.5+', 'illumina-1.8+'。 (默认 "sanger")
-G, --gap-letters string 间隙字母(默认 "- .")
-h, --help stats的帮助
-e, --skip-err 跳过错误,只显示警告消息
-S, --skip-file-check 在给定文件或文件列表时跳过输入文件检查。
-i, --stdin-label string 用于替换标准输入默认"-"的标签(默认"-")
-T, --tabular 以机器友好的表格格式输出
4.5 sum
用法
计算FASTA/Q文件中所有序列的消息摘要
注意事项:
1. 序列头和质量被跳过,只有序列本身重要。
2. 序列记录的顺序无关紧要。
3. 循环完整基因组支持使用 -c/--circular 标志。
- 相同的双链基因组,即使起始位置不同或在反向互补链上,也不会影响结果。
- 对于单链基因组,如ssRNA基因组,使用 -s/--single-strand。
- 消息摘要会随着k-mer大小的不同而改变。
4. 多个文件并行处理(-j/--threads)。
方法:
1. 将序列转换为小写,可选地移除间隙(-g)。
2. 对所有序列或循环完整基因组的k-mer计算哈希(xxhash)(-c/--circular)。
3. 对所有哈希值排序,以忽略序列的顺序。
4. 从哈希值、序列长度和序列数量计算MD5摘要。
遵循Poly中的seqhash(https://github.com/TimothyStiles/poly/),
我们添加元信息到消息摘要中,格式为:
seqkit.<版本>_<序列类型><序列结构><链>_<kmer大小>_<序列摘要>
<版本>: 摘要版本
<序列类型>: 'D' 表示DNA,'R' 表示RNA,'P' 表示蛋白质,'N' 表示其他
<序列结构>: 'L' 表示线性序列,'C' 表示循环基因组
<链>: 'D' 表示双链,'S' 表示单链
<kmer大小>: 0 表示线性序列,循环基因组的其他值
示例:
seqkit.v0.1_DLS_k0_176250c8d1cde6c385397df525aa1a94 DNA.fq.gz
seqkit.v0.1_PLS_k0_c244954e4960dd2a1409cd8ee53d92b9 蛋白质.fasta
seqkit.v0.1_RLS_k0_0f1fb263f0c05a259ae179a61a80578d 单链RNA.fasta
seqkit.v0.1_DCD_k31_e59dad6d561f1f1f28ebf185c6f4c183 双链循环DNA.fasta
seqkit.v0.1_DCS_k31_dd050490cd62ea5f94d73d4d636b7d60 单链循环DNA.fasta
用法:
seqkit sum [标志]
标志:
-a, --all 显示所有信息,包括序列长度和序列数量
-b, --basename 只输出文件的基本名
-c, --circular 文件包含一个循环的基因组序列
-G, --gap-letters string 使用 -g/--remove-gaps 标志删除的间隙字母(默认 "- \t.*")
-h, --help sum的帮助
-k, --kmer-size int 处理循环基因组的k-mer大小(默认1000)
-g, --remove-gaps 删除在选项 -G/gap-letters 中设置的间隙字符
--rna2dna 将RNA转换为DNA
-s, --single-strand 只考虑循环基因组的正链,例如ssRNA病毒基因组
4.6 faidx
用法
创建FASTA索引文件并提取子序列
这个命令类似于 "samtools faidx",但增加了一些额外的功能:
1. 使用标志 -f 输出完整的头行
2. 支持使用标志 -r 将序列ID作为正则表达式
3. 如果你有大量的ID,你可以使用:
seqkit faidx seqs.fasta -l IDs.txt
注意事项:
1. 推荐使用 -U/--update-faidx 标志以确保 .fai 文件与FASTA文件匹配。
区域的定义是基于1的,并有一些自定义设计。
示例:
1-基索引 1 2 3 4 5 6 7 8 9 10
负数索引 0-9-8-7-6-5-4-3-2-1
序列 A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
用法:
seqkit faidx [标志] <fasta文件> [区域...]
标志:
-f, --full-head 打印完整的头行,而不仅仅是ID。将创建以 .seqkit.fai 结尾的新fasta索引文件
-h, --help faidx的帮助
-i, --ignore-case 忽略大小写
-I, --immediate-output 立即打印输出,不使用写缓冲
-l, --region-file string 包含区域列表的文件
-U, --update-faidx 如果存在,则更新fasta索引文件。如果您不确定fasta文件是否更改,请使用此选项
-r, --use-regexp ID是正则表达式。但这里不支持子序列区域。
4.7 watch
用法
监控和在线直方图显示序列特征
用法:
seqkit watch [标志]
标志:
-B, --bins int 直方图的区间数(默认 -1)
-W, --delay int 在线绘图后暂停多少秒(默认 1)
-y, --dump 将直方图数据打印到stderr而不是绘图
-f, --fields string 目标字段,可用值:ReadLen, MeanQual, GC, GCSkew
(默认 "ReadLen")
-h, --help watch的帮助
-O, --img string 将直方图保存到此PDF/图片文件
-H, --list-fields 列出可用字段
-L, --log 对数值进行log10(x+1)变换
-x, --pass 透传模式(将输入写入stdout)
-p, --print-freq int 每处理这么多记录后打印/报告(-1为在EOF后打印)(默认 -1)
-b, --qual-ascii-base int ASCII基数,33表示Phred+33(默认33)
-Q, --quiet-mode 禁止所有绘制到stderr
-R, --reset 每次报告后重置直方图
-v, --validate-seq 根据字母表验证碱基
-V, --validate-seq-length int 验证序列的长度(0表示整个序列)(默认10000)
4.8 sana
用法
清理损坏的单行FASTQ文件
用法:
seqkit sana [标志]
标志:
-A, --allow-gaps 允许序列中有间隙字符(-)
-i, --format string 输入和输出格式:fastq或fasta(默认 "fastq")
-h, --help sana的帮助
-I, --in-format string 输入格式:fastq或fasta
-O, --out-format string 输出格式:fastq或fasta
-b, --qual-ascii-base int ASCII基数,33表示Phred+33(默认33)
4.9 scat
用法
实时递归连接和流式传输fastx文件
用法:
seqkit scat [标志]
标志:
-A, --allow-gaps 允许序列中有间隙字符(-)
-d, --delta int 触发解析的最小大小增加,以千字节为单位(默认 5)
-D, --drop-time string 通知丢弃间隔(默认 "500ms")
-f, --find-only 连接现有文件并退出
-i, --format string 输入和输出格式:fastq或fasta(fastq)(默认 "fastq")
-g, --gz-only 仅查找gzip压缩的文件(.gz后缀)
-h, --help scat的帮助
-I, --in-format string 输入格式:fastq或fasta(fastq)
-O, --out-format string 输出格式:fastq或fasta
-b, --qual-ascii-base int ASCII基数,33表示Phred+33(默认33)
-r, --regexp string 观察文件的正则表达式,默认根据输入格式猜测
-T, --time-limit string 在这个时间段内无活动后退出
-p, --wait-pid int 在此PID的进程退出后(默认 -1)
4.10 fq2fa
Usage
convert FASTQ to FASTA
Usage:
seqkit fq2fa [flags]
4.11 fa2fq
用法
通过FASTA文件检索相应的FASTQ记录
注意事项:
1. 我们假设FASTA文件来源于FASTQ文件,
因此它们共享序列ID,且FASTA中的序列
应该是FASTQ文件中序列的子序列。
用法:
seqkit fa2fq [标志]
标志:
-f, --fasta-file string FASTA文件
-h, --help fa2fq的帮助
-P, --only-positive-strand 仅在正链上搜索
4.12 fx2tab & tab2fx
用法 (fx2tab)
将FASTA/Q转换为表格格式,并提供各种信息,
如序列长度、GC含量/GC偏移。
注意:
1. 对于FASTA或FASTQ,输出固定的三列(ID, 序列, 质量),除非开启了 -n/--name 标志。这是为了格式兼容性。
用法:
seqkit fx2tab [标志]
标志:
-a, --alphabet 打印字母表
-q, --avg-qual 打印读取的平均质量
-B, --base-content strings 打印碱基含量。(忽略大小写,支持多个值)例如:-B AT -B N
-C, --base-count strings 打印碱基计数。(忽略大小写,支持多个值)例如:-C AT -C N
-I, --case-sensitive 计算区分大小写的碱基含量/序列哈希
-g, --gc 打印GC含量
-G, --gc-skew 打印GC偏移
-H, --header-line 打印标题行
-h, --help fx2tab的帮助
-l, --length 打印序列长度
-n, --name 仅打印名称(不包括序列和质量)
-Q, --no-qual 即使对于FASTQ文件也只输出两列
-i, --only-id 打印ID而不是完整头部
-b, --qual-ascii-base int ASCII基数,33表示Phred+33(默认33)
-s, --seq-hash 打印序列的哈希(MD5)
用法 (tab2fx)
将表格格式(前两/三列)转换为FASTA/Q格式
用法:
seqkit tab2fx [标志]
标志:
-b, --buffer-size string 缓冲区大小,支持的单位:K, M, G。当报告“bufio.Scanner: token too long”错误时,需要增加此值(默认 "1G")
-p, --comment-line-prefix strings 注释行前缀(默认 [#,//])
-h, --help tab2fx的帮助
4.13 convert
用法
在Sanger、Solexa和Illumina之间转换FASTQ质量编码
用法:
seqkit convert [标志]
标志:
-d, --dry-run 试运行
-f, --force 对于Illumina-1.8+ -> Sanger,将分数> 40截断为40
--from string 源质量编码。如果未给出,我们将猜测它
-h, --help convert的帮助
-n, --nrecords int 猜测质量编码的记录数(默认1000)
-N, --thresh-B-in-n-most-common int 在N个最常见的质量中'B'的阈值,用于猜测Illumina 1.5。(默认2)
-F, --thresh-illumina1.5-frac float 在前N条记录中Illumina 1.5比例的阈值(默认0.1)
--to string 目标质量编码(默认"Sanger")
4.14 translate
用法
将DNA/RNA翻译成蛋白质序列(支持含有模糊碱基的密码子)
注意:
1. 该命令支持包含任何模糊碱基的密码子。
请开启标志 -L INT 查看详情。例如,对于标准表:
ACN -> T
CCN -> P
CGN -> R
CTN -> L
GCN -> A
GGN -> G
GTN -> V
TCN -> S
MGR -> R
YTR -> L
翻译表/遗传密码:
# https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes
1: 标准密码
2: 脊椎动物线粒体密码
3: 酵母线粒体密码
4: 真菌、原生动物和腔肠动物线粒体密码以及支原体/螺旋体密码
5: 无脊椎动物线粒体密码
6: 纤毛虫、鞭毛藻和六鞭毛虫核密码
9: 棘皮动物和扁形动物线粒体密码
10: 纤毛虫核密码
11: 细菌、古细菌和植物质体密码
12: 另类酵母核密码
13: 海鞘线粒体密码
14: 另类扁形动物线粒体密码
16: 绿藻线粒体密码
21: 吸虫线粒体密码
22: Scenedesmus obliquus线粒体密码
23: Thraustochytrium线粒体密码
24: Pterobranchia线粒体密码
25: 候选分支SR1和Gracilibacteria密码
26: Pachysolen tannophilus核密码
27: Karyorelict核
28: Condylostoma核
29: Mesodinium核
30: Peritrich核
31: Blastocrithidia核
用法:
seqkit translate [标志]
标志:
-x, --allow-unknown-codon 将未知代码翻译为'X'。并且你可能不会使用标志 --trim,它会移除'X'
-F, --append-frame 在序列ID中添加帧信息
--clean 将所有终止密码子位置从'*'字符更改为'X'(未知残基)
-f, --frame strings 翻译的帧,可用值:1, 2, 3, -1, -2, -3, 以及6代表所有六帧(默认 [1])
-h, --help translate的帮助
-M, --init-codon-as-M 将起始处的初始密码子翻译为'M'
-l, --list-transl-table int 显示翻译表N的详情,0代表所有(默认 -1)
-L, --list-transl-table-with-amb-codons int 显示翻译表N的详情(包括模糊密码子),0代表所有。(默认 -1)
-m, --min-len int 氨基酸序列的最小长度
-s, --out-subseqs 输出由终止符号"*"分隔的单独氨基酸子序列
-T, --transl-table int 翻译表/遗传密码,输入 'seqkit translate --help' 查看更多详情(默认 1)
--trim 从翻译的右端移除所有'X'和'*'字符
4.15 grep
用法
通过ID/名称/序列/序列基序搜索序列,允许不匹配
注意事项:
0. 默认情况下,我们用模式匹配序列ID,使用 "-n/--by-name" 以匹配完整名称而不仅是ID。
1. 与POSIX/GNU grep不同,默认情况下我们将模式与整个目标(ID/完整头部)进行比较。请开启 "-r/--use-regexp" 进行部分匹配。
2. 按序列搜索时,是部分匹配,并且会搜索正负两个链。如果您只想在正链上搜索,请开启 "-P/--only-positive-strand"。使用标志 "-m/--max-mismatch" 允许不匹配,您可以增加 "-j/--threads" 的值以加快处理速度。
3. 标志 -d 也支持模糊碱基/残基,如 "RYMM.."。但不要在正则表达式中使用模糊碱基/残基,您需要将它们转换为正则表达式,例如,将 "N" 或 "X" 更改为 "."。
4. 通过标志 '-p' 提供搜索模式(基序)时,请对包含逗号的模式使用双引号,例如,-p '"A{2,}"' 或 -p "\"A{2,}\""。因为命令行参数解析器接受逗号分隔值(CSV)作为多个值(基序)。文件中的模式不遵循此规则。
5. 结果中的序列顺序与原始文件中的顺序一致,而不是查询模式的顺序。但对于FASTA文件,您可以使用:
seqkit faidx seqs.fasta --infile-list IDs.txt
6. 对于多个模式,您可以多次设置 "-p",即 -p pattern1 -p pattern2,或通过 "-f/--pattern-file" 给出模式文件。
您可以使用标志 -R (--region) 指定搜索序列的区域。区域的定义是基于1的,并有一些自定义设计。
示例:
1-基索引 1 2 3 4 5 6 7 8 9 10
负数索引 0-9-8-7-6-5-4-3-2-1
序列 A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
用法:
seqkit grep [标志]
标志:
-D, --allow-duplicated-patterns 当给出重复模式时,输出记录多次
-n, --by-name 通过完整名称而不是仅ID匹配
-s, --by-seq 在序列上搜索子序列。默认情况下,正负两个链都会被搜索,您可能会使用 -P/--only-positive-strand。使用标志 -m/--max-mismatch 允许不匹配
-c, --circular 循环基因组
-C, --count 仅打印匹配记录的数量。使用 -v/--invert-match 标志,计数不匹配的记录
-d, --degenerate 模式/基序包含模糊碱基
--delete-matched 匹配后立即删除一个模式,这样做可以保留首次匹配的数据并在使用正则表达式时加快速度
-h, --help grep的帮助
-i, --ignore-case 忽略大小写
-I, --immediate-output 立即打印输出,不使用写缓冲
-v, --invert-match 反转匹配意义,以选择不匹配的记录
-m, --max-mismatch int 按序列匹配时的最大不匹配数。对于像人类基因组这样的大基因组,使用映射/比对工具会更快
-P, --only-positive-strand 仅在正链上搜索
-p, --pattern strings 搜索模式(支持多个值。注意:对于包含逗号的模式,请使用双引号,例如,-p '"A{2,}"')
-f, --pattern-file string 模式文件(每行一个记录)
-R, --region string 指定搜索序列的区域。例如1:12表示前12个碱基,-12:-1表示最后12个碱基
-r, --use-regexp 模式是正则表达式
4.17 locate
用法
定位子序列/基序,允许不匹配
注意事项:
1. 基序可以是包含 "ACTGN" 的普通序列,或者是像 "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)" 用于开放阅读框的正则表达式。
2. 标志 -d 也支持模糊碱基/残基,如 "RYMM.."。但不要在正则表达式中使用模糊碱基/残基,您需要将它们转换为正则表达式,例如,将 "N" 或 "X" 更改为 "."。
3. 通过标志 '-p' 提供搜索模式(基序)时,请对包含逗号的模式使用双引号,例如,-p '"A{2,}"' 或 -p "\"A{2,}\""。因为命令行参数解析器接受逗号分隔值(CSV)作为多个值(基序)。文件中的模式不遵循此规则。
4. 使用标志 "-m/--max-mismatch" 允许不匹配,您可以增加 "-j/--threads" 的值以加快处理速度。
5. 使用标志 --circular 时,跨越基因组序列末端的匹配子序列的结束位置会大于序列长度。
用法:
seqkit locate [标志]
标志:
--bed 以BED6格式输出
-c, --circular 循环基因组。输入 "seqkit locate -h" 查看详情
-d, --degenerate 模式/基序包含模糊碱基
--gtf 以GTF格式输出
-h, --help locate的帮助
-M, --hide-matched 不显示匹配的序列
-i, --ignore-case 忽略大小写
-I, --immediate-output 立即打印输出,不使用写缓冲
-s, --max-len-to-show int 对搜索模式或匹配的序列,最多显示X个字符
-m, --max-mismatch int 按序列匹配时的最大不匹配数。对于像人类基因组这样的大基因组,使用映射/比对工具会更快
-G, --non-greedy 非贪婪模式,更快但可能错过与其他基序重叠的基序
-P, --only-positive-strand 仅在正链上搜索
-p, --pattern strings 模式/基序(支持多个值。注意:对于包含逗号的模式,请使用双引号,例如,-p '"A{2,}"')
-f, --pattern-file string 模式/基序文件(FASTA格式)
-F, --use-fmi 使用FM索引,对于大量序列模式的搜索会快得多
-r, --use-regexp 模式/基序是正则表达式
-V, --validate-seq-length int 验证序列的长度(0表示整个序列)(默认10000)
4.18 fish
用法
使用局部对齐在较大序列中查找短序列
注意:
1. 输出坐标是类BED格式的0基,左闭右开。
2. 对齐信息打印到STDERR。
用法:
seqkit fish [标志]
标志:
-a, --all 搜索所有
-p, --aln-params string 对齐参数格式 "<匹配>,<不匹配>,<间隙开放>,<间隙扩展>"(默认 "4,-4,-2,-1")
-h, --help fish的帮助
-i, --invert 打印出与任何查询不匹配的参考序列
-q, --min-qual float 最小映射质量(默认 5)
-b, --out-bam string 将对齐保存到此BAM文件(占用大量内存)
-x, --pass 透传模式(将输入写入stdout)
-g, --print-aln 打印序列对齐
-D, --print-desc 打印完整序列头部
-f, --query-fastx string 查询fasta
-F, --query-sequences string 查询序列
-r, --ranges string 目标范围,例如:":10,30:40,-20:"
-s, --stranded 仅搜索+链
-v, --validate-seq 根据字母表验证碱基
-V, --validate-seq-length int 验证序列的长度(0表示整个序列)(默认10000)
4.19 amplicon
用法
通过引物提取扩增子(或其周围特定区域)。
注意事项:
1. 每对引物只返回一个(最长的)匹配位置。
2. 允许不匹配,但不控制不匹配的位置(5'或3')。您可以增加 "-j/--threads" 的值以加快处理速度。
3. 支持模糊碱基/残基,如 "RYMM.."。但不要在正则表达式中使用模糊碱基/残基,您需要将它们转换为正则表达式,例如,将 "N" 或 "X" 更改为 "."。
示例:
0. 没有给定区域。
F
-----===============-----
F R
-----=====-----=====-----
=============== 扩增子
1. 内部区域 (-r x:y)。
F
-----===============-----
1 3 5 x/y
-5-3-1 x/y
F R
-----=====-----=====----- x:y
=============== 1:-1
======= 1:7
===== 3:7
===== 6:10
===== -10:-6
===== -7:-3
-x:y(无效)
2. 侧翼区域 (-r x:y -f)
F
-----===============-----
-3-1 x/y
1 3 5 x/y
F R
-----=====-----=====-----
===== -5:-1
=== -5:-3
===== 1:5
=== 3:5
================= -1:1
========================= -5:5
x:-y(无效)
用法:
seqkit amplicon [标志]
标志:
--bed 以BED6+1格式输出,扩增子作为第7列
-f, --flanking-region 区域是侧翼区域
-F, --forward string 正向引物(5'-引物-3'),允许模糊碱基
-h, --help amplicon的帮助
-I, --immediate-output 立即打印输出,不使用写缓冲
-m, --max-mismatch int 匹配引物时的最大不匹配数,不允许模糊碱基
-P, --only-positive-strand 仅在正链上搜索
-M, --output-mismatches 附加总不匹配数及5'端和3'端的不匹配数
-p, --primer-file string 3列或2列的制表引物文件,第一列为引物名称
-r, --region string 指定要返回的区域。输入 "seqkit amplicon -h" 查看详情
-R, --reverse string 反向引物(5'-引物-3'),允许模糊碱基
-u, --save-unmatched 同时保存未匹配任何引物的记录
-s, --strict-mode 严格模式,即丢弃未完全匹配(较短)给定区域范围的序列
4.20 duplicate
用法
将序列复制N次
您可能需要使用 "seqkit rename" 来使序列ID唯一。
用法:
seqkit duplicate [标志]
别名:
duplicate, dup
标志:
-h, --help duplicate的帮助
-n, --times int 复制次数(默认为1)
4.21 rmdup
用法
通过ID/名称/序列移除重复序列
注意事项:
1. 按序列比较时,正负两个链都会被比较。只考虑正链,请开启 -P/--only-positive-strand。
2. 对于重复的记录,只保存第一条记录。
用法:
seqkit rmdup [标志]
标志:
-n, --by-name 通过完整名称而不仅是id
-s, --by-seq 通过序列
-D, --dup-num-file string 保存重复序列的数量和ID列表的文件
-d, --dup-seqs-file string 保存重复序列的文件
-h, --help rmdup的帮助
-i, --ignore-case 忽略大小写
-P, --only-positive-strand 按序列比较时仅考虑正链
4.22 common
用法
通过id/name/sequence找出多个文件的共有序列
注意:
1. 'seqkit common' 设计用于支持2个及以上的文件。
2. 按序列比较时,
a) 正负两个链都会被比较。您可以开启 -P/--only-positive-strand 仅考虑正链。
b) 您可以开启 -e/--check-embedded-seqs 来检查嵌入的序列。
例如,对于文件A和B,文件B中CCCC的反向互补序列是文件A中TTGGGGTT的一部分,我们将从文件A中提取并输出GGGG。
如果序列CCC存在于除文件A外的其他文件中,我们将跳过它,因为它是GGGG的嵌入子序列。
建议将最小的文件作为第一个文件,以节省内存使用。
3. 对于2个文件,'seqkit grep' 更快且消耗更少内存:
seqkit grep -f <(seqkit seq -n -i small.fq.gz) big.fq.gz # 按序列ID
但注意,按序列搜索会慢得多,因为它是部分字符串匹配。
seqkit grep -s -f <(seqkit seq -s small.fq.gz) big.fq.gz # 慢得多!!!!
4. 一个文件中的某些记录可能有相同的序列/ID。如果序列/ID在多个文件中共享,它们将全部被检索。
因此,记录数量可能大于最小文件的记录数量。
用法:
seqkit common [标志]
标志:
-n, --by-name 通过完整名称而不仅是id匹配
-s, --by-seq 通过序列匹配
-e, --check-embedded-seqs 检查嵌入的序列,例如,如果一个序列是另一个序列的一部分,我们将保留较短的那个
-h, --help common的帮助
-i, --ignore-case 忽略大小写
-P, --only-positive-strand 按序列比较时仅考虑正链
4.23 split
用法
根据名称ID、给定区域的子序列、部分大小或部分数量将序列分割成文件。
如果您只想按部分或大小分割,请使用 "seqkit split2",
它可以应用于成对和单端FASTQ。
如果您想将一个序列切割成多个片段。
1. 要均匀切割,请使用 'kmcp utils split-genomes'
(https://bioinf.shenwei.me/kmcp/usage/#split-genomes)。
例如,均匀切割成4个等大小的片段,相邻片段之间无重叠:
kmcp utils split-genomes -m 1 -k 1 --split-number 4 --split-overlap 0 input.fasta -O out
2. 要切割成多个固定大小的片段,请使用 'seqkit sliding'。
例如,切割成40 bp的片段,并保留最后一个可能短于40 bp的片段。
seqkit sliding -g -s 40 -W 40 input.fasta -o out.fasta
注意事项:
1. 对于两遍模式 (-2/--two-pass),建议使用标志 -U/--update-faidx
以确保 .fai 文件与FASTA文件匹配。
区域的定义是基于1的,并有一些自定义设计。
示例:
1-基索引 1 2 3 4 5 6 7 8 9 10
负数索引 0-9-8-7-6-5-4-3-2-1
序列 A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
用法:
seqkit split [标志]
标志:
-i, --by-id 根据序列ID分割序列
--by-id-prefix string --by-id的文件前缀
-p, --by-part int 将序列分割成N部分
--by-part-prefix string --by-part的文件前缀
-r, --by-region string 根据给定区域的子序列分割序列。例如1:12表示前12个碱基,-12:-1表示最后12个碱基。输入 "seqkit split -h" 查看更多示例
--by-region-prefix string --by-region的文件前缀
-s, --by-size int 将序列分割成多部分,每部分N个序列
--by-size-prefix string --by-size的文件前缀
-d, --dry-run 试运行,只打印消息,不创建文件。
-e, --extension string 设置输出文件扩展名,例如 ".gz", ".xz", 或 ".zst"
-f, --force 覆盖输出目录
-h, --help split的帮助
-k, --keep-temp 使用两遍模式时保留临时FASTA和.fai文件(仅适用于FASTA格式)
-O, --out-dir string 输出目录(默认值是 $infile.split)
-2, --two-pass 两遍模式读取文件两次以降低内存使用。(仅适用于FASTA格式)
-U, --update-faidx 如果存在,则更新fasta索引文件。如果您不确定fasta文件是否更改,请使用此选项
4.24 split2 比split更好用!
用法
根据部分大小或部分数量将序列分割成文件
此命令支持FASTA和成对或单端FASTQ,具有低内存占用和快速速度。
输出文件的前缀:
1. 对于stdin:stdin
2. 其他:与输入文件相同
3. 通过选项设置:--by-length-prefix, --by-part-prefix, 或 --by-size-prefix
输出文件的扩展名:
1. 对于stdin:.fast[aq]
2. 其他:与输入文件相同
3. 通过选项 -e/--extension 添加额外扩展名,例如,为纯文本输入输出gzip压缩文件:
seqkit split2 -p 2 -O test tests/hairpin.fa -e .gz
如果您想将一个序列切割成多个片段。
1. 要均匀切割,请使用 'kmcp utils split-genomes'
(https://bioinf.shenwei.me/kmcp/usage/#split-genomes)。
例如,均匀切割成4个等大小的片段,相邻片段之间无重叠:
kmcp utils split-genomes -m 1 -k 1 --split-number 4 --split-overlap 0 input.fasta -O out
2. 要切割成多个固定大小的片段,请使用 'seqkit sliding'。
例如,切割成40 bp的片段,并保留最后一个可能短于40 bp的片段。
seqkit sliding -g -s 40 -W 40 input.fasta -o out.fasta
用法:
seqkit split2 [标志]
标志:
-l, --by-length string 将序列分割成 >=N 基的块,支持 K/M/G 后缀
--by-length-prefix string --by-length的文件前缀
-p, --by-part int 将序列分割成N部分
--by-part-prefix string --by-part的文件前缀
-s, --by-size int 将序列分割成多部分,每部分N个序列
--by-size-prefix string --by-size的文件前缀
-e, --extension string 设置输出文件扩展名,例如 ".gz", ".xz", 或 ".zst"
-f, --force 覆盖输出目录
-h, --help split2的帮助
-O, --out-dir string 输出目录(默认值是 $infile.split)
-1, --read1 string (gzip压缩的) read1文件
-2, --read2 string (gzip压缩的) read2文件
4.25 pair
用法
从两个fastq文件中匹配配对末端读取
注意事项:
1. 两个文件中头部的顺序最好相同(未被打乱),否则,它会消耗大量内存来缓存内存中的读取。
2. 不配对的读取是可选输出,使用标志 -u/--save-unpaired。
3. 如果没有给出标志 -O/--out-dir,输出将保存在输入的同一目录中,后缀为 "paired",例如,read_1.paired.fq.gz。否则,在给定的输出目录中保持名称不变。
4. 配对的gzip压缩文件可能比原始文件略大,因为使用了不同的gzip包/库,不用担心。
提示:
1. 支持MGI等平台生成的配对读取文件的 '/1' 和 '/2' 标签。您可以简单地指定用于提取序列ID的正则表达式:
--id-regexp '^(\S+)\/[12]'
用法:
seqkit pair [标志]
标志:
-f, --force 覆盖输出目录
-h, --help pair的帮助
-O, --out-dir string 输出目录
-1, --read1 string (gzip压缩的) read1文件
-2, --read2 string (gzip压缩的) read2文件
-u, --save-unpaired 如果存在不配对的读取,则保存
4.26 sample
用法
通过数量或比例抽取序列。
注意事项:
1. 不要在大的FASTQ文件上使用 '-n',它会将所有序列加载到内存中!
相反,使用 'seqkit sample -p 0.1 seqs.fq.gz | seqkit head -n N'!
用法:
seqkit sample [标志]
标志:
-h, --help sample的帮助
-n, --number int 通过数量抽取(结果可能不完全匹配),不要在大的FASTQ文件上使用。
-p, --proportion float 通过比例抽取
-s, --rand-seed int 随机种子。对于成对端数据,使用相同的种子跨fastq文件抽取相同的读取对(默认11)
-2, --two-pass 两遍模式读取文件两次以降低内存使用。当从stdin读取时不允许
4.27 head
Usage
print first N FASTA/Q records
For returning the last N records, use:
seqkit range -r -N:-1 seqs.fasta
Usage:
seqkit head [flags]
Flags:
-n, --number int print first N FASTA/Q records (default 10)
4.28 head-genome
用法
打印名称中有共同前缀的第一个基因组的序列
对于包含多个菌株的contigs的FASTA文件(见下面的例子),
没有可用的ID列表来检索特定菌株的序列,
而每个菌株的描述共享相同的前缀。
此命令用于检索第一个菌株的序列,
即 "Vibrio cholerae strain M29"。
>NZ_JFGR01000001.1 Vibrio cholerae strain M29 Contig_1, whole genome shotgun sequence
>NZ_JFGR01000002.1 Vibrio cholerae strain M29 Contig_2, whole genome shotgun sequence
>NZ_JFGR01000003.1 Vibrio cholerae strain M29 Contig_3, whole genome shotgun sequence
>NZ_JSTP01000001.1 Vibrio cholerae strain 2012HC-12 NODE_79, whole genome shotgun sequence
>NZ_JSTP01000002.1 Vibrio cholerae strain 2012HC-12 NODE_78, whole genome shotgun sequence
注意事项:
1. 文件中的序列应该组织得很好。
用法:
seqkit head-genome [标志]
标志:
-h, --help head-genome的帮助
-m, --mini-common-words int 最小共享前缀词数(默认4)
4.29 range
用法
打印位于范围内的FASTA/Q记录(起始:结束)
示例:
1. 前100条记录(head -n 100)
seqkit range -r 1:100
2. 最后100条记录(tail -n 100)
seqkit range -r -100:-1
3. 移除前100条记录(tail -n +101)
seqkit range -r 101:-1
4. 其他范围:
seqkit range -r 10:100
seqkit range -r -100:-10
用法:
seqkit range [标志]
标志:
-h, --help range的帮助
-r, --range string 范围。例如,1:12表示前12条记录(head -n 12),-12:-1表示最后12条记录(tail -n 12)
4.30 replace
用法
通过正则表达式替换名称/序列。
注意替换支持捕获变量。
例如,$1代表第一个子匹配的文本。
注意:在*nix操作系统中使用单引号而不是双引号。
示例:给所有碱基添加空格。
seqkit replace -p "(.)" -r '$1 ' -s
或使用转义字符 \。
seqkit replace -p "(.)" -r "\$1 " -s
更多信息:http://bioinf.shenwei.me/seqkit/usage/#replace
特殊替换符号(仅适用于替换名称而非序列):
{nr} 记录编号,从1开始
{kv} 键值文件中键(捕获变量$n)对应的值,n可以通过标志 -I (--key-capt-idx) 指定(默认:1)
特殊情况:
1. 如果替换内容包含'$',
a). 如果使用'{kv}',你需要使用'$$$$'代替单个'$':
-r '{kv}' -k <(sed 's/\$/$$$$/' kv.txt)
b). 如果不是,使用'$$':
-r 'xxx$$xx'
筛选记录以进行编辑:
您可以使用类似于 "seqkit grep" 中的标志来选择部分记录进行编辑。
用法:
seqkit replace [标志]
标志:
-s, --by-seq 替换序列(仅FASTA)
--f-by-name [目标筛选] 通过完整名称而不仅是ID匹配
--f-by-seq [目标筛选] 在序列上搜索子序列,正负两个链都会被搜索,并且使用标志 -m/--max-mismatch 允许不匹配
--f-ignore-case [目标筛选] 忽略大小写
--f-invert-match [目标筛选] 反转匹配意义,选择不匹配的记录
--f-only-positive-strand [目标筛选] 仅在正链上搜索
--f-pattern strings [目标筛选] 搜索模式(支持多个值。注意:对于包含逗号的模式,请使用双引号,例如,-p '"A{2,}"')
--f-pattern-file string [目标筛选] 模式文件(每行一个记录)
--f-use-regexp [目标筛选] 模式是正则表达式
-h, --help replace的帮助
-i, --ignore-case 忽略大小写
-K, --keep-key 当键没有对应的值时保留键作为值(仅适用于序列名称)
-U, --keep-untouch 当键没有对应的值时不做任何改变(仅适用于序列名称)
-I, --key-capt-idx int 键的捕获变量索引(基于1)(默认 1)
-m, --key-miss-repl string 没有对应值的键的替换内容
-k, --kv-file string 用于替换键为值的制表符分隔键值文件,当在 -r (--replacement) 中使用 "{kv}" 时(仅适用于序列名称)
--nr-width int 标志 -r/--replacement 中{nr}的最小宽度。例如,通过 --nr-width 3 将 "1" 格式化为 "001"(默认 1)
-p, --pattern string 搜索正则表达式
-r, --replacement string 替换内容。支持捕获变量。例如,$1代表第一个子匹配的文本。注意:对于*nix OS,使用单引号而不是双引号或使用转义字符。当给出{kv}时,使用${1}代替$1!
4.31 rename
用法
重命名重复的ID
注意事项:
1. 该命令仅将 "_N" 附加到重复的序列ID上以使它们唯一。
2. 使用 "seqkit replace" 通过正则表达式编辑序列ID/头部。
示例:
$ seqkit seq seqs.fasta
>id comment
actg
>id description
ACTG
$ seqkit rename seqs.fasta
>id comment
actg
>id_2 description
ACTG
用法:
seqkit rename [标志]
标志:
-n, --by-name 通过完整名称而不仅是id检查重复
-f, --force 覆盖输出目录
-h, --help rename的帮助
-m, --multiple-outfiles 对多个输入文件的结果写入分离的文件
-O, --out-dir string 输出目录(默认为"renamed")
-1, --rename-1st-rec 同样重命名第一条记录
-s, --separator string 原始ID/名称与计数器之间的分隔符(默认为"_")
-N, --start-num int *重复* ID/名称的起始计数数字,应大于零(默认为2)
4.32 restart
用法
重置圆形基因组的起始位置
示例
$ echo -e ">seq\nacgtnACGTN"
>seq
acgtnACGTN
$ echo -e ">seq\nacgtnACGTN" | seqkit restart -i 2
>seq
cgtnACGTNa
$ echo -e ">seq\nacgtnACGTN" | seqkit restart -i -2
>seq
TNacgtnACG
用法:
seqkit restart [标志]
标志:
-i, --new-start int 新的起始位置(1-base,支持从末尾计数的负值)(默认为1)
4.33 concat
用法
将多个文件中具有相同ID的序列连接起来
注意事项:
1. 默认情况下,只有在所有文件中都出现ID的序列才会被输出。
使用 -f/--full 来输出所有序列。
2. 如果同一ID有多个序列,我们会输出序列的笛卡尔积。
3. 描述也会用分隔符连接起来 (-s/--separator)。
4. 不同ID的序列的顺序是随机的。
用法:
seqkit concat [标志]
别名:
concat, concate
标志:
-f, --full 保留所有序列,类似于全连接/外连接
-h, --help concat的帮助
-s, --separator string 具有相同ID的记录的描述分隔符(默认为"|")
4.34 mutate
用法
编辑序列(点突变,插入,删除)
注意事项:
1. 允许多个点突变 (-p/--point),但只允许单个插入 (-i/--insertion) 或单个删除 (-d/--deletion)。
2. 点突变在插入/删除之前发生。
备注:
1. 您可以使用 'seqkit grep' 中的类似标志来选择要编辑的特定序列。
位置的定义是基于1的,并有一些自定义设计。
示例:
1-基索引 1 2 3 4 5 6 7 8 9 10
负数索引 0-9-8-7-6-5-4-3-2-1
序列 A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
用法:
seqkit mutate [标志]
标志:
-n, --by-name [匹配要编辑的序列] 通过完整名称而不仅是id匹配
-d, --deletion string 删除突变:在范围内删除子序列。例如,-d 1:2 为删除前两个碱基,-d -3:-1 删除最后3个碱基
-h, --help mutate的帮助
-I, --ignore-case [匹配要编辑的序列] 忽略搜索模式的大小写
-i, --insertion string 插入突变:在给定位置后插入碱基,例如,-i 0:ACGT 在开头插入ACGT,-1:* 在末尾添加*
-v, --invert-match [匹配要编辑的序列] 反转匹配意义,以选择不匹配的记录
-s, --pattern strings [匹配要编辑的序列] 搜索模式(支持多个值。注意:对于包含逗号的模式,请使用双引号,例如,-p '"A{2,}"')
-f, --pattern-file string [匹配要编辑的序列] 模式文件(每行一个记录)
-p, --point strings 点突变:在给定位置更改碱基。例如,-p 2:C 将第2个碱基设为C,-p -1:A 将最后一个碱基更改为A
-r, --use-regexp [匹配要编辑的序列] 搜索模式是正则表达式
4.35 shuffle
用法
洗牌序列。
默认情况下,所有记录将被读取到内存中。
对于FASTA格式,使用标志 -2 (--two-pass) 以减少内存使用。FASTQ不支持。
首先,seqkit读取序列ID。如果文件不是纯FASTA文件,
seqkit将把序列写入临时文件,并创建FASTA索引。
其次,seqkit洗牌序列ID并通过FASTA索引提取序列。
注意事项:
1. 对于两遍模式 (-2/--two-pass),建议使用标志 -U/--update-faidx
以确保 .fai 文件与FASTA文件匹配。
用法:
seqkit shuffle [标志]
标志:
-h, --help shuffle的帮助
-k, --keep-temp 使用两遍模式时保留临时FASTA和.fai文件
-s, --rand-seed int 洗牌的随机种子(默认为23)
-2, --two-pass 两遍模式读取文件两次以降低内存使用。(仅适用于FASTA格式)
-U, --update-faidx 如果存在,则更新fasta索引文件。如果您不确定fasta文件是否更改,请使用此选项
4.36 sort
用法
根据id/name/sequence/length对序列进行排序。
默认情况下,所有记录将被读取到内存中。
对于FASTA格式,使用标志 -2 (--two-pass) 以减少内存使用。FASTQ不支持。
首先,seqkit读取序列头和长度信息。
如果文件不是纯FASTA文件,
seqkit将把序列写入临时文件,并创建FASTA索引。
其次,seqkit根据头和长度信息对序列进行排序
并通过FASTA索引提取序列。
注意事项:
1. 对于两遍模式 (-2/--two-pass),建议使用标志 -U/--update-faidx
以确保 .fai 文件与FASTA文件匹配。
用法:
seqkit sort [标志]
标志:
-b, --by-bases 按非空位碱基排序
-l, --by-length 按序列长度排序
-n, --by-name 按完整名称而不仅是id排序
-s, --by-seq 按序列排序
-G, --gap-letters string 间隙字母(默认为"- \t.")
-h, --help sort的帮助
-i, --ignore-case 忽略大小写
-k, --keep-temp 使用两遍模式时保留临时FASTA和.fai文件
-N, --natural-order 按自然顺序排序,当按ID/完整名称排序时
-r, --reverse 反转结果
-L, --seq-prefix-length int 序列前缀长度,seqkit按序列排序时基于此长度(0代表整个序列)(默认10000)
-2, --two-pass 两遍模式读取文件两次以降低内存使用。(仅适用于FASTA格式)
-U, --update-faidx 如果存在,则更新fasta索引文件。如果您不确定fasta文件是否更改,请使用此选项
4.37 bam
用法
监控和在线直方图显示BAM记录特征
用法:
seqkit bam [标志]
标志:
-B, --bins int 直方图的区间数(默认 -1)
-N, --bundle int 将BAM文件划分为位点(-1)或以此最小大小的包
-c, --count string 按参考序列计数读取并保存到此文件
-W, --delay int 绘图后暂停这么多秒(默认 1)
-y, --dump 将直方图数据打印到stderr而不是绘图
-G, --exclude-ids string 排除此文件中包含ID的记录
-e, --exec-after string 报告后执行命令
-E, --exec-before string 报告前执行命令
-f, --field string 目标字段
-g, --grep-ids string 仅保留此文件中包含ID的记录
-h, --help bam的帮助
-C, --idx-count 基于BAM索引的快速按参考序列计数
-i, --idx-stat 基于BAM索引的快速统计
-O, --img string 将直方图保存到此PDF/图片文件
-H, --list-fields 列出所有可用的BAM记录特征
-L, --log 对数值进行log10(x+1)变换
-q, --map-qual int 最小映射质量
-x, --pass 透传模式(将过滤后的BAM转发到输出)
-k, --pretty 某些TSV输出的美化打印
-F, --prim-only 过滤掉非主要对齐记录
-p, --print-freq int 每处理这么多记录后打印/报告(-1为在EOF后打印)(默认 -1)
-Q, --quiet-mode 禁止所有绘制到stderr
-M, --range-max float 丢弃字段(-f)值大于此标志的记录(默认 NaN)
-m, --range-min float 丢弃字段(-f)值小于此标志的记录(默认 NaN)
-R, --reset 每次报告后重置直方图
-Z, --silent-mode 禁止TSV输出到stderr
-s, --stat 打印输入文件的BAM统计信息
-T, --tool string 调用YAML格式的工具箱(见文档)
-@, --top-bam string 将前-?记录保存到此bam文件
-?, --top-size int top模式缓冲区的大小(默认 100)
4.38 merge-slides
用法
合并由seqkit滑动窗口生成的滑动窗口
例如,
ref.contig00001_sliding:454531-454680
ref.contig00001_sliding:454561-454710
ref.contig00001_sliding:454591-454740
ref.contig00002_sliding:362281-362430
ref.contig00002_sliding:362311-362460
ref.contig00002_sliding:362341-362490
ref.contig00002_sliding:495991-496140
ref.contig00044_sliding:1-150
ref.contig00044_sliding:31-180
ref.contig00044_sliding:61-210
ref.contig00044_sliding:91-240
可以合并成
ref.contig00001 454530 454740
ref.contig00002 362280 362490
ref.contig00002 495990 496140
ref.contig00044 0 240
输出(BED3格式):
1. chrom - 染色体名称
2. chromStart - 起始位置(0-based)
3. chromEnd - 结束位置(0-based)
用法:
seqkit merge-slides [标志]
标志:
-b, --buffer-size string 缓冲区大小,支持单位:K, M, G。当报告"bufio.Scanner: token too long"错误时,您需要增加此值(默认 "1G")
-p, --comment-line-prefix strings 注释行前缀(默认 [#,//])
-h, --help merge-slides的帮助
-g, --max-gap int 两个相邻区域起始位置的最大距离,0表示无限制,1表示不合并。
-l, --min-overlap int 两个相邻区域的最小重叠,建议 $sliding_step_size - 1.(默认 1)
-r, --regexp string 用于提取参考名称和窗口位置的正则表达式。(默认 "^(.+)_sliding:(\\d+)\\-(\\d+)")
4.39 genautocomplete
用法
生成shell自动补全脚本
支持的shell:bash|zsh|fish|powershell
Bash:
# 生成补全脚本
seqkit genautocomplete --shell bash
# 如果从未配置过,请配置。
# 如果找不到 "complete" 命令,请安装 bash-completion。
echo "for bcfile in ~/.bash_completion.d/* ; do source \$bcfile; done" >> ~/.bash_completion
echo "source ~/.bash_completion" >> ~/.bashrc
Zsh:
# 生成补全脚本
seqkit genautocomplete --shell zsh --file ~/.zfunc/_seqkit
# 如果从未配置过,请配置
echo 'fpath=( ~/.zfunc "${fpath[@]}" )' >> ~/.zshrc
echo "autoload -U compinit; compinit" >> ~/.zshrc
fish:
seqkit genautocomplete --shell fish --file ~/.config/fish/completions/seqkit.fish
用法:
seqkit genautocomplete [标志]
标志:
--file string 自动补全文件(默认 "/home/shenwei/.bash_completion.d/seqkit.sh")
-h, --help genautocomplete的帮助
--type string 自动补全类型(目前仅支持bash)(默认 "bash")
一共39个,但是split两个算一个,还是38个subcommand。
05 常用命令行
依据自己的需求,查看且在linux多尝试即可,因为这个命令真的很多,唯手熟尔!
#Read and print
#From file:
$ seqkit seq hairpin.fa.gz
>cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
$ seqkit seq reads_1.fq.gz
@HWI-D00523:240:HF3WGBCXX:1:1101:2574:2226 1:N:0:CTGTAG
TGAGGAATATTGGTCAATGGGCGCGAGCCTGAACCAGCCAAGTAGCGTGAAGGATGACTGCCCTACGGG
+
HIHIIIIIHIIHGHHIHHIIIIIIIIIIIIIIIHHIIIIIHHIHIIIIIGIHIIIIHHHHHHGHIHIII
From stdin:
zcat hairpin.fa.gz | seqkit seq
#Sequence types
#By default, seqkit seq automatically detect the sequence type
$ echo -e ">seq\nacgtryswkmbdhvACGTRYSWKMBDHV" | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA DNA 1 28 28 28 28
$ echo -e ">seq\nACGUN ACGUN" | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA RNA 1 11 11 11 11
$ echo -e ">seq\nabcdefghijklmnpqrstvwyz" | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA Protein 1 23 23 23 23
$ echo -e "@read\nACTGCN\n+\n@IICCG" | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTQ DNA 1 6 6 6 6
#You can also set sequence type by flag -t (--seq-type). But this only take effect on #subcommands seq and locate.
$ echo -e ">seq\nabcdefghijklmnpqrstvwyz" | seqkit seq -t dna
[INFO] when flag -t (--seq-type) given, flag -v (--validate-seq) is automatically switched on
[ERRO] error when parsing seq: seq (invalid DNAredundant letter: e)
Only print names
Full name:
$ seqkit seq hairpin.fa.gz -n
cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
cel-lin-4 MI0000002 Caenorhabditis elegans lin-4 stem-loop
cel-mir-1 MI0000003 Caenorhabditis elegans miR-1 stem-loop
Only ID:
$ seqkit seq hairpin.fa.gz -n -i
cel-let-7
cel-lin-4
cel-mir-1
#Custom ID region by regular expression (this could be applied to all subcommands):
$ seqkit seq hairpin.fa.gz -n -i --id-regexp "^[^\s]+\s([^\s]+)\s"
MI0000001
MI0000002
MI0000003
#Only print seq (global flag -w defines the output line width, 0 for no wrap)
$ seqkit seq hairpin.fa.gz -s -w 0
UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAACUAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
AUGCUUCCGGCCUGUUCCCUGAGACCUCAAGUGUGAGUGUACUAUUGAUGCUUCACACCUGGGCUCUCCGGGUACCAGGACGGUUUGAGCAGAU
AAAGUGACCGUACCGAGCUGCAUACUUCCUUACAUGCCCAUACUAUAUCAUAAAUGGAUAUGGAAUGUAAAGAAGUAUGUAGAACGGGGUGGUAGU
Convert multi-line FASTQ to 4-line FASTQ
$ seqkit seq reads_1.fq.gz -w 0
#Reverse comlement sequence
$ seqkit seq hairpin.fa.gz -r -p
>cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
UCGAAGAGUUCUGUCUCCGGUAAGGUAGAAAAUUGCAUAGUUCACCGGUGGUAAUAUUCC
AAACUAUACAACCUACUACCUCACCGGAUCCACAGUGUA
#Remove gaps and to lower/upper case
$ echo -e ">seq\nACGT-ACTGC-ACC" | seqkit seq -g -u
>seq
ACGTACTGCACC
RNA to DNA
$ echo -e ">seq\nUCAUAUGCUUGUCUCAAAGAUUA" | seqkit seq --rna2dna
>seq
TCATATGCTTGTCTCAAAGATTA
Filter by sequence length
$ cat hairpin.fa | seqkit seq | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA RNA 28,645 2,949,871 39 103 2,354
$ cat hairpin.fa | seqkit seq -m 100 | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA RNA 10,975 1,565,486 100 142.6 2,354
$ cat hairpin.fa | seqkit seq -m 100 -M 1000 | seqkit stats
file format type num_seqs sum_len min_len avg_len max_len
- FASTA RNA 10,972 1,560,270 100 142.2 938
06 参考文献
W Shen, S Le, Y Li*, F Hu*. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PLOS ONE. doi:10.1371/journal.pone.0163962.
孟倩.基于高通量测序的短序列生物数据压缩研究[J].计算机应用与软件,2017,34(04):22-27+98.