SeqKit的學習 --20191017

軟件的介紹

SeqKit是一種跨平臺的闯参、極快的，全面的fasta/q處理工具级历。SeqKit為所有的主流操作系統(tǒng)提供了一種可執(zhí)行的雙元文件玩讳，包括Windows熏纯，Linux樟澜，Mac OS X秩贰，并且不依賴于任何的配置或預先配置就可以直接使用。

軟件的安裝

## Install via conda
conda install -c bioconda seqkit

軟件的命令

## 序列和子序列
**seq**  轉(zhuǎn)換序列（序列顛倒丙唧，序列互補想际，提取ID）
**subseq** 從區(qū)域/gtf/bed中獲得序列胡本，包括側(cè)面的序列
**sliding** 滑動序列，支持環(huán)式基因組
**stats**   對FASTA/Q files進行簡單統(tǒng)計
**faidx** 創(chuàng)造fasta索引文件并提取子序列
**watch** 檢測并連線序列特點的柱狀圖
**sana** 清除質(zhì)量不好的單線的fastq文件
## 格式轉(zhuǎn)換
**fx2tab**  將FASTA/Q 文件轉(zhuǎn)變成表格形式 (1th: name/ID, 2nd: sequence, 3rd: quality)
**tab2fx** 轉(zhuǎn)變表格形式為fasta/q格式
**fq2fa** 轉(zhuǎn)變fastq文件為fasta文件
**convert** 在Sanger, Solexa and Illumina中轉(zhuǎn)換fastq的質(zhì)量編碼
**translate** 將DNA/RNA序列轉(zhuǎn)變成蛋白序列（支持模棱兩可的堿基）
## 搜索
**grep** 根據(jù)ID/名稱/序列/序列motif 搜索序列，且允許錯配
**locate** 定位子序列/motif,且允許錯配
**fish** 使用本地比對在較大序列中尋找短序列
**amplicon** 經(jīng)由引物檢索擴增子（或它附近特定的區(qū)域）
## bam文件的處理和監(jiān)視
**bam** 監(jiān)視和連線bam文件記錄特點的直方圖
## 設置參數(shù)
**head** 打印第一個Nfasta/q的記錄
**range** 在一個范圍內(nèi)（start:end）打印fasta/q的記錄
**sample** 通過數(shù)量或比例來體驗序列
**rmdup** 通過id/名稱/序列 來去除復制的序列
**duplicate**  復制N次的序列
**common** 通過id/名稱/序列 發(fā)現(xiàn)多條序列中共有的序列
**split** 通過id/seq region/size/parts (mainly for FASTA) 將序列劈開成文件
**split2** 將序列通過大小或部分 劈開成文件
## 編輯
**replace** 通過規(guī)律表達來代替名字或序列
**rename** 重新命名復制的ID
**restart** 為環(huán)狀基因組重新設置起始位置
**concat** 從多個文件中經(jīng)由相同的ID來連接序列
**mutate** 編輯序列（點突妆档，插入，刪除）
## 排序
**shuffle** 變換序列位置
**sort** 將序列經(jīng)由id/name/sequence 進行排序

軟件命令詳解

Sequence ID
大部分的軟件须板，包括seqkit默認將主導的非空格字母作為ID习瑰。

舉例說明軟件如何使用

##下載參考序列甜奄，一個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 的使用（用來統(tǒng)計序列的信息）
Usage: seqkit fx2tab [flags]
Flags:
-a, --alphabet 打印字母表字母
-q, --avg-qual 打印read的平均質(zhì)量
-B, --base-content strings 輸出指定的堿基含量
-g, --gc 輸出GC含量
-H, --header-line 打印標題列
-l, --length 打印序列的長度
-n, --name 只打印名字，而沒有序列或者質(zhì)量
-i, --only-id 只打印基因的ID

從fastq/a中根據(jù)名字和ID提取序列子集

