? ? ? ? 我們對一個新物種進行基因組denovo測序時，一般需要先對該物種進行基因組survey，估計基因組大小、雜合率、重復(fù)序列比例、GC含量和倍性等信息，這些結(jié)果會直接影響我們后續(xù)測序的數(shù)據(jù)量及估計組裝難度。公司一般會先用二代測個100x的數(shù)據(jù)量來做這一步，現(xiàn)在二代已經(jīng)很便宜了，測50X和100X也差不了多少錢(當然也不是越多越好)，但在測序之前，建議通過核型分析和流式等方法預(yù)估一下基因組大小是多少再安排送測序。綜上，通過survey，我們能對所測物種的基因組的基本信息有一個初步了解，并得到高質(zhì)量的二代clean data 為后面的三代數(shù)據(jù)組裝后的scaffold進行糾錯。

? ? ? ? 在得到下機數(shù)據(jù)*.fq.gz后，先通過jellyfish計算k-mer分布頻率，然后通過Genomescope2估計基因組大小、重復(fù)和雜合度。首先，我們需要先知道什么叫K-mer

一、K-mer 分析原理：

1.??K-mer定義：

? ? ? ? 從一段連續(xù)序列中迭代地選取長度為K個堿基的序列，若每條序列的長度為L，則K-mer長度為K，那么可以得到L-K+1個K-mer，后面我取K=21來進行分析(針對長度為K的K-mer，對于A,T,C,G四種堿基類型，一共能產(chǎn)生的K-mer種類數(shù)為4的K次冪個，選擇K=21是為了產(chǎn)生足夠多的K-mer種類數(shù)去覆蓋整個基因組，當然，你也可以設(shè)置不同的K值進行探索，21是比較常用的)。

?????? 此外，針對一些物種可能會有較高的雜合率及多倍性，K-mer分布圖會有不同的特征，比如：雜合K-mer，重復(fù)K-mer，需要我們進一步了解掌握。

?

對于不同的基因組雜合度，Kmer分布如下：

簡單基因組：

圖中有一個明顯的主峰，沒有其它峰值，根據(jù)Kmer分布圖，可初步判斷該物種基因組為簡單基因組。

高重復(fù)基因組

圖中在49x深度和104x深度位置分別有1個峰，其中49x的位置為主峰即基因組正常期望深度，104x位置為重復(fù)峰位置。根據(jù)kmer 分布圖，可初步判斷該物種基因組為高重復(fù)基因組。

2. 基因組大小預(yù)估

? ? ? ?我們需要了解K-mer深度分布曲線含義：K-mer深度—K-mer深度頻率。從所有測序Reads中逐堿基取K-mer，并統(tǒng)計K-mer深度頻數(shù)分布，從而獲得K-mer深度分布曲線。根據(jù)曲線獲得K-mer深度估計值，用于估計基因組大小。

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??基因組大小=(基因切割的Kmer數(shù)目)/(主峰深度)

? ? ? ? K-mer分析適用于分析唯一主峰區(qū)域所占比例較大的基因組，當基因組雜合非常高或者重復(fù)序列比例非常大時，其影響可能導(dǎo)致無法通過K-mer分析正確估計基因組大小。將K-mer深度等于1的情況認為是錯誤情況，計算錯誤率，并用于修正基因組大小。

二、K-mer實戰(zhàn)

2.1 軟件安裝：

? ? ? ?Jellyfish是CBCB(Center for Bioinformatics and Computational Biology)的Guillaume Mar?ais 和Carl Kingsford 研發(fā)的一款計數(shù)DNA 的 k-mers 的軟件。該軟件運用 Hash 表來存儲數(shù)據(jù)，同時能多線程運行，速度快，內(nèi)存消耗小。該軟件只能運行在64位的Linux系統(tǒng)下。其文章于2011年發(fā)表在雜志Bioinformatics 上。在其cbcb官網(wǎng)上有pdf版本使用手冊，有詳細介紹。jellyfish的功能有：kmer計數(shù)；融合二進制的Hash結(jié)果；統(tǒng)計Hash結(jié)果；通過Hash結(jié)果來畫直方圖；將Hash結(jié)果輸出成文本格式；查詢指定k-mer的數(shù)目。

安裝jellyfish，按照官網(wǎng)介紹逐步安裝就好了，非常好安。

$ wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-2.3.0.tar.gz

$ tar zxvf jellyfish-2.3.0.tar.gz

$ mkdir jellyfish

$ cd jellyfish-2.3.0

$ ./configure --prefix=$PWD/../jellyfish

