seqkit的使用说明1

Seqkit

承豆

6090人浏览 · 2022-04-05 18:56:10

承豆 · 2022-04-05 18:56:10 发布

不得不吹一个软件，这是我目前见过所有软件最牛的一个软件，其他软件跟它比较都是弱爆了。序列处理的巅峰软件seqkit（个人建议，至少我目前没见过比这个更厉害的软件）

seqkit 的官网网站：https://bioinf.shenwei.me/seqkit/

SeqKit - 用于 FASTA/Q 文件操作的跨平台和超快工具包

特征

易于安装（下载）
- 为多个平台（Linux/Windows/macOS、amd64/arm64）提供可执行二进制文件
- 轻量级开箱即用，无依赖，无需编译，无需配置
使用方便
- 超快
- 无缝解析 FASTA 和 FASTQ 格式
- 支持（gzip/压缩）STDIN xz/ zstdSTDOUT和输入/输出文件，轻松集成在管道中
- sample可重现的结果（和中的可配置 rand 种子shuffle）
- 通过正则表达式支持自定义序列 ID
- 支持Bash/Zsh 补全

多功能命令（用法和示例）

37个子命令支持的实用功能

## 序列和子序列
**seq**  转换序列（序列颠倒，序列互补，提取ID）
**subseq** 从区域/gtf/bed中获得序列，包括侧面的序列
**sliding** 滑动序列，支持环式基因组
**stats**   对FASTA/Q files进行简单统计
**faidx** 创造fasta索引文件并提取子序列
**watch** 检测并连线序列特点的柱状图
**sana** 清除质量不好的单线的fastq文件
## 格式转换
**fx2tab**  将FASTA/Q 文件转变成表格形式 (1th: name/ID, 2nd: sequence, 3rd: quality)
**tab2fx** 转变表格形式为fasta/q格式
**fq2fa** 转变fastq文件为fasta文件
**convert** 在Sanger, Solexa and Illumina中转换fastq的质量编码
**translate** 将DNA/RNA序列转变成蛋白序列（支持模棱两可的碱基）
## 搜索
**grep** 根据ID/名称/序列/序列motif 搜索序列，且允许错配
**locate** 定位子序列/motif,且允许错配
**fish** 使用本地比对在较大序列中寻找短序列
**amplicon** 经由引物检索扩增子（或它附近特定的区域）
## bam文件的处理和监视
**bam** 监视和连线bam文件记录特点的直方图
## 设置参数
**head** 打印第一个Nfasta/q的记录
**range** 在一个范围内（start:end）打印fasta/q的记录
**sample** 通过数量或比例来体验序列
**rmdup** 通过id/名称/序列 来去除复制的序列
**duplicate**  复制N次的序列
**common** 通过id/名称/序列 发现多条序列中共有的序列
**split** 通过id/seq region/size/parts (mainly for FASTA) 将序列劈开成文件
**split2** 将序列通过大小或部分 劈开成文件
## 编辑
**replace** 通过规律表达来代替名字或序列
**rename** 重新命名复制的ID
**restart** 为环状基因组重新设置起始位置
**concat** 从多个文件中经由相同的ID来连接序列
**mutate** 编辑序列（点突，插入，删除）
## 排序
**shuffle** 变换序列位置
**sort** 将序列经由id/name/sequence 进行排序

全局参数

Global Flags:
--alphabet-guess-seq-length int seqkit根据它猜测序列类型的第一个FASTA记录的序列前缀长度(0表示整个seq)(默认10000)
--id-ncbi                  FASTA序列命名格式ncbi格式的>gi|110645304|ref|NC_002516.2| 
--id-regexp string         解析ID的正则表达式(default "^(\\S+)\\s?")
--infile-list string       file of input files list 
-w, --line-width int       输出FASTA格式时的线宽(default 60, 0表示不需要换行)
-o, --out-file string      输出文件 (.gz for gzipped out) (default "-")
--quiet                    不要显示额外的信息
-t, --seq-type string      序列类型 (dna|rna|protein|unlimit|auto)(default "auto", up to the first sequence )
-j, --threads int          CPU个数 (default 4)

安装

