SeqKit的学习 --20191017

软件的介绍

SeqKit是一种跨平台的、极快的，全面的fasta/q处理工具。SeqKit为所有的主流操作系统提供了一种可执行的双元文件，包括Windows，Linux，Mac OS X，并且不依赖于任何的配置或预先配置就可以直接使用。

软件的安装

## Install via conda
conda install -c bioconda seqkit

软件的命令

## 序列和子序列
**seq**  转换序列（序列颠倒，序列互补，提取ID）
**subseq** 从区域/gtf/bed中获得序列，包括侧面的序列
**sliding** 滑动序列，支持环式基因组
**stats**   对FASTA/Q files进行简单统计
**faidx** 创造fasta索引文件并提取子序列
**watch** 检测并连线序列特点的柱状图
**sana** 清除质量不好的单线的fastq文件
## 格式转换
**fx2tab**  将FASTA/Q 文件转变成表格形式 (1th: name/ID, 2nd: sequence, 3rd: quality)
**tab2fx** 转变表格形式为fasta/q格式
**fq2fa** 转变fastq文件为fasta文件
**convert** 在Sanger, Solexa and Illumina中转换fastq的质量编码
**translate** 将DNA/RNA序列转变成蛋白序列（支持模棱两可的碱基）
## 搜索
**grep** 根据ID/名称/序列/序列motif 搜索序列，且允许错配
**locate** 定位子序列/motif,且允许错配
**fish** 使用本地比对在较大序列中寻找短序列
**amplicon** 经由引物检索扩增子（或它附近特定的区域）
## bam文件的处理和监视
**bam** 监视和连线bam文件记录特点的直方图
## 设置参数
**head** 打印第一个Nfasta/q的记录
**range** 在一个范围内（start:end）打印fasta/q的记录
**sample** 通过数量或比例来体验序列
**rmdup** 通过id/名称/序列 来去除复制的序列
**duplicate**  复制N次的序列
**common** 通过id/名称/序列 发现多条序列中共有的序列
**split** 通过id/seq region/size/parts (mainly for FASTA) 将序列劈开成文件
**split2** 将序列通过大小或部分 劈开成文件
## 编辑
**replace** 通过规律表达来代替名字或序列
**rename** 重新命名复制的ID
**restart** 为环状基因组重新设置起始位置
**concat** 从多个文件中经由相同的ID来连接序列
**mutate** 编辑序列（点突，插入，删除）
## 排序
**shuffle** 变换序列位置
**sort** 将序列经由id/name/sequence 进行排序

软件命令详解

Sequence ID
大部分的软件，包括seqkit默认将主导的非空格字母作为ID。

FASTA header	ID
>123456 gene name	123456
>longname	longname
>gi|110645304|ref|NC_002516.2| Pseudomona	gi|110645304|ref|NC_002516.2|

举例说明软件如何使用

##下载参考序列，一个fastq文件，两个fasta文件
wget http://data.biostarhandbook.com/reads/duplicated-reads.fq.gz
wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz
wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.protein.faa.gz

对fastq文件进行一个概括浏览

$ seqkit stat *.gz
file                      format  type     num_seqs    sum_len  min_len   avg_len  max_len
duplicated-reads.fq.gz    FASTQ   DNA        15,000  1,515,000      101       101      101
viral.1.1.genomic.fna.gz  FASTA   DNA             6    195,842    5,386  32,640.3  154,675
viral.1.protein.faa.gz    FASTA   Protein       112     55,855       38     498.7    3,122

在fasta/q文件中获取每条序列的GC含量

GC content
$ seqkit fx2tab --name --only-id --gc viral*.fna.gz
NC_001798.2         70.38
NC_030692.1         50.10
NC_027892.1         40.57
NC_029642.1         39.88
NC_001607.1         50.07
NC_001422.1         44.76

Custom bases (A, C and A+C)
$ seqkit fx2tab -H -n -i -B a -B c -B ac viral.1.1.genomic.fna.gz
#name   seq qual    a   c   ac
NC_001798.2         14.87   35.03   49.90
NC_030692.1         25.03   25.25   50.28
NC_027892.1         29.68   19.44   49.11
NC_029642.1         30.14   19.22   49.36
NC_001607.1         25.30   25.54   50.84
NC_001422.1         23.97   21.48   45.45

