作为生信入门训练,我们常用perl python等脚本语言实现对基因组文件的处理,练习常规的文本文件处理。最近再做单细胞方面的技术开发,需要自己手动合成一个嵌合两种物种的基因组,且要配套对应的gtf文件。用于评估单细胞技术的双胞率。
于是把seqkit这个由重庆第三军医大学Wei Shen 使用 go 语言开发用于处理 fa 和 fq 文件的利器拿出来温习一下,发现我的需求基本都可以满足。同等于处理bed文件的万能工具bedtools,感慨一套处理工具软件包可以取代一批生信工程师,这话确实如此。
一、序列操作:
1.取反向序列
seqkit seq test.fa -r > test_re.fa
2.取互补序列
seq test.fa -p > test_com.fa
3.取反向互补序列
seqkit seq test.fa -r -p > test_re_com.fa
4.DNA序列转换为RNA序列
seqkit seq test.fa --nda2rna > test_rna.fa
5.RNA序列转换为DNA序列
seqkit seq test.fa rna2dna > test_dna.fa
6.将序列以小写字母的形式输出
seqkit seq test.fa -l > test_lower.fa
7.将序列以大写字母的形式输出
seqkit seq test.fa -u > test_upper.fa
8.指定每行序列的输出长度(为0的话,代表为一整行,默认的输出 长度是60个碱基)
seqkit seq test.fa -w 50 > test_50.fa (指定序列的长度为50)
9.将多行序列转换为一行序列
seqkit seq test.fa -w 0 > test_w.fa
10.只输出序列
seqkit seq test.fa -s -w 0 > test_seq.fa
11.将只输出的序列的,指定每行输出的碱基数
seqkit seq test_seq.fa -s -w 40 > test_seq40.fa
12输出参数-o
seqkit seq test.fa -s -w 20 -o test_20.fa
二、Fasta/q之间以及与tab格式互换
1.将fataq文件转化为fasta格式.
seqkit seq fq2fa test.fq -o test.fa
2.将fasta格式转化为tab格式
seqkit fx2tab test.fa > test_tab.fa (没有seq参数)
三、序列信息统计
1.序列碱基含量
seqkit fx2tab -l -g -n -i -H test.fa (这些参数组合起来比较好看)
2.序列长度的整体分布统计
seqkit stat test.fa
seqkit grep [flags] #这部分内容用的比较广泛
参数:
-n, --by-name
匹配整个序列的名字,包含deion部分,而不是序列id。
-s, --by-seq
匹配序列
-d, --degenerate
pattern/motif 包含简并碱基
-i, --ignore-case
忽略大小写
-v, --invert-match
输出不匹配此模式的内容
-p,
匹配模式,支持连续写多个模式,匹配任一模式即输出。如-p ^ATG -p TAA$。注意该功能仅能正向匹配,不能实现对互补链匹配。
-f, --pattern-file string
支持匹配模式写到一个文件中,如要提取的序列ID。
-R, --region string
匹配位置选择。e.g 1:12 for first 12 bases, -12:-1 for last 12 bases
-r, --use-regexp
使用正则表达式,必须加入此参数,如^匹配首端。同-p联合使用。
举例:
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##匹配限定到某区域
五、motif定位
对grep的拓展,可以正反链同时匹配,输出匹配的位置。
seqkit locate [flags]
参数:
-d, --degenerate
pattern/motif contains degenerate base
-i, --ignore-case
ignore case
-P, --only-positive-strand
only search at positive strand
-p, --pattern value
search pattern/motif
-f, --pattern-file string
pattern/motif file (FASTA format)
举例:
seqkit locate -i -d -p AUGGACUN test.fa
输出结果:
seqID
patternName
pattern
strand
start
end
matched
cel-mir-58a
AUGGACUN
AUGGACUN
81
88
AUGGACUG
ath-MIR163
AUGGACUN
AUGGACUN
122
129
AUGGACUC
六、多个序列文件比较寻找相同的序列或者ID相同的序列
seqkit common [flags]
参数:
-n, --by-name
匹配整个序列的名字,包含deion部分,而不是序列id
-s, --by-seq
match by sequence
-i, --ignore-case
ignore case
-m, --md5
use MD5 reduce memory usage
举例:
1、By ID (default,>后面,空格之前的名字)输出ID名字相同的。
seqkit common test1.fa test2.fa -o common.fasta
2、By full name(整个序列的名字,包含deion部分)。输出序列名字相同的。
seqkit common test1.fa test2.fa -n -o common.fasta
3、输出要比较的文件中序列相同的序列
seqkit common test1.fa test2.fa -s -i -o common.fasta
4、输出要比较的文件中序列相同的序列 (for large sequences)
seqkit common test1.fa test2.fa -s -i -o common.fasta --md5
七、提取部分序列
如随机抽取10000条FASTQ序列做NT污染评估。同时他也可以对FASTA序列提取
seqkit sample [flags]
参数:
-n, --number int
sample by number (result may not exactly match)
-p, --proportion float
sample by proportion(按比例提)
-s, --rand-seed int
rand seed for shuffle (default 11)
-2, --two-pass
2-pass modelower memory
举例:随机抽取序列
seqkit sample -n 10000 -s 11 test1_1.fq -o sample.fq
seqkit sample -p 0.1 -s 11 test1_1.fq -o sample.fq
八、排序输出命令
seqkit sort [flags]
参数:
-l, --by-length
按照序列长度排序
-n, --by-name
by full name
-s, --by-seq
按照序列排序
-i, --ignore-case
按序列排序时忽略大小写
-r, --reverse
反向排序
-2, --two-pass
对于FASTA序列排序可以减少内存
举例:
seqkit sort -ltest.fa
九、文件切割
seqkit split [flags]
参数:
-i, --by-id
split squences according to sequence ID
-p, --by-part int
将一个文件分割成N 份
-s, --by-size int
将一个文件按照N 条序列一个文件进行分割
-O, --out-dir string
output directory (default value is infile.split)
-2, --two-pass
two-pass mode to lower memory usage(only FAST)
举例:
seqkit split hairpin.fa.gz -p 4
参考材料:
1.https://github.com/shenwei356/seqkit
2.https://www.cnblogs.com/huangyinger/p/10421805.html