(全文約1520字)
【推薦】用Smudgeplot評估物種倍性後，用組合jellyfish+GenomeScope1.0做二倍體物種的基因組調查，用組合KMC+GenomeScope2.0做多倍體物種的基因組調查。
1. k-mer進行基因組調查的軟件
k-mer進行基因組調查分為k-mer頻數統計和基因組特徵評估兩步。
- jellyfish可以實現第一步k-mer頻數統計。
- jellyfish的結果sample.histo可以用在GenomeScope上，實現第二步基因組特徵評估。
2. jellyfish 簡介
jellyfish是Center for Bioinformatics and Computational Biology在2011年研發的一款對DNA的k-mers計數的軟件，用Hash表儲存數據，能多線程運行。
3. jellyfish 安裝
- conda安裝
conda install -c bioconda jellyfish#安裝的是v2.2.10- github安裝
- 在github：jellyfish上通過源碼安裝。
4. jellyfish 運行
一般先用jellyfish count進行k-mer計數，然後用jellyfish histo對結果進行統計，獲得k-mer的頻數分佈直方表sample.histo。
4.1. count —— k-mer計數
- 命令
jellyfish count -m 17 -s 10G -t 12 -C -o sample.jf <(zcat sample_1.fq.gz) <(zcat sample_2.fq.gz)
- 參數
-
sample_1.clean.fq sample_2.clean.fq
使用的PE reads，不支持壓縮格式*.fq.gz輸入文件，如果不解壓縮，也可以用
<(zcat sample_1.fq.gz) <(zcat sample_2.fq.gz)代替sample_1.fq sample_2.fq; 或者使用這種形式zcat *fq.gz | jellyfish count /dev/fd/0，其中/dev/fd/0是進程輸入標誌，代表管道前結果傳遞。 -
-m 17
k-mer長度設置為17bp。如果基因組大小為G(單位是bp)，k-mer長度推薦設置成log(200*G)/log(4)。500Mbp的基因組對應約為17，1Gbp的19，10Gbp的21。
-
-s 1000M
存儲用的hash表大小為1000M，這個參數識別單位M(Mbp)和G(Gbp)。若該值不夠大，則會生成多個hash文件，以數字區分文件名。最好設置的值大於總的獨特的(distinct)k-mer數，這樣生成的文件只有一個。如果基因組大小為G，每個reads有一個錯誤，總共有n條reads，則該值可以設置為[(G + n)/0.8]。
-
-t 12
線程12
-
-C
對DNA正負鏈都進行統計，表示考慮DNA正義與反義鏈，遇到反義kmer時，計入正義kmer頻數中。如果是雙端測序reads，需要這個參數。
-
-o sample.jf
結果文件名為sample.jf，會生成k-mer計數文件sample.jf，是hash的二進制文件。
-
c 7
k-mer的計數結果所佔的最大比特數，默認支持的最大數字是2^7=128。該值最大，消耗內存越大。
-
-out-counter-len=4
輸出的二進制hash文件中的計數結果所佔的字節數，一個字節是8比特，則默認支持的最大數字是2^32=4.3G。
- 不推薦用-Q，會將低質量的鹼基替換成N。
-
-L
不輸出低於此值的k-mer
-
-U
不輸出高於此值的k-mer
- 輸出
-
sample.jf
hash格式儲存的k-mer頻數文件
4.2. histo —— 統計k-mer頻率
- 命令
jellyfish histo -t 12 sample.jf > sample.histo
統計k-mer計數(sample.jf)得到k-mer頻數分佈直方表(sample.histo)。
- 參數
-
-t 12
線程12。
-
-l 1
x的最小值，默認是1。結果會將小於此值的所有的k-mer的數目作為(x‐1)的值總結到一行。
-
-h 10000
x的最大值，默認是10000。結果會將大於此值的所有的k-mer的數目作為(x+1)的值總結到一行。
-
-i 1
x軸取值間隔，每隔該數值取值，默認為1。
- 結果
- k-mer頻數分佈直方表(sample.histo)包含空格分隔的兩列數據。
- 第一列代表k值出現的次數x(x=1,2,3...)，第二列是出現了x次的kmer的種類的數量y。
- sample.histo的兩列即是kmer分佈頻率直方圖的x和y軸的值。
4.3. merge 合併【按需選擇】
如果jellyfish count模塊輸出結果的二進制hash文件有多個，需要將多個hash文件合併，合併到merge.jf。
jellyfish merge sample_hash1.jf sample_hash2.jf sample_hash3.jf -o merge.jf
4.4. stats 統計【可選】
jellyfish stats sample.jf -o counts_stats.txt
可以用stats模塊來統計出k-mer總數（Total），特異的k-mer數目（Distinct），只出現過一次的k-mer數量（Unique），頻數最高的k-mer數量（Max_count）等信息。
5. 基因組特徵評估
獲得k-mer頻數分佈表sample.histo後，推薦用GenomeScope1.0或者GenomeScope2.0或者GenomeScope的R腳本來做基因組特徵評估和畫圖。也可直接用R繪製sample.histo的頻率分佈直方圖/頻率分佈曲線。
5.1. GenomeScope 網頁版
5.1.1. GenomeScope1.0 網頁版 —— 適用於二倍體物種
- 在GenomeScope1.0 網頁版上傳前一步獲得的k-mer頻數分佈表sample.histo文件。
- 設置參數k-mer length為第一步選擇的k-mer長度值，這裏是17；參數Read length為序列讀長，一般為150；最後一個參數Max kmer coverage建議修改成更大的10000，以統計更多的k-mers。
- 結果顯示預估的基因組大小，雜合度，重複率等信息。
5.1.2. GenomeScope2.0 網頁版 —— 適用於多倍體物種
GenomeScope2.0 網頁版也是類似的步驟。
5.2. R繪製
R繪製k-mer頻數分佈曲線初步查看基因組特徵。
獲得kmer_plot.png為頻數分佈曲線，可根據曲線峯值對基因組大小進行計算和預估。
#R 腳本示例
kmer <- read.table('sample.histo')
kmer <- subset(kmer, V1 >=5 & V1 <=500) #對頻數範圍5-500的數據進行繪製
Frequency <- kmer$V1
Number <- kmer$V2
png('kmer_plot.png')
plot(Frequency, Number, type = 'l', col = 'blue')
dev.off()6. references
- jellyfish paper：https://academic.oup.com/bioi...
- jellyfish github：https://github.com/gmarcais/J...
- jellyfish參數推薦：https://www.bilibili.com/read...
- chenlianfu blog: jellyfish參數推薦：http://www.chenlianfu.com/?p=806