$ make -j 8

$ make install

也可以在github上下載，https://github.com/gmarcais/Jellyfish，里面有更多的詳細介紹。

2.2 Jellyfish count

? ? ? ?我們要用的是count 命令來進行K-mer計數(shù)，使用fastq文件在默認參數(shù)上和fasta文件沒有區(qū)別，生成的hash結(jié)果為*.jf為后綴的二進制文件。

? ? ? ?以count基本命令為例：

? ? ? ? jellyfish?count?-m?21?-s?100M?-t?10?-C?reads.fasta

? ? ? ? 參數(shù)含義可以通過jellyfish count --help來進行理解，(大部分參數(shù)默認就可以)：

? ? 其中我覺得可以根據(jù)自己需求進行設(shè)置的有：

? ? -m使用的k-mer的長度。如果基因組大小為G，則k-mer長度選擇為: k ~= log(200G)/log(4)；

? ? -s ?Hash 的大小，最好設(shè)置的值大于總的獨特的(distinct)k-mer數(shù),這樣生成的文件只有一個。若該值不夠大，則會生成多個hash文件，以數(shù)字區(qū)分文件名。如果基因組大小為G，每個reads有一個錯誤，總共有n條reads，則該值可以設(shè)置為[(G + n)/0.8],該值識別M 和 G；

? ? -t |--threads=<num>? default: 1? 使用的CPU線程數(shù)；

? ? -o |--output=<string>? default:mer_counts ?輸出的結(jié)果文件前綴；

? ? -c |--counter-len=<num>? default:7， k-mer的計數(shù)結(jié)果所占的比特數(shù),默認支持的最大數(shù)字是2^7=128。對于基因組測序覆蓋度為N，則要使設(shè)置的該值要大于N。該值越大，消耗內(nèi)存越大；

?????? 值得注意的是，我們的下機數(shù)據(jù)通常都是*fq.gz壓縮文件，這個時候通過zcat命令及管道符“|”進行計算，一定要保留/dev/fd/0 , 這個是進程輸入標志，也就是代表管道符前邊的結(jié)果傳遞，其它參照上述參數(shù)設(shè)置。

2.3. stats命令對hash結(jié)果進行統(tǒng)計

? ? ? ? k-mer的結(jié)果以hash的二進制文件結(jié)果給出，需要統(tǒng)計只出現(xiàn)過一次的kmer數(shù);特異的k-mer數(shù)目;總的k-mer數(shù);出現(xiàn)了最多的k-mer的數(shù)目.

$ jellyfish stats?*.jf? ?-o?kmer_count

?示例結(jié)果cat kmer_count：

? ? Unique:?? ?106239? ?#只出現(xiàn)過一次的k-mer的數(shù)目

? ? Distinct:? ??153852? #特異性的k-mer數(shù)目

? ? Total:???? ?2974220? #總的k-mer數(shù)目

? ? Max_count: 1085? ?#同一個k-mer出現(xiàn)最多的數(shù)目

2.4 histo輸出k-mer頻率直方圖數(shù)據(jù)，如：

? ? ? ? 對k-mer的計數(shù)結(jié)果有個直觀的認識，則需要統(tǒng)計出現(xiàn)了x(x=1,2,3…)次的kmer的數(shù)目y,以x，y為橫縱坐標畫出直方圖。使用 histo 命令能給出 x 和 y 對應(yīng)的值，將結(jié)果默認輸出到標準輸出。

????????得到*.histo后綴的文件后，就可以通過Genomescope2進行K-mer分布結(jié)果的可視化作圖了，并得到預(yù)估的基因組大小、雜合率、重復(fù)序列比例。

三、附：

? ? ? ? 普通基因組：單倍體、純合二倍體或者雜合度<0.5%，且重復(fù)序列含量<50%，GC含量為35%到65%之間的二倍體。

? ? ? ??復(fù)雜基因組：雜合度>0.5%，重復(fù)序列含量>50%，多倍體，GC含量處于異常的范圍（GC含量<35%或者GC含量>65%的二倍體）。

? ? ? ??二倍體復(fù)雜基因組進一步細分為:微雜合基因組（0.5%<雜合率<=0.8%)；高雜合基因組（雜合率>0.8%)；高重復(fù)基因組（重復(fù)序列比例>50%）。

參考資料：

https://github.com/gmarcais/Jellyfish/

http://www.chenlianfu.com/?p=806

https://www.jianshu.com/p/94da86093843

http://www.360doc.com/content/21/0714/12/76149697_986499807.shtml

?

?

本文使用 文章同步助手 同步