linux中fq格式转fa,fasta/fq文件处理万能工具——Seqkit学习记录
shenwei爪哥开发的处理Fasta/Fastq文件的万能工具。之前处理fq/fa文件时花时间写的一些脚本发现在seqkit里直接能一行命令就解决。实在是提升效率，整合流程中十分好的工具。本文是对Seqkit官方介绍(https://bioinf.shenwei.me/seqkit/usage/)的学习，参考学习的过程中可以对照着官方文档中的例子进行操作学习。
熟练的运用关键还是需要多练习，搭建分析流程时多多回顾记得使用此工具。
序列和子序列(Sequence and subsequence)
1. seq 序列操作
-p, -r 反向互补序列
-w 指定输出每列的序列长度。例如-w 0即为一行格式序列输出。
-n 输出序列的名字
-m 500 过滤序列长度小于500的序列
-n , -i 输出序列名字id
-i --id-regexp 根据正则匹配对序列名称进行操作(仅匹配正则匹配括号里的)。
>cel-mir-1 MI0000003 Caenorhabditis elegans miR-1 stem-loop
ATAAGCGCGCGCGCG
seqkit seq hairpin.fa.gz -i --id-regexp "^[^\s]+\s([^\s]+)\s"
2. subseq根据区域/gtf/bed文件提取序列，以1为开始。
-r 1:12 取每条序列的前12bp
--gtf a.gtf b.fas 根据gtf文件提取序列
--gtf --feature cds -u 1000提取cds序列以及上游1000bp启动子区序列。
--gtf a.gtf -u 1000 -f根据gtf文件仅选取上游1000bp的启动子序列。
seqkit subseq --gtf t.gtf t.fa -u 3 -f
3. sliding根据滑窗取序列
4. stats对序列fa/fq文件进行基本统计
-a *.fq.gz包括所有的统计信息
-T 输出\t分割的文件，可接下来进行管道操作
##接管道csvtk进行操作
seqkit stats *.f{a,q}.gz -T | csvtk pretty -t
## 转为markdown文件格式
seqkit stats *.f{a,q}.gz -T | csvtk csv2md -t
5. faidx创建类似于samtools faidx的index文件。
可用于提取某一序列的指定区域序列，并且可以根据正则匹配来匹配序列姓名
##提取某一序列20~30bp区域的序列。
seqkit faidx tests/hairpin.fa hsa-let-7a-1:20-30
格式的转换(Format conversion)
1. fq2fa如其名，fastq文件转换为fasta文件
2. fx2tab 将每条序列fa/fq转换为tab分割的一行格式，也可输出每条序列的一些基本信息。
-n 只输出name行(去序列及质量)，-i输出name行起使的id名
-l 输出序列的长度，-g 输出gc含量
-H 输出文件的header表头
管道接csvtk可以接 csvtk inter多个文件中的相同序列，csvtk join
## 根据序列长度排序，也可seqkit sort -l
$ zcat hairpin.fa.gz | seqkit fx2tab -l | sort -t"`echo -e '\t'`" -n -k4,4 | seqkit tab2fx
>cin-mir-4129 MI0015684 Ciona intestinalis miR-4129 stem-loop
UUCGUUAUUGGAAGACCUUAGUCCGUUAAUAAAGGCAUC
>mmu-mir-7228 MI0023723 Mus musculus miR-7228 stem-loop
UGGCGACCUGAACAGAUGUCGCAGUGUUCGGUCUCCAGU
3. convert FASTQ的质量值的转换。
4. translate 翻译DNA/RNA序列为amino acid序列。
--trim去除翻译终止密码子的*
-T 不同的国际翻译规则，有31种规则。
-f , --frame -3~3 不同的翻译框。
-l 1显示标准密码子翻译表。AAA->K... -L 1显示ambiguous模式，包含一些简并碱基的翻译。
序列的搜索和定位(Searching：grep and locate)
1. grep 根据id名(可正则匹配)or特定的序列模式(motifs)来 搜索/提取序列。
-n 根据序列全名来提取序列。
-p -r 根据一定的正则匹配pattern来提取序列。
-v反向提取序列，即去除匹配的序列。
-f 根据id list来批量提取序列，很好用。
-s -p atcgatcg 根据某段序列motif来提取序列。
-d 支持简并碱基模式。
-R 1:30 只对1~30位的序列进行匹配。
$ zcat hairpin.fa.gz | seqkit grep -r -p ^hsa
>hsa-let-7a-1 MI0000060 Homo sapiens let-7a-1 stem-loop
UGGGAUGAGGUAGUAGGUUGUAUAGUUUUAGGGUCACACCCACCACUGGGAGAUAACUAU
ACAAUCUACUGUCUUUCCUA
## 对hmmsearch结果进行整理，提取其中的id，并获得拟南芥中HDZIP3的cds序列，其中利用 -f -
grep -v "#" ath.out |cut -d '|' -f 1 |seqkit grep -r -f - Ath_cds.fas.gz
### 根据small.fq文件中的id名来掉出big.fq.gz文件中的序列
seqkit grep -f
##Extract sequences starting with AGGCG
seqkit grep -s -r -i -p ^aggcg hairpin.fa
##根据简并碱基来提取序列
seqkit grep -s -d -i -p TTSAA
seqkit grep -s -r -i -p TT[CG]AA
2. locate 定位并返回短序列(motifs)的具体位置(location)，支持简并碱基/正则匹配。
--bed输出bed6格式，--gtf输出gtf格式文件；
-p 具体搜寻的motif序列。
-d 支持简并碱基。-i -d -p AWDTTT
-P 仅搜寻正义链。
-m 1 允许mismatch。
-f 输入motif序列的fasta文件进行批量定位。
## 输出每条序列种含有的ORF序列 和 具体位置(正则匹配)
zcat hairpin.fa.gz | seqkit locate -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)"
### 定位motif的具体位置(简并碱基的使用)。ChIP-seq里会用到。
zcat hairpin.fa.gz | seqkit locate -i -d -p AUGGACUN
序列的集合运算操作(set operations)
1. duplicate 对序列复制重复-n次
## 取序列进行重复2次，并进行重新命名
cat hairpin.fa | seqkit head -n 1 | seqkit duplicate -n 2 | seqkit rename
2. rmdup一个文件内根据id/name/seq来去除重复序列。
-s 根据相同的序列来去除重复
-n 根据相同的姓名去除重复
-d file.fas 保存重复的序列。
-D file_detail.txt 重复序列的一些基本信息。
3. common多个文件中根据id/name/seq来找相同的序列。(如果只是两个文件，仅用seqkit grep可快速实现)
-n 根据相同的全名name，取共同的序列
-s 根据相同的序列来提取。
-i ignore case
--md5 处理大文件时根据md5值
#####输出ID名字相同的:
seqkit common test1.fa test2.fa -o common.fasta
#####输出序列名字相同的:
seqkit common test1.fa test2.fa -n -o common.fasta
#####输出要比较的文件中序列相同的序列:
seqkit common test1.fa test2.fa -s -i -o common.fasta
4. split 对fa文件根据name/id/subseq/区域进行分割/split2支持fq文件，更小内存，更快速度。
-s 10000 分割出每10000条seq一个子文件。
-p 4 分割出4个子文件。--two-pass减少内存
-f force，强制分割。
-i 根据id分割文件。
-r 1:3 根据seq的前3位序列进行分割。
-1 -2 split2对fastq文件分割操作
5. sample 随机取一定数量n或比例的序列。
-n 100 随机约位100个子序列。
-p 0.1 取大约为10%的子序列。
-s 设置随机种子。
6. head 取前-n个fq/fa子序列； range 取指定范围-r -10:-1的序列
Fastq/Fasta的编辑修改(Edit)
1. replace 可根据正则对name/seq进行修改。
-s 对序列进行更改。默认是对名字name进行修改
-p 要匹配的正则pattern "\s.*"。
-r 要替换成的pattern。 用单引号''。支持捕获式匹配capture matching '$1'
-k 哈希表类似的key-value对应表文件tsv文件。需要用到-r {kv} --keep key 如果没有对应的key-value对应，则保留原来的名字key不做修改。
###1. 去除id空格\s后的内容
$ echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p " .+"
>seq1
ACGT-ACGT
###2. 修改id中的-改为=
$ echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p "\-" -r '='
>seq1 abc=123
ACGT-ACGT
####3. 在每个字符后加空格，(capture matching)
echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p "(.)" -r '$1 ' -s
#### 4. 将每个序列名字改位seq_1,seq_2，利用{nr}
echo -e ">abc\nACTG\n>123\nATTT" | seqkit replace -p ".+" -r "seq_{nr}" --nr-width 5
#### 5. 根据已知的key-value的对应表文件来修改姓名。
seqkit replace -p ' (.+)$' -r ' {kv}' -k alias.txt test.fa --keep-key
2. rename 对姓名重复的id名进行修改。
3. restart start position for circular genome.
4. concat 根据相同的id名连接序列。
###连接序列的前两位bp和后两位bp
seqkit concat
5. mutate 对每条序列的特定位置进行突变，包括point mutation, insertion, deletion。
-p 1:X 对每条序列的第一位bp编辑改为X。(位置为1-base)
-d 1:2 对每一条序列的前两位进行删除。
-i 0:xx 对每条序列的第0位之后插入xx，即在每条序列之间插入xx
-s chr1,chr2仅对chr1,chr2两条序列进行编辑。
-r -s chr支持正则匹配所有chr的进行编辑
-s chr1,chr2 -v 反向匹配chr1,chr2
###对每条序列的多个位点进行点突变
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" | seqkit mutate -p 1:x -p -1:x --quiet
### 对每条序列的最后3位进行删除
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" | seqkit mutate -d -3:-1 --quiet
### 仅对chr1，chr2两条序列进行编辑
cat tests/hsa.fa | seqkit mutate -p -1:X -s chr1,chr2
排序(ordering)
1. shuffle 对fasta文件序列随机重排。
2. sort 根据id/name/sequence/length重新对序列进行排序。
-n 按照全名进行排序
-N, --natural-order 按照数字——字母的方式进行排序。
-i 忽略大小写。
-l 根据序列长度进行排序。
-s 根据序列进行排序。
-r 反向输出结果。
-2, --two-pass 减少内存使用。
--quiet 不输出中间信息。
### 根据数字——字母的方式进行排序
echo -e ">3\na\n>1\na\n>Y\na\n>x\na\n>Mt\na\n>11\na\n>2\na\n" | seqkit sort -N -i -2
### 根据序列id进行排序，忽略大小写
echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" | seqkit sort --quiet -s -i
### 根据序列长度进行排序
echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAAnnn\n>seq3\nacgt" | seqkit sort --quiet -l