劉小澤寫于18.12.27
參考biostar handbook以及Wei Shen的SeqKit（處理fa/fq領(lǐng)域還比較出名）http://bioinf.shenwei.me/seqkit/，它可以處理壓縮的或者解壓的fa/fq
下載測試數(shù)據(jù)
wget http://data.biostarhandbook.com/reads/duplicated-reads.fq.gz
wget ftp://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.protein.faa.gz
小操作們
先看一下整體的數(shù)據(jù)統(tǒng)計(jì)
$ 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 7 197,450 1,608 28,207.1 154,675
viral.1.protein.faa.gz FASTA Protein 113 56,244 38 497.7 3,122
關(guān)于GC含量
使用fx2tab 將fq/fa的統(tǒng)計(jì)信息轉(zhuǎn)為制表符分割的表格，(fasta默認(rèn)包括ID、sequence；fastq默認(rèn)包括ID、sequence、quality）；另外可以自己指定length、GC等信息
# 默認(rèn)輸出：ID、sequence
$ seqkit fx2tab viral.1.1.genomic.fna.gz
# 指定輸出(輸出name、GC、length)
$ seqkit fx2tab --name --only-id --gc --length
viral.1.1.genomic.fna.gz
NC_024015.1 1608 44.22
NC_001798.2 154675 70.38
NC_030692.1 8908 50.10
NC_027892.1 8957 40.57
NC_029642.1 9006 39.88
NC_001607.1 8910 50.07
NC_001422.1 5386 44.76
# [--name可以用-n替代，--only-id可以用-i替代]
當(dāng)然，除了GC含量，任何堿基的含量也可以得到
$ seqkit fx2tab -H -n -i -B a -B c -B ac viral.1.1.genomic.fna.gz
#name seq qual a c ac
NC_024015.1 27.18 26.43 53.61
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
如何根據(jù)ID得到fa/fq的子集？
這個(gè)還真是一個(gè)高頻問題，許多時(shí)候就想要一個(gè)fastq/fasta中的某些ID的信息，利用seqkit grep可以方便解決
# 先做一個(gè)測試數(shù)據(jù)，包括一些ID信息(就是從duplicated-reads.fq.gz中隨機(jī)取樣)
$ seqkit sample --proportion 0.001 duplicated-reads.fq.gz | seqkit seq --name --only-id >id.txt
SRR1972739.2996
SRR1972739.3044
SRR1972739.3562
SRR1972739.4162
SRR1972739.5975
# 然后根據(jù)ID抽取其中的序列
$ seqkit grep --pattern-file id.txt duplicated-reads.fq.gz > dup_sub.fq.gz
# 用這個(gè)檢查一下是不是這些ID
$ seqkit seq --name --only-id dup_sub.fq.gz
找到質(zhì)量差的堿基（非ATCG）并定位它們
先用conda安裝好csvtk
$ seqkit fx2tab -n -i -a xxx.fq.gz | csvtk -H -t -f 4 -r -i grep -p "[^ATCG]"
# 因?yàn)檫@里的實(shí)例數(shù)據(jù)都是不錯(cuò)的，沒有低質(zhì)量(可能之前的序列結(jié)果又更新了)?？梢宰约赫乙粋€(gè)包括N堿基的序列試試，替換掉xxx.fq.gz就好
# csvtk參數(shù)表示：
# -H: --no-header-row
# -t: --tabs
# -f: --fileds
# -r: --use-regexp
# -i: --ignore-case
# grep -p: grep --pattern
假入找到了低質(zhì)量堿基，我們想過濾掉
# 先找到低質(zhì)量序列ID
$ seqkit fx2tab -n -i -a xxx.fq.gz | csvtk -H -t -f 4 -r -i grep -p "[^ATCG]" | csvtk -H -t cut -f 1 > id2.txt
# 然后找到這些ID并去除
$ seqkit grep --pattern-file id2.txt --invert-match xxx.fq.gz >clean.fa
或者我們想看看這些低質(zhì)量堿基所在序列的信息（包括基因ID、正負(fù)鏈、起始終止位點(diǎn)等）
$ seqkit grep --pattern-file id2.txt xxx.fq.gz | seqkit locate --ignore-case --only-positive-strand --pattern N+
# 其中--pattern可以多設(shè)置幾個(gè)，就找出包含多種堿基的序列，比如可以再加一個(gè) --pattern K+
# 這個(gè)+表示含有一個(gè)或者多個(gè)前面的堿基
想要移除重復(fù)的序列
$ seqkit rmdup --by-seq --ignore-case duplicated-reads.fq.gz > rmdup.fq.gz
# 如果序列非常大的話，可以用 -m/--md5，利用md5信息代替原始序列信息，來減少內(nèi)存占用
# 如果想根據(jù)序列名進(jìn)行去重，可以使用參數(shù) --by-name
序列定位motif/subsequence/enzyme digest sites
假入現(xiàn)在有一個(gè)文件存放的是motifs / enzyme digest sites
$ cat enzymes.fa
>EcoRI
GAATTC
>MmeI
TCCRAC
>SacI
GAGCTC
>XcmI
CCANNNNNNNNNTGG #表示我們這里就像匹配到CCA和TGG，中間保證9個(gè)堿基就好
我們?yōu)榱吮ＷC上面的XcmI這種情況也能成功進(jìn)行匹配，需要定義一個(gè)-d或者--degenerate參數(shù)
seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz
#或者
seqkit locate -d -i -f enzymes.fa viral.1.1.genomic.fna.gz
將fasta序列排序
可以按照序列長度
$ seqkit sort --by-length viral.1.1.genomic.fna.gz > sorted.fa
# 對于大型的fasta文件，可以使用-2或者--two-pass減少內(nèi)存使用
以上都是一些較為常用的小操作，對于更復(fù)雜的比如“拆分fasta”、“合并fasta“等，用到可以再去搜索
歡迎關(guān)注我們的公眾號～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com