附seqkit fx2tab 的使用（用来统计序列的信息）
Usage: seqkit fx2tab [flags]
Flags:
-a, --alphabet 打印字母表字母
-q, --avg-qual 打印read的平均质量
-B, --base-content strings 输出指定的碱基含量
-g, --gc 输出GC含量
-H, --header-line 打印标题列
-l, --length 打印序列的长度
-n, --name 只打印名字，而没有序列或者质量
-i, --only-id 只打印基因的ID

从fastq/a中根据名字和ID提取序列子集

$ seqkit sample --proportion 0.001  duplicated-reads.fq.gz \
    | seqkit seq --name --only-id > id.txt  ##管道符前面的命令是随机取总文件0.1%的序列，管道符后面的是提取前面的符合要求的序列的ID。
## 查看ID list内容
$ head id.txt 
SRR1972739.2996
SRR1972739.3044
SRR1972739.3562
##通过ID list文件来搜索序列
$ seqkit grep --pattern-file id.txt duplicated-reads.fq.gz > duplicated-reads.subset.fq.gz

seqkit sample 的使用
-n, --number int 通过数量随机提取序列（结果也许并不完全匹配）
-p, --proportion float 通过比例随机提取序列
-s, --rand-seed int 随机函数
-2, --two-pass 减少内存占用
举例：随机抽取序列

seqkit sample -n 10000 -s 11 test1_1.fq -o sample.fq

seqkit sample -p 0.1 -s 11 test1_1.fq -o sample.fq

从fasta/q序列中找到合并碱基并找到它的位置(这个仿佛有点难度，不报错也不打印内容到屏幕)？？

$ seqkit fx2tab -n -i -a viral.1.1.genomic.fna.gz \
    | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]"
## 定位简并碱基，e.g, N and K，这个仿佛也有问题，不报错也不打印内容到屏幕？？？
seqkit grep --pattern-file id2.txt viral.1.1.genomic.fna.gz \
    | seqkit locate --ignore-case --only-positive-strand --pattern K+ --pattern N+

移去相同序列中重复的fasta/q记录

$ seqkit rmdup --by-seq --ignore-case duplicated-reads.fq.gz > duplicated-reads.uniq.fq.gz

[INFO] 7172 duplicated records removed

seqkit rmdup 的使用（用来通过id/名称/序列来去除重复的序列）
Usage: seqkit rmdup [flags]
Flags:
-n, --by-name 通过全名而不是id
-s, --by-seq 通过序列
-D, --dup-num-file string 保存数量的文件并列出重复的seqs
-d, --dup-seqs-file string 保存重复seqs的文件
-i, --ignore-case 忽视字母大小写

在fastq/a序列中定位motif/子序列/酶切位点

$ cat enzymes.fa
>EcoRI
GAATTC
>MmeI
TCCRAC
>SacI
GAGCTC
>XcmI
CCANNNNNNNNNTGG
$ seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz
输出的内容
seqID                 patternName   pattern           strand   start   end     matched
gi|526245010|ref|NC_021865.1|   MmeI   TCCRAC            +     1816    1821    TCCGAC
gi|526245010|ref|NC_021865.1|   SacI   GAGCTC            +     19506   19511   GAGCTC
gi|526245010|ref|NC_021865.1|   XcmI   CCANNNNNNNNNTGG   +     2221    2235 CCATATTTAGTGTGG
seqkit locate 的使用
Usage:
  seqkit locate [flags]
  -d, --degenerate 包含简并碱基模式和motif
  --gtf 输出为GTF格式
  -i, --ignore-case 忽视字母大小写
  -m, --max-mismatch int 通过序列匹配时允许的最大错配
  -G, --non-greedy 非贪婪模式，更快但是可能错过与其他重叠的motif
  -P, --only-positive-strand 只搜索正链
  -f, --pattern-file 模式或motif文件（fasta格式）
  -p, --pattern strings 搜索pattern或motif
 
seqkit locate -i -d -p AUGGACUN test.fa

输出结果：

seqID      patternName pattern strand start end matched

cel-mir-58a AUGGACUN  AUGGACUN  +     81 88   AUGGACUG

