在組裝未知基因組時(shí)，往往需要利用重測序數(shù)據(jù)提前進(jìn)行基因組調(diào)查，以獲取其基因組規(guī)模，雜合率，重復(fù)序列比例，GC含量等信息。從而更好地?cái)M定后繼測序策略。基因組調(diào)查可以采用kmers方法。kmers基因組調(diào)查分為kmers頻數(shù)統(tǒng)計(jì)和基因組評(píng)估兩步。原理已經(jīng)有大佬講得很清楚啦：http://www.itdecent.cn/p/94da86093843
這里以獼猴桃基因組hongyang為例，具體使用kmc+genomescope軟件進(jìn)行基因組調(diào)查。kmc會(huì)對(duì)重測序reads數(shù)據(jù)進(jìn)行kmers庫的構(gòu)建，genomescope則根據(jù)kmers庫進(jìn)行基因組特征評(píng)估。我們來嘗試復(fù)刻文章中的基因組調(diào)查結(jié)論。
1.下載基因組重測序數(shù)據(jù)
我們以獼猴桃hongyang為例，首先獲取獼猴桃基因組重測序數(shù)據(jù)。進(jìn)入NCBI網(wǎng)址即可查看具體信息：https://trace.ncbi.nlm.nih.gov/Traces/index.html?run=SRR9329821
找到對(duì)應(yīng)的重測序數(shù)據(jù)編號(hào)SRR9329821。
服務(wù)器安裝sra-tools后可以直接使用prefetch下載。
$ prefetch SRR9329821
2.sra文件轉(zhuǎn)化為fastq文件
下載得到的數(shù)據(jù)約為11.3G，我們接下來使用fastq-dump將sra文件處理為fastq文件。
$ fastq-dump --gzip --split-3 SRR9329821/
這里對(duì)于不知道是單端測序還是上端測序的sra文件，可一律采用--split-3，程序會(huì)自動(dòng)識(shí)別。這一步較為耗時(shí)，建議nohup掛載到后臺(tái)運(yùn)行。 運(yùn)行結(jié)束后產(chǎn)生SRR9329821_1.fastq.gz和SRR9329821_2.fastq.gz兩個(gè)雙端文件。
3.fq文件質(zhì)控
采用fastp軟件進(jìn)行文件質(zhì)控。
此處-l代表過濾長度在36bp以下的reads，-w設(shè)定滑窗規(guī)模為6，以此標(biāo)準(zhǔn)計(jì)算平均reads質(zhì)量，-q代表過濾平均質(zhì)量Q20以下的reads，--compression表示壓縮程度（1-9），越大代表壓縮過程速度最慢，越小代表速度快。
$ fastp -i SRR9329821_1.fastq.gz -I SRR9329821_2.fastq.gz -o SRR9329821_1.fastq.cleandata.gz -O SRR9329821_2.fastq.cleandata.gz -l 36 -q 20 -n 6 -w 6 --compression=6
至此，重測序數(shù)據(jù)準(zhǔn)備工作已經(jīng)完畢，下面開始制備kmer數(shù)據(jù)庫以及基因組調(diào)查。
需要用到軟件KMC，smudgeplot和genomescope，請(qǐng)?zhí)崆跋螺d并安裝。
三者都可以用conda下載：
$ conda install -c bioconda smudgeplot
$ conda install -c bioconda genomescope2
$ conda install -c bioconda kmc
安裝完畢后可執(zhí)行后繼步驟。
4.構(gòu)建重測序kmers庫
在剛剛質(zhì)控后的存有雙端測序*fq.gz文件的文件夾下，利用雙端測序數(shù)據(jù)構(gòu)建kmer庫。
$ mkdir tmp #建立臨時(shí)文件夾
$ ls SRR9329821_1.fastq.cleandata.gz SRR9329821_2.fastq.cleandata.gz > KiWi
$ kmc -k21 -t32 -m64 -ci1 -cs10000 @KiWi test_kiwi tmp #運(yùn)行kmc以構(gòu)建kmer庫，設(shè)定kmer構(gòu)建長度為21，線程占用數(shù)為32，內(nèi)存為64G。
$ kmc_dump -ci50 -cx3000 test_kiwi kmer21.dump #提取所有覆蓋度在50X-3000X的kmers
kmc執(zhí)行reads建庫命令，運(yùn)行完kmers建庫后會(huì)生成kmers數(shù)據(jù)庫，其數(shù)據(jù)分別保存于.kmc_pre和.kmc_suf兩個(gè)文件。kmc_dump命令會(huì)產(chǎn)生.dump文件，里面存放有kmers的list。
5.匹配雜合kmers對(duì)
為了后繼估計(jì)基因組雜合程度，需要提前統(tǒng)計(jì)雜合的kmers對(duì)。將上述結(jié)果.dump文件做為輸入，運(yùn)行后的結(jié)果保存于kmer_hongyang_coverages.tsv和kmer_hongyang_sequences.tsv兩文件中。-o相當(dāng)于指定兩個(gè)tsv結(jié)果文件的標(biāo)題。
$ smudgeplot.py hetkmers -o kmer_hongyang < kmer21.dump
至此為止所有數(shù)據(jù)準(zhǔn)備工作已經(jīng)結(jié)束啦，下面正式開始進(jìn)行基因組特征調(diào)查。
6.smudgeplot評(píng)估基因組倍型
smudgeplot軟件可以以kmer_hongyang_coverages.tsv為輸入可視化倍型，這里直接上代碼。
$ smudgeplot.py plot -t "hongyang" -q 0.99 kmer_hongyang_coverages.tsv
7.genomescope評(píng)估基因組特征
Genomescope是2017年發(fā)表于bioinformatic的一個(gè)工具，以處理一些高復(fù)雜度的基因組調(diào)查。第一個(gè)版本僅能預(yù)測二倍體基因組，第二個(gè)版本可以預(yù)測多倍體基因組特征。Genomescope利用kmer頻數(shù)統(tǒng)計(jì)結(jié)果，即KMC結(jié)果.hist文件進(jìn)行基因組評(píng)估，過程如下：
$ kmc_tools transform test_kiwi histogram kmer21_hongyang.hist -cx10000
$ genomescope2 -i kmer21_hongyang.hist -k 21 -p 2 -o . -n hongyang_genomescope
得到的峰圖可以查看kmers頻數(shù)分布圖，也可以得到summary.txt以查看基因組大小，雜合度，重復(fù)片段比例等詳細(xì)信息。
結(jié)果可見基因組預(yù)估大小為617M，雜合度為1.12%，重復(fù)率0.69。跟原文的結(jié)論還是非常接近的，只不過原文用的是19-mer進(jìn)行的基因組調(diào)查，本文用的21-mer，可能造成一些結(jié)果出入。另外，genomescope作圖的美觀性還是不錯(cuò)的，推薦大家試一試。
