寫在前面
Emmm，前面推送了幾期好基友-比較基因組小王子的 WGDI 軟件使用和比較基因組圖表解讀。后臺常常受到消息，大體意思“快更新！”。原本，我和小王子商量的是一周一更。但這兩天小王子就要結(jié)婚了，接下一小段時間估計也忙事情。索性，今天我們把WGDI目前已成稿的內(nèi)容一次推出，以饗讀者。
開展批量化的 ks 計算
有比較基因組分析經(jīng)驗，更或者多少接觸過相關(guān)領(lǐng)域文章的朋友應(yīng)對ks有了解。ks值的大小可以反映一定時間尺度的演化時間。wgdi中ka、ks批量計算，通過paml中的yn00實現(xiàn)（相對于其他實現(xiàn)，往往更容易得到認(rèn)可）。
當(dāng)然，這需要安裝
pal2nal.pl和paml，以及需要 在wgdi安裝包目錄（site-package/wgdi）的conf.ini 中配置路徑。我們在前面的《安裝 | 比較基因組系列之一 - WGDI 軟件安裝與配置》一文中已詳細(xì)介紹安裝和配置方式。
一切準(zhǔn)備就緒。wgdi的操作仍然簡單。
首先是獲取參數(shù)
wgdi -ks ? >> total.conf
隨后編輯配置文件，調(diào)整參數(shù)即可
[ks]
cds_file = Grape.cds.fa
pep_file = Grape.pep.fa
align_software = muscle
pairs_file = collinearity file
ks_file = Grape.ks
其中，
- align_software 支持 muscle 和 mafft
- pairs_file 自動支持 colinearscan 和 mcscanx 的輸出結(jié)果，當(dāng)然也可以是 tab 分隔的 基因?qū)ξ谋疚募?/li>
- ks_file 即輸出文件
隨后，直接運行即可
wgdi -ks total.conf
計算需要一點時間，畢竟是大量基因?qū)?。大體結(jié)果如下。
注意到，
wgdi同時給出了ng86和yn00的計算結(jié)果，用戶可以自己選擇。如果基因?qū)λ悴怀鰜韐s值，那么默認(rèn)為 -1。使用
wgdi好處之一是支持?jǐn)帱c續(xù)跑。只要用戶保持輸出文件名字不變化，那么即使運行中斷，下一次運行，wgdi會從斷點開始計算，節(jié)省時間。針對來源不同的輸入基因?qū)ξ募鄳?yīng)計算結(jié)果都會輸出至用戶指定文件。
精細(xì)的共線性區(qū)塊分析 - BlockInfo
共線性和ks值計算結(jié)果信息量極大，事實上，并不容易展示。此外，嚴(yán)謹(jǐn)?shù)谋容^基因組分析項目上，我們常常需要對共線性結(jié)果進行篩選。因此，將相關(guān)信息匯總成表，會便于刪選。
wgdi中直接提供這一整合信息的功能。
首先獲取功能參數(shù)
wgdi -bi ? >> total.conf
調(diào)整對應(yīng)參數(shù)
vim total.conf
[blockinfo]
blast = Grape.blastp
gff1 = Grape.gff
gff2 = Grape.gff
lens1 = Grape.len
lens2 = Grape.len
collinearity = collinearity file
score = 100
evalue = 1e-5
repeat_number = 20
position = order
ks = Grape.ks
ks_col = ks_NG86
savefile = block_information.csv
運行即可
wgdi -bi total.conf
輸出結(jié)果解讀，直接在本地使用 Excel 打開即可
結(jié)果文件中內(nèi)容較多，但我們不需要被嚇到。因為更全面的信息，可以幫助我們做更準(zhǔn)確的判斷。以下，我們逐列說明具體含義：
-
id即共線性的結(jié)果的唯一標(biāo)識 -
chr1,start1,end1即參考基因組（點圖的左邊）的共線性范圍 -
chr2,start2,end2即參考基因組（點圖的上邊）的共線性范圍 -
pvalue即共線性結(jié)果評估，常常認(rèn)為小于0.01的更合理些 -
length即共線性片段長度 -
ks_median即共線性片段上所有基因?qū)?code>ks的中位數(shù)（主要用來評判ks分布的） -
ks_average即共線性片段上所有基因?qū)?code>ks的平均值 -
homo1,homo2,homo3,homo4,homo5與multiple參數(shù)有關(guān)，即一共有homo+multiple列
主要規(guī)則是：基因?qū)θ绻邳c圖中為紅色，賦值為1，藍色賦值為0，灰色賦值為-1（顏色參考前述wgdi點圖 相關(guān)推文）。共線性片段上所有基因?qū)x值后求平均值，這樣就賦予該共線性一個-1,1的值。如果共線性的點大部分為紅色，那么該值接近于1?？梢宰鳛楣簿€性片段的篩選。 -
block1,block2分別為共線性片段上基因order的位置。 -
ks共線性片段的上基因?qū)Φ?code>ks值 -
density1,density2共線性片段的基因分布密集程度。值越小表示稀疏
神之一手 - 使用 ks 對 dotplot 著色
純粹的點圖，我們只能模糊地看到片段是否復(fù)制更或者復(fù)制了多少次。如果將 ks 和 共線性 在點圖中展示，那么就可以清晰看到共線性上ks值的變化（Ks值，事實上又對應(yīng)了一定尺度上的演化時間）。
老規(guī)矩，第一步，獲取參數(shù)配置文件
wgdi -bk ? >> total.conf
修改配置文件即可
vim total.conf
[blockks]
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
blockinfo = block_information.csv
pvalue = 0.01
tandem = true
tandem_length = 20
markersize = 1
area = 0,2
block_length = 5
figsize = 8,8
savefig = Grape.ks.dotplot.pdf
開始繪制
wgdi -bk total.conf
由于
repeat_number 的影響，近期加倍的片段影響較小，如果對古老多倍化的某些共線性感興趣，強烈建議修改lens文件(去除一些染色體，保留感興趣的染色體)，針對性的計算dotplot和ks dotplot。從
ks dotplot可以清楚看到，共線性上ks值的變化。你會發(fā)現(xiàn)，近期加倍產(chǎn)生的共線性片段，ks的波動一般很?。ㄓ袝r候，觀測到波動，那么就說明是染色體結(jié)構(gòu)變異區(qū)）?；谶@個現(xiàn)象，我們在判定ks 峰值的時候，采用共線性的中位數(shù)要相比其他方式要準(zhǔn)確地多。另外，有幸發(fā)現(xiàn)神奇的現(xiàn)象，如水稻11,12號就屬于中獎了。
正確地計算 Ks 峰值
在判斷ks峰值的時候，目前仍有不少分析報道直接采用所有旁系同源基因ks值進行計算。但這存在不少問題。比如串聯(lián)重復(fù)基因積累，會影響Ks峰值分布?；谏厦娴?code>ks dotplot，我們可以直接選取近期加倍產(chǎn)生的共線性片段來計算ks峰值。
首先，獲取配置文件參數(shù)
wgdi -kp ? >> total.conf
編輯參數(shù)文件
vim total.conf
[kspeaks]
blockinfo = block_information.csv
pvalue = 0.05 # 可以用來簡單篩選
tandem_length = 20 # 過濾串聯(lián)重復(fù)
tandem = true # 同ks dotplot
block_length = 5 # 如調(diào)整為 10
ks_area = 0,10 # ks值的取值范圍
multiple = 1
homo = -1,1 # 在`blockinfo`介紹過，可以快速篩選不同的共線性片段，如 0-1 或 0.3-1
fontsize = 9
area = 0,3 # 橫軸的范圍
figsize = 10,6.18
savefig = Grape.ks_median.distri.pdf
savefile = Grape.ks_median.distri.csv # 共線性中位數(shù)的數(shù)據(jù)集合，用來繪制ks峰的分布圖，這個文件可以重新回到上一步 點圖 繪制，確定是否提取到感興趣的共線性區(qū)塊
得到這個圖之后，我們也可以通過點圖，重新確認(rèn)是否提取準(zhǔn)確。不斷修改參數(shù)，更或者手動篩選，直到準(zhǔn)確為止。
[blockks]
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
blockinfo = Grape.ks_median.distri.csv
pvalue = 0.01
tandem = true
tandem_length = 20
markersize = 1
area = 0,2
block_length = 5
figsize = 8,8
savefig = Grape.ks.dotplot.filter.pdf
舉個例子，簡單修改參數(shù)，就可以得到一個穩(wěn)定的分布，如
block_length設(shè)定為10，homo設(shè)定為0,1，area范圍是0,2。因此，葡萄近期加倍的峰值應(yīng)該為1.25。多次調(diào)整參數(shù)，可以測試峰值是否穩(wěn)定。
注意，Ks值峰值接近于0的時候，要小心，因為共線性片段中有很多等于0的值。這些等于0的值如何處理，如何判定分化時間，我個人沒確鑿的把握，建議重視。
寫在最后
比較基因組分析，事實上是一個非常精細(xì)的工作。軟件的使用固然很重要，但具體的分析思維更重要。wgdi是我個人基于近幾年在比較基因組工作的實戰(zhàn)經(jīng)驗開發(fā)出來的一款瑞士軍刀式軟件。幾乎覆蓋了目前絕大多數(shù)基礎(chǔ)的比較基因組分析所需功能。當(dāng)然，更多的功能，如祖先染色體重構(gòu)的，目前還在進一步測試，相信很快會與大家見面。
總的來說，希望通過這一系列教程，能讓大家更為全面地認(rèn)識wgdi，做出原則上無誤的比較基因組數(shù)據(jù)分析。
