Seqtk
安裝 # Conda也可
git clone https://github.com/lh3/seqtk
cd seqtk
make
1.將fastq 文件轉(zhuǎn)換成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 根據(jù)name.list（不帶>符號）提取子序列 -l可設(shè)定輸出的每行長度
seqtk subseq -l 40 B1_NR100nl.fasta name.list > out.fa
5.隨機抽取序列,按照比例或者數(shù)量 -s設(shè)定隨機種子，便于重復(fù)
seqtk sample -s 10 test.fq 0.4 #比例
seqtk sample -s 10 test.fq 100 #數(shù)量
6.重命名 會將序列id變?yōu)閺?到n....
seqtk rename in.fa <前綴> > out.fa
7..fastq轉(zhuǎn)換為fasta，支持壓縮格式
seqtk seq -a in.fq.gz > out.fa.gz
8.使用Phred算法從兩端修剪低質(zhì)量的堿基:
seqtk trimfq in.fq > out.fq
9.從每個讀數(shù)的左端修剪5bp，從右端修剪10bp
seqtk trimfq -b 5 -e 10 in.fa > out.fa
10.合并兩個序列,通常我們Illunima雙端測序會得到兩個文件R1.fq.gz和R2.fq.gz豺妓，這個命令就是幫助怎么完美實現(xiàn)兩兩配對
$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，現(xiàn)有的輪子就不用再去造了，只記錄一些常用的小命令琳拭，后續(xù)用到再進行更新训堆，大家可以去原網(wǎng)頁學習完整教程，點Here
1.seq子命令 讀白嘁，換, 查看類型,統(tǒng)計基本信息
seqkit seq hairpin.fa.gz #查看
cat in.fa | seqkit stats #自動檢測類型并給出統(tǒng)計信息
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 #可以統(tǒng)計多個fa/fq信息坑鱼，結(jié)果形式更加友好
# -j 參數(shù)會并行運算更快速度
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中內(nèi)容(適用于所有的子命令)
seqkit seq hairpin.fa.gz -s -w 0 #只打印序列(全局標志-w定義輸出行寬，0表示不換行)
seqkit seq reads_1.fq.gz -w 0 #轉(zhuǎn)換多行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 #過濾序列長度并進行統(tǒng)計
- 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文件得到序列并移除重復(fù)序列
seqkit subseq --bed Homo_sapiens.GRCh38.84.bed.gz --chr 1 hsa.fa | seqkit rmdup > chr1.bed.rmdup.fa
- 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 輸出相應(yīng)位置序列
5.fq轉(zhuǎn)換fa
seqkit fq2fa reads_1.fq.gz -o reads_1.fa.gz
6.將FASTA/Q轉(zhuǎn)換為表格格式权薯，并提供各種信息姑躲，
如序列長度 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
#兩種形式轉(zhuǎn)換
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
- grep序列
zcat hairpin.fa.gz | seqkit grep -r -p ^hsa #提取ID開頭為hsa的reads -v取想反
zcat hairpin.fa.gz | seqkit grep -f list > new.fa #根據(jù)list取子集
cat hairpin.fa.gz | seqkit grep -s -i -p aggcg #提取序列里有AGGCG的reads -m 允許誤配的數(shù)量
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 #重復(fù)序列2次
9.rmdup
移除重復(fù)序列 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
- 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 #按照數(shù)量
- rename
cat in.fa | less #和seqtk中rename的區(qū)別是前者會從1到n重新排序，后者是對后來重復(fù)的內(nèi)容加_2到_n的后綴
>a comment
acgt
>b comment of b
ACTG
>a_2 a comment
aaaa
OK,相信大家應(yīng)該體會到這倆個軟件的強大之處了盟蚣，學好會省去你不少的time快去練習吧~