# 安装seqkit
conda install seqkit
# 补充bedops(gtf2bed)和csvtk工具(可选)
conda install  bedops -y
conda install  csvtk -y

数据准备

# miRBase的RNA序列1.5M/785K/110K
wget -c https://www.mirbase.org/ftp/CURRENT/hairpin.fa.gz # Fasta格式的所有miRNA发夹序列
wget -c https://www.mirbase.org/ftp/CURRENT/mature.fa.gz # Fasta可格式化所有成熟miRNA序列
wget -c https://www.mirbase.org/ftp/CURRENT/miRNA.diff.gz # 上一个版本和这个版本之间的变化

# 下载拟南芥基因组数据
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.gff.gz

seq 序列处理文件

转换序列(提取ID，按长度过滤，删除间隙)

Flags: 
-k, --color            colorize sequences - to be piped into "less -R" 
-p, --complement       complement sequence, flag '-v' is recommended to switch on 
   --dna2rna           DNA转化为RNA 
-G, --gap-letters string     gap letters (default "- \t.") 
-h, --help             help for seq 
-l, --lower-case       小写打印序列 
-M, --max-len int      只打印小于最大长度的序列 (-1 for no limit) (default -1) 
-R, --max-qual float   只打印平均质量小于这个限制的序列(-1 for no limit) (default -1) 
-m, --min-len int      只打印大于最小长度的序列(-1 for no limit) (default -1) 
-Q, --min-qual float   只打印平均质量qreater或等于这个限制的序列 (default -1) 
-n, --name             只打印序列名称 
-i, --only-id          只打印序列名称的ID
-q, --qual            only print qualities 
-b, --qual-ascii-base int    ASCII BASE, 33 for Phred+33 (default 33) 
-g, --remove-gaps      删除间隙
 -r, --reverse          reverse sequence 
   --rna2dna          RNA转换为DNA 
-s, --seq            只打印序列
-u, --upper-case     大写打印序列
-v, --validate-seq   根据字母表验证基     
-V, --validate-seq-length int  要验证的序列长度(0表示整个seq)(默认10000)

示例

seqkit seq hairpin.fa.gz #序列结果展示
>aae-mir-11912 MI0037955 Aedes aegypti miR-11912 stem-loop
AGCAGUGUGUGAACCGUUGGCGGCGCAACGCCUGAUGGUAUUGAUGCUGCU

seqkit seq hairpin.fa.gz # 仅展示序列名称
gga-mir-1784b MI0041072 Gallus gallus miR-1784b stem-loop
mdo-mir-7385g-1 MI0041073 Monodelphis domestica miR-7385g-1 stem-loop
mdo-mir-7385g-2 MI0041074 Monodelphis domestica miR-7385g-2 stem-loop

seqkit seq hairpin.fa.gz -n -i  #展示序列ID
smc-mir-12460
smc-mir-12461
mdo-mir-7385g-2

seqkit seq hairpin.fa.gz -n -i --id-regexp "^[^\s]+\s([^\s]+)\s"  #使用正则匹配序列名
MI0000001
MI0000002
MI0000003

# RNA转化为DNA并输出 
seqkit seq hairpin.fa.gz --rna2dna > hairpinDNA.fasta 
# 反转录序列
seqkit seq hairpin.fa.gz -r -p  
# 去除序列gap 并大写碱基
echo -e ">seq\nACGT-actgc-ACC" | seqkit seq -g -u

subseq

**目的：**按区域/gtf/床获取子序列，包括侧翼序列

**注意事项:**1：使用“seqkit grep”提取序列子集。“seqtk subseq seqs.fasta id.txt” 等于 “seqkit grep -f id.txt seqs.fasta”

Flags:
    --bed string      字符串由制表符分隔的bed文件
    --chr strings     使用——gtf或——bed(支持多值，大小写忽略)时，使用序列id选择序列
-d, --down-stream int 下游的长度
    --feature strings 选择有限的特性类型(支持多值，忽略大小写，只适用于GTF)
    --gtf string      by GTF (version 2.2) file
    --gtf-tag string  将这个标记作为序列注释输出(default "gene_id")
