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

11. 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_{快去练习吧}~