Seqtk、Seqkit两个处理fa/fq神器的学习记录
Seqtk
安装 # Conda也可
git clone https://github.com/lh3/seqtk cd seqtk make
1.将fastq 文件转换成fasta 文件
seqtk seq -A input.fastq > output.fasta
2.得到反向互补序列
seqtk seq -Ar input.fastq > output.fasta
3.seqtk comp: 得到fastq/fasta 文件的碱基组成
（输出格式：序列id 序列长度 A C G T ）
seqtk comp in.fa > out.fa
4.subseq 根据name.list（不带>符号）提取子序列 -l可设定输出的每行长度
seqtk subseq -l 40 B1_NR100nl.fasta name.list > out.fa
5.随机抽取序列,按照比例或者数量 -s设定随机种子，便于重复
seqtk sample -s 10 test.fq 0.4 #比例 seqtk sample -s 10 test.fq 100 #数量
6.重命名 会将序列id变为从1到n....
seqtk rename in.fa <前缀> > out.fa
7..fastq转换为fasta，支持压缩格式
seqtk seq -a in.fq.gz > out.fa.gz
8.使用Phred算法从两端修剪低质量的碱基:
seqtk trimfq in.fq > out.fq
9.从每个读数的左端修剪5bp，从右端修剪10bp
seqtk trimfq -b 5 -e 10 in.fa > out.fa
10.合并两个序列,通常我们Illunima双端测序会得到两个文件R1.fq.gz和R2.fq.gz，这个命令就是帮助怎么完美实现两两配对
$seqtk mergefa Usage: seqtk mergefa [options] <in1.fa> <in2.fa># 合并两个的FASTA/Q files Options: -q INT quality threshold [0] -i take intersection#取交集 -m convert to lowercase when one of the input base is N -r pick a random allele from het -h suppress hets in the input
- cutN 去掉序列中的N/或者查看N所在的位置信息（-n 指定N的最小长度） -g 只打印出位置信息
Usage: seqtk cutN [options] <in.fa> Options: -n INT min size of N tract [1000] -p INT penalty for a non-N [10] -g print gaps only, no sequence
Seqkit
这个软件真的feel good，现有的轮子就不用再去造了，只记录一些常用的小命令，后续用到再进行更新，大家可以去原网页学习完整教程，点Here
1.seq子命令 读，换, 查看类型,统计基本信息
seqkit seq hairpin.fa.gz #查看 cat in.fa | seqkit stats #自动检测类型并给出统计信息 file format type num_seqs sum_len min_len avg_len max_len - FASTA RNA 1 11 11 11 11 seqkit stats *.f{a,q}.gz -T | csvtk pretty -t #可以统计多个fa/fq信息，结果形式更加友好 # -j 参数会并行运算更快速度 seqkit seq hairpin.fa.gz -n #打印序列ID全名 seqkit seq hairpin.fa.gz -n -i #只打印ID seqkit seq hairpin.fa.gz -n -i --id-regexp ... #通过正则去打印ID中内容(适用于所有的子命令) seqkit seq hairpin.fa.gz -s -w 0 #只打印序列(全局标志-w定义输出行宽，0表示不换行) seqkit seq reads_1.fq.gz -w 0 #转换多行fq为单行fq seqkit seq hairpin.fa.gz -r -p #反向互补 echo -e ">seq\nACGT-ACTGC-ACC" | seqkit seq -g -u #移除gap cat hairpin.fa | seqkit seq -m 100 -M 1000 | seqkit stats #过滤序列长度并进行统计
2. subseq子命令
前12碱基 $ zcat hairpin.fa.gz | seqkit subseq -r 1:12 后12碱基 zcat hairpin.fa.gz | seqkit subseq -r -12:-1 过滤前后12碱基 zcat hairpin.fa.gz | seqkit subseq -r 13:-13 通过gtf文件得到序列 seqkit subseq --gtf t.gtf t.fa 通过bed文件得到序列并移除重复序列 seqkit subseq --bed Homo_sapiens.GRCh38.84.bed.gz --chr 1 hsa.fa | seqkit rmdup > chr1.bed.rmdup.fa
3. sliding
sliding sequences, circular genome supported Usage: seqkit sliding [flags] Flags: -C, --circular-genome circular genome. -g, --greedy greedy mode, i.e., exporting last subsequences even shorter than windows size -s, --step int step size -W, --window int window size
4.faidx （类似samtools中的faidx）
Usage: seqkit faidx [flags] <fasta-file> [regions...] Flags: -f, --full-head print full header line instead of just ID. New fasta index file ending with .seqkit.fai will be created -h, --help help for faidx -i, --ignore-case ignore case -r, --use-regexp IDs are regular expression. But subseq region is not suppored here. seqkit faidx tests/hairpin.fa hsa-let-7a-1 hsa-let-7a-2 #提取指定序列 -f 输出ID全名 1:10 输出相应位置序列
5.fq转换fa
seqkit fq2fa reads_1.fq.gz -o reads_1.fa.gz
6.将FASTA/Q转换为表格格式，并提供各种信息，
如序列长度 GC
$ seqkit fx2tab hairpin.fa.gz -l -g -n -i -H | head -n 4 | csvtk -t -C '&' pretty #name seq qual length GC cel-let-7 99 43.43 cel-lin-4 94 54.26 cel-mir-1 96 40.62 #两种形式转换 zcat hairpin.fa.gz | seqkit fx2tab | seqkit tab2fx #按照长度排列序列 seqkit sort -l hairpin.fa.gz #得到前1000条reads seqkit fx2tab hairpin.fa.gz | head -n 1000 | seqkit tab2fx
7. grep序列
zcat hairpin.fa.gz | seqkit grep -r -p ^hsa #提取ID开头为hsa的reads -v取想反 zcat hairpin.fa.gz | seqkit grep -f list > new.fa #根据list取子集 cat hairpin.fa.gz | seqkit grep -s -i -p aggcg #提取序列里有AGGCG的reads -m 允许误配的数量 zcat hairpin.fa.gz | seqkit grep -s -r -i -p TT[CG]AA #带有模糊碱基的序列匹配 -R 1:30取前30 个碱基
8.duplicate
cat tests/hairpin.fa | seqkit head -n 1 \ | seqkit duplicate -n 2 #重复序列2次
9.rmdup
移除重复序列 by id/name/sequence Usage: seqkit rmdup [flags] Flags: -n, --by-name by full name instead of just id -s, --by-seq by seq -D, --dup-num-file string file to save number and list of duplicated seqs -d, --dup-seqs-file string file to save duplicated seqs -h, --help help for rmdup -i, --ignore-case ignore case
10. sample
zcat hairpin.fa.gz | seqkit sample -p 0.1 -o sample.fa.gz #按照比例取序列 zcat hairpin.fa.gz | seqkit sample -n 1000 -o sample.fa.gz #按照数量
11. rename
cat in.fa | less #和seqtk中rename的区别是前者会从1到n重新排序，后者是对后来重复的内容加_2到_n的后缀 >a comment acgt >b comment of b ACTG >a_2 a comment aaaa
OK,相信大家应该体会到这俩个软件的强大之处了，学好会省去你不少的time快去练习吧~
