前言
点阵图（dotplot）作为比较基因组最基础的分析，常用的分析软件有JCVI、mummer等，这里介绍一种新的点阵图方法，即WGDI画点阵图。软件安装见小编另一篇博客WGDI软件（一）：安装与配置。
小编这里选用雷公藤基因组作为研究示例对象，雷公藤除了古老的γ-WGD事件，最近还发生过一次三倍化（WGT）事件，有兴趣参考文献原文Genome of Tripterygium wilfordii and identification of cytochrome P450 involved in triptolide biosynthesis。这里主要参考了该软件开发大神的博客，感谢，详细见 点阵图 | 比较基因组分析之基石 - 手把手入门教程 。
准备文件
（1）基因组和gff3注释文件 #ascp下载工具见NCBI数据下载工具：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往往只需要研究各个染色体，所以小编将23条染色体外的scaffold都去除了，并将染色体ID都简化为Chr01-Chr23，最后，基因组解压命名为Twi.genome.fa，gff3解压命名为Twi.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' Twi.genome.fa >Twi.chr.length
然后统计每条染色体上的基因数目（gff文件中有gene或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
然后合并两个文件
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，第二列为染色体长度，第三列为基因数目
图1.Twi.len
然后准备处理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.Twi.gff
（2）blast比对文件
提取雷公藤的蛋白序列文件
gffread Twi.gff3 -g Twi.genome.fa -y Twi.pep.fa #gffread为cufflinks软件的一个工具
perl -i -lpe 's/\.$// unless /^>/' Twi.pep.fa # 清除终止密码子
使用 DIAMOND软件进行BLASTP比对，该软件比对非常快，且可获得与BLAST比较一致的结果
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参数输出文件即为blast软件的m8输出文件
dotplot图绘制
（1）WGDI 所有的分析，从一个配置文件开始。对于点阵图，我们可以通过运行以下脚本得到配置文件模板
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）运行
wgdi -d total.conf #运行脚本
脚本输出如下就表面运行成功
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.
主要参数解读
evalue、score：balstp比对结果过滤阈值，对应.blastp文件中的倒数第二、倒数第一列
repeat_number：根据blastp结果，每个基因最多取多少个同源基因，取的越多图片得到的背景点越多
savefig：添加后缀可以设置成.png、.pdf、.svg等格式
markersize：点的大小
figsize：图片大小，可根据染色体多少调整
结果解读
将Twi.dot.png下载到本地，结果展示如下图3
Twi.dot
WGDI 点图规则：根据左侧基因组的基因，最优同源（相似度最高）的点为红色，次好的四个基因为蓝色，剩余部分（同源基因有限制个数）为灰色。
结论
（1）物种自身比对，需横向看，可以看到红色点很集中，每条染色体基本上能找到两段明显红色点线，可以认为每个片段与两个片段很相似，也就是有三个相似片段，可以判定该物种近期发生了WGT事件，与文献一致
（2）仔细观察还可以发现很多蓝色点线，较红色点相对不集中，这个是由古老的γ-WGD事件引起的
（3）对角线上，本不应该出现片段。自身比自身显然是最好匹配，这些点（WGDI）已经去掉了。所以，其由串联重复造成的，即该基因临近位置的基因相似度很高，为同源基因，打在了对角线附近。可以通过计算这些串联重复的ks值，估算重复片段的爆发时间
该博客完全参考博客 点阵图 | 比较基因组分析之基石 - 手把手入门教程 ，有兴趣可以查阅
网友评论