seqkit 对fasta文件的各种操作

根据ID提取序列

#不包含空格
seqkit grep -f id.txt input.fa > output.fa
#可能包含空格
seqkit grep -n -f id.txt input.fa > output.fa

计算fasta序列长度

seqkit fx2tab --length --name --header-line test.fa > output.fa

反向互补fasta序列

#反向
seqkit seq test.fa -r > test_re.fa
#互补
seq test.fa -p > test_com.fa
#反向互补
seqkit seq test.fa -r -p > test_re_com.fa

DNA与RNA互转

#DNA转为RNA
seqkit seq test.fa --nda2rna > test.rna.fa
#RNA转为DNA
seqkit seq test.fa --rna2dna > test.dna.fa

序列以小写或者大写字母输出

#小写
seqkit seq test.fa -l > test_lower.fa
#大写
seqkit seq test.fa -u > test_upper.fa

过滤fasta序列长度

#长度大于60
seqkit seq -m 60 test.fa
#长度小于600
seqkit seq -M 600 test.fa
#长度大于60且小于600
seqkit seq -m 60 -M 600 test.fa

提取fasta的ID

#打印名称全行
seqkit seq test.fa -n > test.id
#仅仅打印ID
seqkit seq test.fa -n -i > test.id

删除序列中的gap（-）

seqkit seq -g test.fa > out.fa

截取指定区域的序列

#提取1-100个碱基
seqkit subseq -r 1:100 test.fa > out.fa
#提取后100个碱基
seqkit subseq -r -100:-1 test.fa > out.fa

根据gtf提取目标染色体上目标区域的序列

#提取1号染色体上cds序列
seqkit subseq --gtf test.gtf --chr 1 --feature cds  test.fa  > test_chr1_cds.fa
#提取序列多加上游的3个碱基
seqkit subseq --gtf test.gtf test.fa -u 3  > out.fa
#仅提取上游序列，如提取启动子区2k：-f仅定位不输出位置序列，-u输出上游序列3bp
seqkit subseq --gtf  test.gtf test.fa -u 3 -f > out.fa

滑窗提取序列

#步长为3，序列长度6
seqkit sliding -s 3 -W 6 test.fa > out.fa
#步长为3，序列长度6，最后一个不足6也取
seqkit sliding -s 3 -W 6 -g test.fa > out.fa
#环状DNA
seqkit sliding -s 3 -W 6 -C test.fa > out.fa

fastq转为fasta

seqkit fq2fa test.fq.gz -o test.fa

翻译DNA/RNA为蛋白质序列

#转为蛋白序列
seqkit translate test.fa > out.fa
# 去除'X' 和 '*'
seqkit translate test.fa
seqkit translate test.fa --trim > out.fa

提取包含指定碱基的fa

#以AGC开头的序列
seqkit grep -r -p ^AGC test.fa > out.fa
#移除包含AGC mmu的序列
seqkit grep -r -p ^AGC -p ^mmu -v test.fa > out.fa
#包含aggcg的序列
seqkit grep -s -i -p aggcg test.fa > out.fa
#去除 包含aggcg的序列
seqkit grep -s -r -i -p ^aggcg test.fa > out.fa

定位指定序列位置

seqkit locate -i -d -p AUGGACUN --bed
seqkit locate -i -d -p AUGGACUN --gtf

引物检索

#acgcccactgaaatga为扩增产物，指定正向引物序列为“ccc”，指定反向引物序列为“ttt”
echo -ne ">seq\nacgcccactgaaatga\n"  seqkit amplicon -F ccc -R ttt test.fa > out.fa
#输出格式为bed 
echo -ne ">seq\nacgcccactgaaatga\n"  seqkit amplicon -F ccc -R ttt --bed test.fa > out.bed

删除重复序列

seqkit rmdup -s test.fa -o out.fa

拆分fasta

#按照10000个序列拆分，-2 表示指定输出文件的前缀为 2
seqkit split2 test.fa -s 10000 -2
#按照前10个碱基拆分
seqkit split2 test.fa -r 1:3 -2
#拆双端文件，-p 拆分为两个文件，-O out 表示输出文件的目录。这里输出目录设置为 out，表示输出的文件将保存在名为 out 的文件夹中，-f 表示覆盖已存在的输出文件。
seqkit split2 -1 C1_1.fq.gz -2 C1_2.fq.gz -p 2 -O out -f

fastq随机抽样

#随机抽取10%的reads
seqkit sample -p 0.1  test.fq.gz -o test_10.fq
#双端需要设置种子数，保证抽出来的id一致
seqkit sample -s 123 -p 0.1  test.fq1.gz -o test_10_1.fq
seqkit sample -s 123 -p 0.1  test.fq2.gz -o test_10_2.fq
#抽取1000条reads
seqkit sample -n 1000 test.fq.gz -o test_sample.fq