寫在前面
前述，我們通過《安裝 | 比較基因組系列之一 - WGDI 軟件安裝與配置》一文，帶大伙安裝和配置了我們 WGDI 軟件。接下來，我們直切主題，從實例數(shù)據(jù)開始，手把手帶大家進行比較基因組分析，并做具體結(jié)果解讀。
比較基因組學的分析工作常常都是從一張點圖開始的。一張清晰的點圖能反應出來非常多的信息。
初探 WGDI
WGDI 所有的分析，從一個配置文件開始。對于點陣圖，我們可以通過運行
wgdi -d ? >> total.conf
進而查看配置文件模板
cat total.conf
[dotplot]
blast = blast file
gff1 = gff1 file
gff2 = gff2 file
lens1 = lens1 file
lens2 = lens2 file
genome1_name = Genome1 name
genome2_name = Genome2 name
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 20
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = savefile(.png,.pdf)
下面，逐個參數(shù)給大家解讀。
blast = blast file，即設置輸入的 blast 文件，一般使用物種A蛋白序列全集 比對到 物種B蛋白序列全集，得到blastp結(jié)果文件（后續(xù)給大家演示）-
gff1 = gff1 file，即基因位置信息文件，記錄每個基因的具體信息，gff 文件使用 Tab鍵 分割，分別為 chr，id，start，end，stand，order，old_id。其中，order是每條染色體從1開始的排序，old_id 和后面列的信息不讀取。
-
lens1 = lens1 file，染色體長度信息，記錄每條染色體的長度。三列分別為染色體號，染色體bp長度，染色體基因個數(shù)。scaffold或contig 可以等效于chr
其他參數(shù)，不重要。主要是blastp結(jié)果過濾以及出圖參數(shù)。
準備示例數(shù)據(jù)
從上述可以看出，輸入數(shù)據(jù)事實上可以直接從兩個文件開始準備。
- 基因組序列文件
- 基因結(jié)構注釋信息
文件可直接從公開數(shù)據(jù)庫 Phytozome 上下載，此處為v10的
ls -1
total.conf
Vitis_vinifera.genome.fna
Vitis_vinifera.genome.gff3
首先，準備染色體長度信息
perl -0076 -ane '@F=map{s/[>\r\n]//gr}@F;$id=shift @F;print $id,qq{\t},length (join q{},@F),qq{\n} if $id' Vitis_vinifera.genome.fna > Grape.chr.length
隨后，統(tǒng)計每條染色體上的基因數(shù)目（注意，此處統(tǒng)計的是 gene 的數(shù)目，如果你的注釋信息文件沒有 gene 的數(shù)目，那么你可能要統(tǒng)計 mRNA 的數(shù)目，并注意是否有可變剪切本等）
perl -lane 'next if /^#/;$count{$F[0]}++ if $F[2] eq "gene";END{print join qq{\t},$_,$count{$_} for sort keys %count}' Vitis_vinifera.genome.gff3 > Grape.gene.counts
合并兩個文件，得到 WGDI 所需的 len 文件
perl -lane 'if($flag){print join qq{\t},$_,$count{$F[0]}}else {$count{$F[0]}=$F[1]}$flag=1 if eof(ARGV)' Grape.gene.counts Grape.chr.length |sort -k1,1V > Grape.len
準備 gff 文件
Emmm，perl one-liner 嘛...
# 此處直接使用 mRNA，葡萄注釋信息中每個基因只對應了一個mRNA
perl -lane 'next unless $F[2] eq "mRNA";/ID=([^;]+)/;push @geneInfo,[$F[0],$1,$F[3],$F[4],$F[6]];END{$preChr="";for(sort {$a->[0] cmp $b->[0] || $a->[2] <=> $b->[2]} @geneInfo){if($preChr ne $_->[0]){$c=0;$preChr=$_->[0]};print join qq{\t},@{$_},++$c}}' Vitis_vinifera.genome.gff3 > Grape.gff
準備 BLAST 文件
本例中，我們做的是葡萄內(nèi)部的，所以只需要提取葡萄的蛋白序列文件，比對到自身即可
gffread Vitis_vinifera.genome.gff3 -g Vitis_vinifera.genome.fna -y Grape.pep.fa
# 清除終止密碼子
perl -i -lpe 's/\.$// unless /^>/' Grape.pep.fa
比對到自身，推薦 BLASTP，因為準確度也很重要。此處使用 DIAMOND，主要是為了加速
./diamond makedb --in Grape.pep.fa --db Grape.pep.fa
./diamond blastp --db Grape.pep.fa --query Grape.pep.fa --outfmt 6 --evalue 1e-5 --max-target-seqs 20 --out Grape.blastp
繪制點陣圖
修改配置文件，設置輸入的兩個文件
[dotplot]
blast = Grape.blastp
gff1 = Grape.gff
gff2 = Grape.gff
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 20
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = Grape.dot.png
隨后運行即可
wgdi -d total.conf
blast = Grape.blastp
gff1 = Grape.gff
gff2 = Grape.gff
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 20
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = Grape.dot.png
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.
查看當前目錄，可以看到輸出 png 文件Grape.dot.png，直接打開或者下載到本地打開即可
上述，
position = order 參數(shù)的意思是，基因按照排序位置可視化，我們也可以直接按照具體染色體上的物理位置進行可視化，使用 position = end
那么為什么要用物理位置可視化？這個與后續(xù) 核型分析，祖先染色體重構等相關，繼續(xù)埋雷，感興趣的就等后續(xù)推送....
與此類似，還有mutiple 參數(shù)，表示最優(yōu)匹配格式，同樣與更高維度的比較基因組分析相關。還是埋雷....哈哈
回到主題，任意修改 lens1 和 lens2 的染色體的數(shù)量和順序，可以得到任意染色體間的點圖。
比如，我們可以直接去掉 random 染色體碎片
perl -i -lne 'print unless /random/' Grape.len
點圖解讀
首先，明確 WGDI 點圖規(guī)則：根據(jù)左側(cè)基因組的基因，最優(yōu)同源（相似度最高）的點為紅色，次好的四個基因為藍色，剩余部分（同源基因有限制個數(shù)）為灰色。
- 點圖需要橫向看和縱向看：這個點圖是物種自身比對，所以只需橫向看。片段深度應該為 2，但最好同源個數(shù)為1，即紅色點沒有集中分布在某個片段上。常常可以認兩個片段很相似，再加上自身。所以，認定為最近發(fā)生的加倍為三倍化。除此之外，藍色或灰色的片段很少有，表明更古老的加倍很不明顯。
- 對角線上，本不應該出現(xiàn)片段。自身比自身顯然是最好匹配，這些點（WGDI）已經(jīng)去掉了。所以，其由串聯(lián)重復造成的，即該基因臨近位置的基因相似度很高，為同源基因，打在了對角線附近。可以通過計算這些串聯(lián)重復的ks值，估算重復片段的爆發(fā)時間。
- 最后，事實上，點圖是后續(xù)跑共線性的基石。score， evalue， repeat_number 判定的同源點的數(shù)量是共線性打分矩陣的來源。
寫在最后
往往，越是高階的分析人員，使用的工具卻越是簡陋。點圖，看似簡單 和 粗略。但其蘊含的才是真正全面的比較基因組信息。
更多內(nèi)容，讓我們一起期待后續(xù)教程。