$ seqkit sample --proportion 0.001  duplicated-reads.fq.gz \
    | seqkit seq --name --only-id > id.txt  ##管道符前面的命令是隨機取總文件0.1%的序列唉擂，管道符后面的是提取前面的符合要求的序列的ID楔敌。
## 查看ID list內(nèi)容
$ 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 通過數(shù)量隨機提取序列（結(jié)果也許并不完全匹配）
-p, --proportion float 通過比例隨機提取序列
-s, --rand-seed int 隨機函數(shù)
-2, --two-pass 減少內(nèi)存占用
舉例：隨機抽取序列

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序列中找到合并堿基并找到它的位置(這個仿佛有點難度，不報錯也不打印內(nèi)容到屏幕)胜臊？黑忱？

$ 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抚吠，這個仿佛也有問題楷力，不報錯也不打印內(nèi)容到屏幕萧朝？检柬？厕吉？
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 保存數(shù)量的文件并列出重復的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
輸出的內(nèi)容
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

輸出結(jié)果：

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 來減少內(nèi)存的使用磁餐，不支持fastq格式的文件
Usage:seqkit sort [flags]
Flags: 
-l, --by-length  通過序列長度
-n, --by-name 通過全名而不是id
-s, --by-seq 通過序列
-i, --ignore-case 忽視大小寫
-r, --reverse 反向排序
-2, --two-pass 雙流程模式讀取文件來降低內(nèi)存使用

根據(jù)標題信息來拆分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轉(zhuǎn)換成新的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

## 根據(jù)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 的使用
對序列進行轉(zhuǎn)換（顛倒怔蚌，互補桦踊，提取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 只打印品質(zhì)？声离？术徊？
-s, --seq 只打印序列
--id-regexp string 對于解析ID的正則表達式赠涮，(default "^(\\S+)\\s?")

seqkit split 的使用
usage：seqkit split [flags]
-i, --by-id 根據(jù)序列ID切割序列
-p, --by-part int 將一份文件分割成N份
-s, --by-size int 將一個文件按照N條序列進行分割
-O, --out-dir string 輸出路徑
-2, --two-pass 降低內(nèi)存的使用

從一個text文件已知的字符串中搜索并替換fasta標題

$ seqkit grep --by-name --use-regexp --ignore-case --pattern hypothetical viral.1.protein.faa.gz \
    | seqkit head -n 3 | seqkit seq --name  ## 此命令也是有報錯笋除，說明沒有執(zhí)行好垃它？国拇？酱吝？

$ 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  輸出不匹配此模式的內(nèi)容
-m, --max-mismatch int 通過序列匹配時允許的最大錯配
-p, --pattern strings 匹配模式辰妙，支持連續(xù)寫多個模式密浑，匹配任一模式即輸出肴掷。
-f, --pattern-file string 支持匹配模式寫到一個文件中呆瞻，如要提取的序列ID痴脾。
-R, --region string 搜索特定區(qū)域序列赞赖，e.g 1:12 for first 12 bases, -12:-1 for last 12 bases
-r, --use-regexp 使用正則表達式辕近，必須加入此參數(shù),如^匹配首端移宅。同-p聯(lián)合使用漏峰。
舉例：
seqkit grep -s -r -i -p ^atg cds.fa＃選取正鏈中有起始密碼子的序列

seqkit grep -f ID.txt test.fa > new.fa＃根據(jù)ID提取序列

seqkit grep -s -d -i -p TTSAA＃簡并堿基使用浅乔。S 代表C or G.

seqkit grep -s -R 1:30 -i -r -p GCTGG＃＃匹配限定到某區(qū)域

usage:seqkit replace(通過正則表達式來代替名字或序列)
-s, --by-seq 代替seq
-i, --ignore-case 忽視大小寫
-K, --keep-key 
-p, --pattern string 搜索正則表達式
-r, --replacement string 替換物

從兩個配對端讀數(shù)的文件提取配對的reads

## 首先創(chuàng)造兩個不平衡的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