-h, --help            help for subseq
-f, --only-flank      只返回上行/下行流序列ce
-r, --region string   by region. e.g 1:12 for first 12 bases, -12:-1 for last 12 bases, 13:-1 for cutting first 12 bases. type "seqkit subseq -h" for more examples
-u, --up-stream int     上游的长度

stat

Flags:
-a, --all                 所有的统计信息, quartiles of seq length, sum_gap, N50
-b, --basename            仅输出文件的basename
-E, --fq-encoding string  fq-encoding字符串fastq质量编码。可用值:'sanger'， 'solexa'， 'illumina-1.3+'， 'illumina-1.5+'， 'illumina-1.8+' (default "sanger")
-G, --gap-letters string   gap letters (default "- .")
-h, --help                 help for stats
-e, --skip-err             skip-err skip错误，只显示警告信息
-i, --stdin-label string   label for replacing default "-" for stdin (default "-")
-T, --tabular             以机器友好的表格格式输出

示例：可以计算N 50，还可以显示序列的1,2,3分节

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ILAAdNvQ-1649156114936)(C:\Users\Administrator\AppData\Roaming\Typora\typora-user-images\image-20220405173213908.png)]

faidx

创建FASTA索引文件并提取子序列（主要目的是提取子序列的名称）

这个命令类似于“samtools faidx”，但有一些额外的要求:

输出带有-f标志的完整标题行
支持正则表达式为序列ID，标志为-r
如果你有大量的id，你可以使用: seqkit faidx seqs.fasta -l IDs.txt

Flags: 
-f, --full-head  full-head打印整个标题行，而不仅仅是ID。新的fasta索引文件以.seqkit结尾。将创建Fai
-h, --help       help for faidx 
-i, --ignore-case 忽略大小写     
-I, --immediate-output    立即打印输出，不使用写缓冲区
-l, --region-file string  字符串包含区域列表 
-r, --use-regexp          id为正则表达式。但是这里不支持subseq区域。

示例

seqkit faidx hairpin1.fa 	    # 提取子序列的基因信息 
smc-mir-12461	87	6132246	60	61
hsa-mir-9902-2	93	6132395	60	61
gga-mir-1784b	58	6132549	58	59

seqkit faidx hairpin1.fa -f 	# 子序列信息加上标题

convert

转换FASTQ质量编码之间的 Sanger, Solexa and Illumina

Flags: 
-d, --dry-run            dry run 
-f, --force              从 Illumina-1.8+ 转化到 Sanger, truncate scores > 40 to 40 
    --from string        source quality encoding. if not given, we'll guess it 
-h, --help               help for convert 
-n, --nrecords int       猜测编码质量的记录数(默认1000)
-N, --thresh-B-in-n-most-common int  猜Illumina 1.5的最常见品质前N的“B”阈值。(默认2)
-F, --thresh-illumina1.5-frac float  前导N记录中Illumina 1.5的派系阈值(默认为0.1)
   --to string            目标质量编码 (default "Sanger")

seqkit convert tests/Illimina1.8.fq.gz | seqkit head -n 1  # 默认转换fq文件质量值到1.8+
seqkit convert tests/Illimina1.8.fq.gz --to Illumina-1.5+ | seqkit head -n 1

grep

通过ID/名称/序列/序列特征搜索序列，不匹配允许

注意：

默认情况下，我们使用模式匹配序列ID，使用“-n/——By -name”来匹配全名，而不是只匹配ID。
与POSIX/GNU grep不同，默认情况下，我们将模式与整个目标(ID/full header)进行比较。请打开“-r/——use-regexp”以进行部分匹配。
当按序列搜索时，它是部分匹配的，并且正链和负链都要搜索。使用标志“-m/——max-mismatch”允许不匹配，可以增加“-j/–threads”的值来加速处理。
像“RYMM…”这样的退化碱基/残基也由标志-d支持。但是不要在正则表达式中使用退化基/残基，需要将它们转换为正则表达式，例如将“N”或“X”改为“.”。
当通过标志’-p’提供搜索模式(主题)时，请对包含逗号的模式使用双引号，如 -p ‘“A{2,}”’ or -p ““A{2,}”” 。因为命令行参数解析器接受逗号分隔值(CSV)用于多个值(主题)。文件中的模式不遵循此规则。
结果中的序列顺序与原始文件中的顺序一致，而不是查询模式的顺序。但是对于FASTA文件，你可以使用: seqkit faidx seqs. fasta --infile-list IDs.txt
对于多个模式，您可以多次设置“-p”，即-p pattern1 -p pattern2，或者通过“-f/——pattern-file”给出一个模式文件。

