前言
點(diǎn)陣圖（dotplot）作為比較基因組最基礎(chǔ)的分析亲茅，常用的分析軟件有JCVI颁股、mummer等，這里介紹一種新的點(diǎn)陣圖方法征绸，即WGDI畫點(diǎn)陣圖久橙。軟件安裝見小編另一篇博客WGDI軟件（一）：安裝與配置。
小編這里選用雷公藤基因組作為研究示例對(duì)象管怠，雷公藤除了古老的γ-WGD事件淆衷，最近還發(fā)生過一次三倍化（WGT）事件，有興趣參考文獻(xiàn)原文Genome of Tripterygium wilfordii and identification of cytochrome P450 involved in triptolide biosynthesis渤弛。這里主要參考了該軟件開發(fā)大神的博客祝拯，感謝，詳細(xì)見 點(diǎn)陣圖 | 比較基因組分析之基石 - 手把手入門教程 她肯。
準(zhǔn)備文件
（1）基因組和gff3注釋文件 #ascp下載工具見NCBI數(shù)據(jù)下載工具：aspera的安裝與使用
ascp -i ~/asperaweb_id_dsa.openssh -QTr -l200m anonftp@ftp.ncbi.nlm.nih.gov:genomes/all/GCA/013/401/445/GCA_013401445.1_ASM1340144v1/GCA_013401445.1_ASM1340144v1_genomic.fna.gz ./
ascp -i ~/asperaweb_id_dsa.openssh -QTr -l200m anonftp@ftp.ncbi.nlm.nih.gov:genomes/all/GCA/013/401/445/GCA_013401445.1_ASM1340144v1/GCA_013401445.1_ASM1340144v1_genomic.gff.gz ./
下載的基因組和gff3文件有23條染色體和很多其他scaffold佳头，dotplot往往只需要研究各個(gè)染色體，所以小編將23條染色體外的scaffold都去除了晴氨，并將染色體ID都簡(jiǎn)化為Chr01-Chr23康嘉，最后，基因組解壓命名為Twi.genome.fa籽前，gff3解壓命名為Twi.gff3亭珍，然后做一下處理
首先計(jì)算每條染色體長(zhǎng)度信息
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' Twi.genome.fa >Twi.chr.length
然后統(tǒng)計(jì)每條染色體上的基因數(shù)目（gff文件中有g(shù)ene或mRNA信息）
perl -lane 'next if /^#/;$count{$F[0]}++ if $F[2] eq "gene";END{print join qq{\t},$_,$count{$_} for sort keys %count}' Twi.gff3 > Twi.gene.counts
然后合并兩個(gè)文件
perl -lane 'if($flag){print join qq{\t},$_,$count{$F[0]}}else {$count{$F[0]}=$F[1]}$flag=1 if eof(ARGV)' Twi.gene.counts Twi.chr.length|sort -k1,1V >Twi.len
合并文件如下圖1所示，第一列為染色體ID枝哄，第二列為染色體長(zhǎng)度肄梨，第三列為基因數(shù)目
然后準(zhǔn)備處理gff文件
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}}' Twi.gff3 > Twi.gff
處理完的格式如下圖2所示，分別為 chr挠锥，id峭范，start，end瘪贱，stand纱控，order辆毡。其中，order是每條染色體從1開始的排序
（2）blast比對(duì)文件
提取雷公藤的蛋白序列文件
gffread Twi.gff3 -g Twi.genome.fa -y Twi.pep.fa #gffread為cufflinks軟件的一個(gè)工具
perl -i -lpe 's/\.$// unless /^>/' Twi.pep.fa # 清除終止密碼子
使用 DIAMOND軟件進(jìn)行BLASTP比對(duì)甜害，該軟件比對(duì)非巢耙矗快，且可獲得與BLAST比較一致的結(jié)果
wget https://github.com/bbuchfink/diamond/releases/download/v0.9.24/diamond-linux64.tar.gz #下載diamond
tar -pzxvf diamond-linux64.tar.gz #解壓即安裝
./diamond makedb --in Twi.pep.fa --db Twi.pep.fa
./diamond blastp --db Twi.pep.fa --query Twi.pep.fa --outfmt 6 --evalue 1e-5 --max-target-seqs 20 --out Twi.blastp # --outfmt 6參數(shù)輸出文件即為blast軟件的m8輸出文件
dotplot圖繪制
（1）WGDI 所有的分析尔店，從一個(gè)配置文件開始眨攘。對(duì)于點(diǎn)陣圖，我們可以通過運(yùn)行以下腳本得到配置文件模板
wgdi -d ? >> total.conf
total.conf配置文件修改成如下所示
[dotplot]
blast = Twi.blastp
gff1 = Twi.gff
gff2 = Twi.gff
lens1 = Twi.len
lens2 = Twi.len
genome1_name = Twi
genome2_name = Twi
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 20
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.4
figsize = 15,15
savefig = Twi.dot.png
（2）運(yùn)行
wgdi -d total.conf #運(yùn)行腳本
腳本輸出如下就表面運(yùn)行成功
blast = Twi.blastp
gff1 = Twi.gff
gff2 = Twi.gff
lens1 = Twi.len
lens2 = Twi.len
genome1_name = Twi
genome2_name = Twi
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 20
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.4
figsize = 15,15
savefig = Twi.dot.png
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.
主要參數(shù)解讀
evalue嚣州、score：balstp比對(duì)結(jié)果過濾閾值鲫售，對(duì)應(yīng).blastp文件中的倒數(shù)第二、倒數(shù)第一列
repeat_number：根據(jù)blastp結(jié)果该肴，每個(gè)基因最多取多少個(gè)同源基因情竹，取的越多圖片得到的背景點(diǎn)越多
savefig：添加后綴可以設(shè)置成.png、.pdf匀哄、.svg等格式
markersize：點(diǎn)的大小
figsize：圖片大小秦效，可根據(jù)染色體多少調(diào)整
結(jié)果解讀
將Twi.dot.png下載到本地，結(jié)果展示如下圖3
WGDI 點(diǎn)圖規(guī)則：根據(jù)左側(cè)基因組的基因涎嚼，最優(yōu)同源（相似度最高）的點(diǎn)為紅色阱州，次好的四個(gè)基因?yàn)樗{(lán)色，剩余部分（同源基因有限制個(gè)數(shù)）為灰色法梯。
結(jié)論
（1）物種自身比對(duì)苔货，需橫向看，可以看到紅色點(diǎn)很集中立哑，每條染色體基本上能找到兩段明顯紅色點(diǎn)線蒲赂，可以認(rèn)為每個(gè)片段與兩個(gè)片段很相似，也就是有三個(gè)相似片段刁憋，可以判定該物種近期發(fā)生了WGT事件滥嘴，與文獻(xiàn)一致
（2）仔細(xì)觀察還可以發(fā)現(xiàn)很多藍(lán)色點(diǎn)線，較紅色點(diǎn)相對(duì)不集中至耻，這個(gè)是由古老的γ-WGD事件引起的
（3）對(duì)角線上若皱，本不應(yīng)該出現(xiàn)片段。自身比自身顯然是最好匹配尘颓，這些點(diǎn)（WGDI）已經(jīng)去掉了走触。所以，其由串聯(lián)重復(fù)造成的疤苹，即該基因臨近位置的基因相似度很高互广，為同源基因，打在了對(duì)角線附近”怪澹可以通過計(jì)算這些串聯(lián)重復(fù)的ks值像樊，估算重復(fù)片段的爆發(fā)時(shí)間
該博客完全參考博客 點(diǎn)陣圖 | 比較基因組分析之基石 - 手把手入門教程 ，有興趣可以查閱