ath-MIR163 AUGGACUN  AUGGACUN  -     122 129  AUGGACUC

怎样通过长度来对大量的fasta文件进行排序

$ seqkit sort --by-length viral.1.1.genomic.fna.gz > viral.1.1.genomic.sorted.fa

seqkit sort 的使用
通过id/名称/序列/长度来排序，对于fasta格式的文件，使用flag -2 来减少内存的使用，不支持fastq格式的文件
Usage:seqkit sort [flags]
Flags: 
-l, --by-length  通过序列长度
-n, --by-name 通过全名而不是id
-s, --by-seq 通过序列
-i, --ignore-case 忽视大小写
-r, --reverse 反向排序
-2, --two-pass 双流程模式读取文件来降低内存使用

根据标题信息来拆分fasta序列

## 先对fasta文件的标题进行概括浏览
$ seqkit head -n 3 viral.1.protein.faa.gz | seqkit seq --name --only-id

YP_009137150.1
YP_009137151.1
YP_009137152.1

## 将默认的ID转换成新的ID
$ seqkit head -n 3 viral.1.protein.faa.gz \
    | seqkit seq --name --only-id --id-regexp "\[(.+)\]"  这一步并不懂？？？

Human alphaherpesvirus 2
Human alphaherpesvirus 2
Human alphaherpesvirus 2

## 根据header进行拆分,会生成一个文件夹
$ seqkit split --by-id --id-regexp "\[(.+)\]" viral.1.protein.faa.gz

   
seqkit head 的使用（首先打印第N条fasta/q记录）
Usage:seqkit head [flags]
Flags:
-n, --number int 先打印N个fasta/q的记录（默认是10）

seqkit seq 的使用
对序列进行转换（颠倒，互补，提取ID等）
Usage:seqkit seq [flags]
Flags:
-p, --complement 取互补序列
--dna2rna DNA到RNA
--rna2dna RNA到DNA
-G, --gap-letters string  gap letters (default "- \t.")
-l, --lower-case 用小写字母打印序列
-M, --max-len int 只打印短于最大长度的的序列
-n, --name 只打印name
-g, --remove-gaps 移去组装序列中的gap
-r, --reverse 取反向序列
-i, --only-id 只打印ID而不是全名
-q, --qual 只打印品质？？？
-s, --seq 只打印序列
--id-regexp string 对于解析ID的正则表达式，(default "^(\\S+)\\s?")

seqkit split 的使用
usage：seqkit split [flags]
-i, --by-id 根据序列ID切割序列
-p, --by-part int 将一份文件分割成N份
-s, --by-size int 将一个文件按照N条序列进行分割
-O, --out-dir string 输出路径
-2, --two-pass 降低内存的使用

从一个text文件已知的字符串中搜索并替换fasta标题

$ seqkit grep --by-name --use-regexp --ignore-case --pattern hypothetical viral.1.protein.faa.gz \
    | seqkit head -n 3 | seqkit seq --name  ## 此命令也是有报错，说明没有执行好？？？

$ seqkit replace --kv-file changes.tsv --pattern "^([^ ]+ )(\w+) " \
    --replacement "\${1}{kv} " --key-capt-idx 2 --keep-key viral.1.protein.faa.gz > renamed.fa


seqkit grep 的使用（通过ID/名称/序列/序列motif来搜索序列，允许错配）
Usage:seqkit grep [flags]
Examples:
 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n
-n, --by-name 通过全名来匹配而不是ID
-s, --by-seq 搜索seq中的子seq，只搜索正链，通过-m/--max-mismatch 来允许错配
-d, --degenerate 包含简并碱基的pattern和motif（简并碱基：一个符号代替某两个或者更多碱基，如编译丙氨酸的可以有4个密码子：GCU\GCC\GCA\GCG,这时生物学上为了方便，用字母N指代UCAG四个碱基，故说编译丙氨酸的密码子是GCN，其中N就是简并碱基。）
-i, --ignore-case 忽视大小写
-v, --invert-match  输出不匹配此模式的内容
-m, --max-mismatch int 通过序列匹配时允许的最大错配
-p, --pattern strings 匹配模式，支持连续写多个模式，匹配任一模式即输出。
-f, --pattern-file string 支持匹配模式写到一个文件中，如要提取的序列ID。
-R, --region string 搜索特定区域序列，e.g 1:12 for first 12 bases, -12:-1 for last 12 bases
-r, --use-regexp 使用正则表达式，必须加入此参数,如^匹配首端。同-p联合使用。
举例：
seqkit grep -s -r -i -p ^atg cds.fa＃选取正链中有起始密码子的序列