Flags:
-n, --by-name           匹配全名而不仅仅是ID
-s, --by-seq            在seq上搜索subseq，搜索正链和负链，使用标志-m/——max-mismatch允许不匹配
-c, --circular          环形基因组
-C, --count             只输出匹配记录的计数。使用-v/——reverse -match标志，计数不匹配的记录
-d, --degenerate        pattern/motif contains degenerate base
    --delete-matched    在匹配后立即删除一个模式，这将保留第一个匹配的数据，并在使用正则表达式时加速
-h, --help              help for grep
-i, --ignore-case       忽略大小写
-I, --immediate-output  立即打印输出，不使用写缓冲区
-v, --invert-match      反转匹配意义，选择不匹配的记录
-m, --max-mismatch int  通过seq匹配时最大的不匹配。对于像人类基因组这样的大型基因组，使用绘图/对齐工具将会更快
-P, --only-positive-strand   只在正链上搜索
-p, --pattern strings        字符串搜索模式(支持多个值) 
-f, --pattern-file string    字符串模式文件(每行一个记录)
-R, --region string 指定搜索序列区域（1:12 for first 12 bases, -12:-1 for last 12 bases）
-r, --use-regexp              正则表达式

示例

zcat hairpin.fa.gz | seqkit grep -r -p ^hsa-mir-807 # 正则匹配序列名称
zcat hairpin.fa.gz | seqkit grep -r -p ^hsa -p ^mmu -v  #2个条件并取反
zcat hairpin.fa.gz | seqkit grep -f list > new.fa  #将需要提取的序列名放在list中
zcat hairpin.fa.gz | seqkit grep -s -r -i -p ^aggcg  #正则匹配序列碱基，-i 忽略大小写
seqkit grep -s -R 1:30 -i -r -p GCTGG  #-R 在前30个碱基中正则匹配
seqkit grep -s -r -i -p ^atg cds.fa ＃选取有起始密码子的序列
seqkit grep -f list test.fa > new.fa ＃根据ID提取序列
seqkit grep -s -d -i -p TTSAA ＃简并碱基使用。S 代表C or G.
seqkit grep -s -R 1:30 -i -r -p GCTGG ＃匹配限定到某区域

rmdup

去除重复序列

Flags:
-n, --by-name               使用全名而不是id
-s, --by-seq                序列
-D, --dup-num-file string   字符串文件用于保存重复序列的数量和列表
-d, --dup-seqs-file string  字符串文件保存重复的seqs
-h, --help                  help for rmdup
-i, --ignore-case           忽略大小写
-P, --only-positive-strand  在按顺序比较时只考虑正链

seqkit rmdup hairpin1.fa -s -o clean.fa.gz # 去除重复序列
[INFO] 3967 duplicated records removed
seqkit rmdup hairpin1.fa -s -i -o clean.fa.gz -d duplicated.fa.gz #-d输出重复序列

common

通过id/name/sequence查找多个文件的公共序列

Flags:
-n, --by-name               使用全名而不是id
-s, --by-seq                序列
-h, --help                  help for common
-i, --ignore-case           忽略大小写
-P, --only-positive-strand  按序列比较时只考虑正链

示例

# 寻找两个文件的共有序列
seqkit common test.fa test1.fa -o common.fasta
# 根据序列名称
seqkit common test.fa test1.fa -n -o common.fasta
# 根据序列
seqkit common test.fa test1.fa -s -o common.fasta

split

分割序列文件的名称ID，子序列的给定区域，部分大小或数量的部分。如果您只是想按零件或尺寸分割，请使用“seqkit split2”，它也适用于成对和单端FASTQ。区域的定义是基于1的，并带有一些自定义设计。

