寫在前面
Emmm誓竿，前面推送了幾期好基友-比較基因組小王子的 WGDI 軟件使用和比較基因組圖表解讀磅网。后臺(tái)常常受到消息，大體意思“快更新筷屡！”知市。原本，我和小王子商量的是一周一更速蕊。但這兩天小王子就要結(jié)婚了，接下一小段時(shí)間估計(jì)也忙事情娘赴。索性规哲，今天我們把WGDI目前已成稿的內(nèi)容一次推出，以饗讀者诽表。
開展批量化的 ks 計(jì)算
有比較基因組分析經(jīng)驗(yàn)唉锌，更或者多少接觸過相關(guān)領(lǐng)域文章的朋友應(yīng)對(duì)ks有了解。ks值的大小可以反映一定時(shí)間尺度的演化時(shí)間竿奏。wgdi中ka袄简、ks批量計(jì)算，通過paml中的yn00實(shí)現(xiàn)（相對(duì)于其他實(shí)現(xiàn)泛啸，往往更容易得到認(rèn)可）绿语。
當(dāng)然，這需要安裝
pal2nal.pl和paml，以及需要 在wgdi安裝包目錄（site-package/wgdi）的conf.ini 中配置路徑吕粹。我們?cè)谇懊娴摹栋惭b | 比較基因組系列之一 - 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 自動(dòng)支持 colinearscan 和 mcscanx 的輸出結(jié)果稳其，當(dāng)然也可以是 tab 分隔的 基因?qū)ξ谋疚募?/li>
- ks_file 即輸出文件
隨后驶赏，直接運(yùn)行即可
wgdi -ks total.conf
計(jì)算需要一點(diǎn)時(shí)間，畢竟是大量基因?qū)染稀４篌w結(jié)果如下煤傍。
注意到，
wgdi同時(shí)給出了ng86和yn00的計(jì)算結(jié)果损趋，用戶可以自己選擇患久。如果基因?qū)λ悴怀鰜韐s值，那么默認(rèn)為 -1浑槽。使用
wgdi好處之一是支持?jǐn)帱c(diǎn)續(xù)跑蒋失。只要用戶保持輸出文件名字不變化，那么即使運(yùn)行中斷桐玻，下一次運(yùn)行篙挽，wgdi會(huì)從斷點(diǎn)開始計(jì)算，節(jié)省時(shí)間镊靴。針對(duì)來源不同的輸入基因?qū)ξ募晨ǎ鄳?yīng)計(jì)算結(jié)果都會(huì)輸出至用戶指定文件。
精細(xì)的共線性區(qū)塊分析 - BlockInfo
共線性和ks值計(jì)算結(jié)果信息量極大偏竟，事實(shí)上煮落，并不容易展示。此外踊谋，嚴(yán)謹(jǐn)?shù)谋容^基因組分析項(xiàng)目上蝉仇，我們常常需要對(duì)共線性結(jié)果進(jìn)行篩選。因此殖蚕，將相關(guān)信息匯總成表轿衔，會(huì)便于刪選。
wgdi中直接提供這一整合信息的功能睦疫。
首先獲取功能參數(shù)
wgdi -bi ? >> total.conf
調(diào)整對(duì)應(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
運(yùn)行即可
wgdi -bi total.conf
輸出結(jié)果解讀害驹，直接在本地使用 Excel 打開即可
結(jié)果文件中內(nèi)容較多，但我們不需要被嚇到蛤育。因?yàn)楦娴男畔⑼鸸伲梢詭椭覀冏龈鼫?zhǔn)確的判斷葫松。以下，我們逐列說明具體含義：
-
id即共線性的結(jié)果的唯一標(biāo)識(shí) -
chr1,start1,end1即參考基因組（點(diǎn)圖的左邊）的共線性范圍 -
chr2,start2,end2即參考基因組（點(diǎn)圖的上邊）的共線性范圍 -
pvalue即共線性結(jié)果評(píng)估摘刑，常常認(rèn)為小于0.01的更合理些 -
length即共線性片段長度 -
ks_median即共線性片段上所有基因?qū)?code>ks的中位數(shù)（主要用來評(píng)判ks分布的） -
ks_average即共線性片段上所有基因?qū)?code>ks的平均值 -
homo1,homo2,homo3,homo4,homo5與multiple參數(shù)有關(guān)进宝，即一共有homo+multiple列
主要規(guī)則是：基因?qū)θ绻邳c(diǎn)圖中為紅色，賦值為1枷恕，藍(lán)色賦值為0党晋，灰色賦值為-1（顏色參考前述wgdi點(diǎn)圖 相關(guān)推文）。共線性片段上所有基因?qū)x值后求平均值徐块，這樣就賦予該共線性一個(gè)-1,1的值未玻。如果共線性的點(diǎn)大部分為紅色，那么該值接近于1胡控“饨耍可以作為共線性片段的篩選。 -
block1,block2分別為共線性片段上基因order的位置昼激。 -
ks共線性片段的上基因?qū)Φ?code>ks值 -
density1,density2共線性片段的基因分布密集程度庇绽。值越小表示稀疏
神之一手 - 使用 ks 對(duì) dotplot 著色
純粹的點(diǎn)圖，我們只能模糊地看到片段是否復(fù)制更或者復(fù)制了多少次橙困。如果將 ks 和 共線性 在點(diǎn)圖中展示瞧掺，那么就可以清晰看到共線性上ks值的變化（Ks值，事實(shí)上又對(duì)應(yīng)了一定尺度上的演化時(shí)間）凡傅。
老規(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 的影響夏跷，近期加倍的片段影響較小哼转，如果對(duì)古老多倍化的某些共線性感興趣，強(qiáng)烈建議修改lens文件(去除一些染色體槽华，保留感興趣的染色體)壹蔓，針對(duì)性的計(jì)算dotplot和ks dotplot。從
ks dotplot可以清楚看到猫态，共線性上ks值的變化庶溶。你會(huì)發(fā)現(xiàn)，近期加倍產(chǎn)生的共線性片段懂鸵，ks的波動(dòng)一般很小（有時(shí)候行疏，觀測(cè)到波動(dòng)匆光，那么就說明是染色體結(jié)構(gòu)變異區(qū)）∧鹆基于這個(gè)現(xiàn)象终息，我們?cè)谂卸╧s 峰值的時(shí)候夺巩，采用共線性的中位數(shù)要相比其他方式要準(zhǔn)確地多。另外周崭，有幸發(fā)現(xiàn)神奇的現(xiàn)象柳譬，如水稻11,12號(hào)就屬于中獎(jiǎng)了。
正確地計(jì)算 Ks 峰值
在判斷ks峰值的時(shí)候续镇，目前仍有不少分析報(bào)道直接采用所有旁系同源基因ks值進(jìn)行計(jì)算美澳。但這存在不少問題。比如串聯(lián)重復(fù)基因積累摸航，會(huì)影響Ks峰值分布制跟。基于上面的ks dotplot酱虎，我們可以直接選取近期加倍產(chǎn)生的共線性片段來計(jì)算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峰的分布圖恢暖，這個(gè)文件可以重新回到上一步 點(diǎn)圖 繪制排监，確定是否提取到感興趣的共線性區(qū)塊
得到這個(gè)圖之后，我們也可以通過點(diǎn)圖胀茵，重新確認(rèn)是否提取準(zhǔn)確社露。不斷修改參數(shù)，更或者手動(dòng)篩選琼娘，直到準(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
舉個(gè)例子，簡單修改參數(shù)脱拼，就可以得到一個(gè)穩(wěn)定的分布瞒瘸，如
block_length設(shè)定為10，homo設(shè)定為0,1熄浓，area范圍是0,2情臭。因此，葡萄近期加倍的峰值應(yīng)該為1.25赌蔑。多次調(diào)整參數(shù)俯在，可以測(cè)試峰值是否穩(wěn)定。
注意娃惯，Ks值峰值接近于0的時(shí)候跷乐，要小心，因?yàn)楣簿€性片段中有很多等于0的值趾浅。這些等于0的值如何處理愕提，如何判定分化時(shí)間馒稍，我個(gè)人沒確鑿的把握，建議重視浅侨。
寫在最后
比較基因組分析纽谒，事實(shí)上是一個(gè)非常精細(xì)的工作。軟件的使用固然很重要如输，但具體的分析思維更重要鼓黔。wgdi是我個(gè)人基于近幾年在比較基因組工作的實(shí)戰(zhàn)經(jīng)驗(yàn)開發(fā)出來的一款瑞士軍刀式軟件。幾乎覆蓋了目前絕大多數(shù)基礎(chǔ)的比較基因組分析所需功能挨决。當(dāng)然请祖，更多的功能，如祖先染色體重構(gòu)的脖祈，目前還在進(jìn)一步測(cè)試肆捕，相信很快會(huì)與大家見面。
總的來說盖高，希望通過這一系列教程慎陵，能讓大家更為全面地認(rèn)識(shí)wgdi，做出原則上無誤的比較基因組數(shù)據(jù)分析喻奥。