seqkit grep -f ID.txt test.fa > new.fa＃根据ID提取序列

seqkit grep -s -d -i -p TTSAA＃简并碱基使用。S 代表C or G.

seqkit grep -s -R 1:30 -i -r -p GCTGG＃＃匹配限定到某区域

usage:seqkit replace(通过正则表达式来代替名字或序列)
-s, --by-seq 代替seq
-i, --ignore-case 忽视大小写
-K, --keep-key 
-p, --pattern string 搜索正则表达式
-r, --replacement string 替换物

从两个配对端读数的文件提取配对的reads

## 首先创造两个不平衡的PE reads文件，整个过程不报错，但是不懂具体含义？？？
$ seqkit rmdup duplicated-reads.fq.gz \
    | seqkit replace --pattern " .+" --replacement " 1" \
    | seqkit sample --proportion 0.9 --rand-seed 1 --out-file read_1.fq.gz
$ seqkit rmdup duplicated-reads.fq.gz \
    | seqkit replace --pattern " .+" --replacement " 2" \
    | seqkit sample --proportion 0.9 --rand-seed 2 --out-file read_2.fq.gz
## overview
$ seqkit stat read_1.fq.gz read_2.fq.gz

file          format  type  num_seqs  sum_len  min_len  avg_len  max_len
read_1.fq.gz  FASTQ   DNA      9,033  912,333      101      101      101
read_2.fq.gz  FASTQ   DNA      8,965  905,465      101      101      101
## sequence headers
$ seqkit head -n 3 read_1.fq.gz | seqkit seq --name

SRR1972739.1 1
SRR1972739.3 1
SRR1972739.4 1

$ seqkit head -n 3 read_2.fq.gz | seqkit seq --name

SRR1972739.1 2
SRR1972739.2 2
SRR1972739.3 2

## 首先提取两个文件序列的ID，并计算它们的交集
$ seqkit seq --name --only-id read_1.fq.gz read_2.fq.gz \
    | sort | uniq -d > id.txt

$ # number of IDs
wc -l id.txt
8090 id.txt
## 然后用id.txt来提取reads
$ seqkit grep --pattern-file id.txt read_1.fq.gz -o read_1.f.fq.gz
$ seqkit grep --pattern-file id.txt read_2.fq.gz -o read_2.f.fq.gz
## 用md5sum来检查两个文件的IDs是否是相同的
$ seqkit seq --name --only-id read_1.f.fq.gz > read_1.f.fq.gz.id.txt
$ seqkit seq --name --only-id read_2.f.fq.gz > read_2.f.fq.gz.id.txt
$ md5sum read_*.f.fq.gz.id.txt
537c57cfdc3923bb94a3dc31a0c3b02a  read_1.f.fq.gz.id.txt
537c57cfdc3923bb94a3dc31a0c3b02a  read_2.f.fq.gz.id.txt

##sort 一下
$ gzip -d -c read_1.f.fq.gz \
    | seqkit fx2tab \
    | sort -k1,1 -T . \
    | seqkit tab2fx \
    | gzip -c > read_1.f.sorted.fq.gz
$ gzip -d -c read_2.f.fq.gz \
    | seqkit fx2tab \
    | sort -k1,1 -T . \
    | seqkit tab2fx \
    | gzip -c > read_2.f.sorted.fq.gz

将两个fasta文件连接成一个

$ cat 1.fa
>seq1
aaaaa
>seq2
ccccc
>seq3
ggggg

$ cat 2.fa
>seq3
TTTTT
>seq2
GGGGG
>seq1
CCCCC
一行命令
$ seqkit concat 1.fa 2.fa
>seq1
aaaaaCCCCC
>seq2
cccccGGGGG
>seq3
gggggTTTTT