Flags:
-i, --by-id              根据序列ID分割序列
-p, --by-part int        将序列拆分为N个部分
-r, --by-region string   字符串根据给定区域的子序列分割序列。如：1:12 for first 12 bases
-s, --by-size int        将切割成每个含N个序列
-d, --dry-run            dry-run，只打印信息，不会创建任何文件。
-f, --force              强制覆盖输出目录
-h, --help               help for split
-k, --keep-temp          在2通道模式下保留临时FASTA和.fai文件
-O, --out-dir string     out-dir字符串输出目录(default value is $infile.split)
-2, --two-pass           两遍模式读取文件两次，以降低内存使用。 (only for FASTA format)

示例

seqkit split hairpin.fa.gz -s 10000  #按序列数分割文件
[INFO] split into 10000 seqs per file
[INFO] write 10000 sequences to file: hairpin.fa.gz.split/hairpin.part_001.fa.gz
[INFO] write 10000 sequences to file: hairpin.fa.gz.split/hairpin.part_002.fa.gz
[INFO] write 10000 sequences to file: hairpin.fa.gz.split/hairpin.part_003.fa.gz
[INFO] write 8589 sequences to file: hairpin.fa.gz.split/hairpin.part_004.fa.gz

seqkit split hairpin.fa.gz -p 5 # 按文件个数分文件
INFO] split into 5 parts
[INFO] read sequences ...
[INFO] read 38589 sequences
[INFO] write 7718 sequences to file: hairpin.fa.gz.split/hairpin.part_001.fa.gz
[INFO] write 7718 sequences to file: hairpin.fa.gz.split/hairpin.part_002.fa.gz
[INFO] write 7718 sequences to file: hairpin.fa.gz.split/hairpin.part_003.fa.gz
[INFO] write 7718 sequences to file: hairpin.fa.gz.split/hairpin.part_004.fa.gz
[INFO] write 7717 sequences to file: hairpin.fa.gz.split/hairpin.part_005.fa.gz

seqkit split hairpin.fa.gz -i --id-regexp "^([\w]+)\-" -2  # 按ID
seqkit split hairpin.fa.gz -r 1:3 -2 #按前三个碱基

range

打印FASTA/Q记录在一个范围内(开始:结束)

Flags: 
-h, --help      help for range 
-r, --range string  范围：1:12 for first 12 records (head -n 12)

seqkit range  hairpin.fa.gz -r 101:150  # 打印指定范围的字符串

sort

按照id/名称/序列/长度对序列进行排序。默认情况下，所有记录都将被读入内存。对于FASTA格式，使用标志-2(——two-pass)来减少内存的使用。FASTQ不受支持的。首先，seqkit读取序列头部和长度信息。如果该文件不是普通的FASTA文件，seqkit将把序列写入临时文件，并创建FASTA索引。其次，seqkit根据头部和长度信息对序列进行排序，并通过FASTA索引提取序列。

Flags:
-b, --by-bases                由非间隙碱基组成
-l, --by-length               序列长度
-n, --by-name                 使用全名而不是id
-s, --by-seq                  序列
-G, --gap-letters string      gap letters (default "- \t.")
-h, --help                    help for sort
-i, --ignore-case             忽略大小写
-k, --keep-temp               在2通道模式下保留临时FASTA和.fai文件
-N, --natural-order           自然顺序排序，当按id /全名排序时
-r, --reverse                 反转结果
-L, --seq-prefix-length int   按序列排序的序列前缀长度(0表示整个序列)(默认10000)l
-2, --two-pass                两遍模式读取文件两次，以降低内存使用。(only for FASTA format)

seqkit sort hairpin1.fa -n # 序列全名称排序

2048 AI社区

有“AI”的1024 = 2048，欢迎大家加入2048 AI社区

更多推荐

基于Spring AI实现多轮对话系统架构设计

2048 AI社区

AI原生应用领域提示工程：构建智能服务体系的关键

2048 AI社区

特斯拉AI布局：从自动驾驶到人形机器人的野心

2048 AI社区

所有评论(0)

查看更多评论

承豆

@yearstime

已为社区贡献3条内容