Ks可視化
我們在上一篇如何用WGDI進行共線性分析(一)得到共線性分析結果和Ks值輸出結果的整合表格文件后晨雳，我們就可以繪制Ks點陣圖和Ks頻率分布圖對共線性區(qū)的Ks進行探索性分析聚请，從而確定可能的Ks峰值漩怎，用于后續(xù)分析。
Ks點陣圖
最初繪制的點圖信息量很大邢隧，基本上涵蓋了歷史上發(fā)生的加倍事件所產(chǎn)生的共線性區(qū)旋廷。我們能大致的判斷片段是否存在復制以及復制了多少次膛虫，至于這些片段是否來自于同一次加倍事件則不太好確定鹃祖。借助于Ks信息 (Ks值可以反應一定尺度內(nèi)的演化時間)，我們就可以較為容易地根據(jù)點陣圖上共線性區(qū)域的顏色來區(qū)分多倍化事件宴凉。
首先誊锭，創(chuàng)建配置文件（這次是 -bk模塊，BlockKs）
wgdi -bk \? >> ath.conf
然后弥锄，修改配置文件
[blockks]
lens1 = ath.len
lens2 = ath.len
genome1_name = A. thaliana
genome2_name = A. thaliana
blockinfo = ath_block_information.csv
pvalue = 0.05
tandem = true
tandem_length = 200
markersize = 1
area = 0,3
block_length = 5
figsize = 8,8
savefig = ath.ks.dotplot.pdf
blockinfo 是前面共線性分析和Ks分析的整合結果丧靡，是繪圖的基礎蟆沫。根據(jù)下面幾個標準進行過濾
- pvalue: 共線性區(qū)的顯著性， 對應blockinfo中的pvalue列
- tandem: 是否過濾由串聯(lián)基因所形成的共線性區(qū)温治，即點陣圖中對角線部分
- tandem_length: 如果過濾饭庞，那么評估tandem的標準就是兩個區(qū)塊的物理距離
- block_length: 一個共線區(qū)的最小基因對的數(shù)量，對應blockinfo中的length列
最后運行wgdi罐盔，輸出圖片但绕。
wgdi -bk ath.conf
輸出圖片中的Ks的取值范圍由參數(shù)area控制救崔。圖中的每個點都是各個共線性區(qū)內(nèi)的基因對的Ks值惶看。不難發(fā)現(xiàn)每個共線區(qū)內(nèi)的Ks的顏色都差不多，意味著Ks值波動不大六孵，基于這個現(xiàn)象纬黎，我們在判定Ks峰值的時候，采用共線性區(qū)的中位數(shù)就比其他方式要準確的多劫窒。在圖中本今，我們大致能觀察到2種顏色，綠色和藍色主巍，對應著兩次比較近的加倍事件冠息。
Ks頻率分布圖
除了用點陣圖展示Ks外，我們還可以繪制Ks頻率的分布情況孕索。假如一個物種在歷史上發(fā)生過多倍化逛艰，那么從那個時間點增加的基因經(jīng)過相同時間的演化后，基因對之間的Ks值應該相差不多搞旭，即散怖，歸屬于某個Ks區(qū)間的頻率會明顯高于其他區(qū)間。
首先肄渗，創(chuàng)建配置文件（-kp, ksPeak）
wgdi -kp ? >> ath.conf
然后镇眷，修改配置文件
[kspeaks]
blockinfo = ath_block_information.csv
pvalue = 0.05
tandem = true
block_length = 5
ks_area = 0,10
multiple = 1
homo = -1,1
fontsize = 9
area = 0,3
figsize = 10,6.18
savefig = ath.ks_median.distri.pdf
savefile = ath.ks_median.distri.csv
這一步除了輸出Ks峰圖(savefig)外，還會輸出一個根據(jù)輸入文件( blockinfo )進行過濾后的輸出文件(savefile)翎嫡。過濾標準除了之前Ks點陣圖提及的 tandem, pvalue, block_length 外欠动，還多了三個參數(shù), ks_area, multiple, homo.
pvalue: 共線性區(qū)的顯著性， 對應blockinfo中的pvalue列惑申，pvalue設置為0.05時會保留看著很不錯的共線性片段翁垂，但是會導致古老片段的減少。
ks_area對應blockinfo中的ks列硝桩，該列記錄了共線區(qū)所有基因對的ks值沿猜。ks_area=0,10 表示只保留ks值在0到10之間的基因對。
multiple和blockinfo中的homo1,homo2,homo3,home4,homo5...列有關碗脊。一般默認為1, 表示選擇homo1列用于后續(xù)的過濾啼肩。如果改成multiple=2, 表示選擇homo2
homo用于根據(jù)共線性中基因對的總體得分(取值范圍為-1到1橄妆，值越高表明最佳匹配的基因對越多)對共線性區(qū)域進行過濾。當multiple=1, homo=-1,1時祈坠，表示根據(jù)homo1進行過濾害碾，只保留取值范圍在-1到1之間的共線性區(qū)。
最后運行程序
wgdi -kp ath.conf
從圖中赦拘，我們能夠更加直觀的觀察到2個peak慌随，基本確定2個多倍化事件。
既然得到了一個Ks的peak圖躺同，我們可以和另一款工具wgd的擬南芥分析結果進行對比阁猜。wgd的Ks值計算流程為，先進行所有蛋白之間的相互比對蹋艺，根據(jù)基因之間的同源性進行聚類剃袍，然后構建系統(tǒng)發(fā)育樹確定同源基因，最后調(diào)用PAML的codeml計算Ks捎谨，對應下圖A的A.thaliana民效。如果存在參考基因組，那么根據(jù)共線性錨點（對應下圖D）對Ks分布進行更新, 對應下圖A的A.thaliana anchors涛救。
同樣是擬南芥的分析畏邢，wgd的點圖（上圖D）信息比較雜亂，存在較多的噪聲點检吆，而WGDI的Ks點陣圖能有效的反應出不同共線性區(qū)域的Ks特點舒萎。wgd的Ks分布圖中的只能看到一個比較明顯的峰，而WGDI的分析結果能明顯的觀察到兩個咧栗。wgd在Ks上存在的問題很大一部分原因是它們是直接采用旁系同源基因計算ks逆甜，容易受到串聯(lián)重復基因積累的影響。而wgd則是基于共線性區(qū)計算Ks致板，結果更加可靠交煞，尤其是后續(xù)還可以通過不斷的調(diào)整參數(shù)，來確保Ks的峰值正確判斷斟或，這也是為什么在繪制Ks頻率分布圖的同時還要生成一個過濾后的文件素征。
題外話: wgd是我在2019年學習WGD相關分析找到的一個流程工具，那個時候雖然也看了文章里面關于軟件的細節(jié)介紹萝挤，但是由于對比較基因組學這一領域并不熟悉御毅，所以也不好評判結果的可靠性。最近在學習wgdi時怜珍，一直和開發(fā)者反復討論軟件的一些參數(shù)細節(jié)端蛆，這才知道這里面的很多細節(jié)。這也警醒我酥泛，不能追求軟件數(shù)量上的多今豆，只求用軟件跑出自己想要的結果發(fā)表文章嫌拣，失去了科研的嚴謹性。
我們會在下一節(jié)介紹如何利用Ks點陣圖和Ks頻率分布圖更可靠的擬合Ks峰值呆躲。
