【比较基因组】共线性分析（WGDI）

如果一个物种在进化过程中发生了多倍体化，那么，在基因组上就会存在一些共线性区域（即两个区域之间的基因是旁系同源基因，其基因的排布顺序基本一致）。多倍化以及后续的基因丢失和二倍化现象存在于大部分的物种中, 是物种进化的重要动力。

所以在基因组相关的paper中，经常可以看到不同物种之间共线性分析的结果，来发现一些不同物种之间保守的区域，以及一些重要的功能基因。

做共线性分析的软件有很多，以前做genome project的时候我用的MUMER。现在就尝试一下其它的一些工具，例如WGDI，jcvi，MCScanX等。这篇笔记主要介绍WGDI。

==软件安装===

雷打不动的，继续选择conda。

# 安装软件  //怕和其它的环境冲突，所以创建了一个wgdi的新环境

conda create -c bioconda -c conda-forge -n wgdi wgdi

# 启动环境

conda activate wgdi

==下载测试数据===

用的是拟南芥的数据作为测试例子。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz

==输入数据准备===

其实我个人觉得，做生信的，有很大一部分时间都在进行不同工具之间格式的转化工作，因为每个文件需求的输入是很不一样的。

wgdi需要三个输入：

分别是BLAST （outfmt为6）

基因的位置信息（ 以tab分隔，分别为chr，id，start，end，strand，order，old_id。（和GFF不太一样））

染色体长度信息和染色体上的基因个数（格式为：chr, length, gene number）

注意：对于每个基因我们只需要一个转录本，通常使用最长的转录本作为该基因的代表。

基因位置信息如下所示：

染色体长度信息和基因数目：

然后通过BLAST进行蛋白之间相互比对，输出格式选择6：

makeblastdb -in ath.pep.fa -dbtype prot

blastp -num_threads 50 -db ath.pep.fa-query ath.pep.fa -outfmt 6 -evalue 1e-5 -num_alignments 20  -out ath.blastp.txt

===共线性分析===

其实绘制点阵图非常简单了，只需要创建配置文件，然后修改配置文件，最后运行wgdi就行。

绘制点阵图

第一步：创建配置文件。

      wgdi -d \? > ath.conf

其中：

multiple = 1   # 最好的同源基因数, 用输出结果中会用红点表示

score = 100     # blast输出的score 过滤

evalue = 1e-5   # blast输出的evalue 过滤

repeat_number = 20  # genome2相对于genome1的最多同源基因数

第二步：修改配置文件。因为这个例子，我们是自我比对，所以genome1和gff1与2是一致的。

最后一步运行wgdi。

       wgdi -d ath.conf

输出的png如下图所示：

注：从图中，我们可以看出有三个颜色（红色，蓝色，灰色）。在WGDI中，红色表示genome1和genome2相似度最高的匹配，次好4个基因的是蓝色，剩余的为灰色。图中对角线出现的片段并非是自身比对，因为WGDI已经过滤掉自身比自身的结果。那么问题来了，这些基因是什么呢？

如果基因组存在加倍事件，那么对于一个基因组的某个区域可能会存在另一个区域与其有着相似的基因（同源基因），并且这些基因的排布顺序较为一致。反应到点阵图上，就是能够观察到有多个点所组成的"线"，但是当你放大之后，你会发现其实还是点而已。

由于进一步鉴定共线性区域就建立在这些"点"的基础上，那么影响这些点是否为同源基因的参数就非常重要了，即配置文件中的score,evalue,repeat_number。

共线性分析

先分清楚两个概念：synteny（同线性）和collinearity（共线性）。同线性指的是两个区域有一定数量的同源基因，对基因顺序排布无要求。共线性是同线性的特殊形式，要求同源基因的排布顺序也相似。

wgdi开发了-icl模块（Improvedversion of ColinearScan）用于进行共线性分析，使用起来和前面的点阵图一样，也是先建立配置文件。

生成配置文件：

       wgdi -icl \? >> ath.V2.conf

修改配置文件：

注：其中的evalue, score和BLAST文件的过滤有关，用于筛选同源基因对。multiple确定最佳比对基因，repeat_number表示最多允许多少基因是潜在的同源基因, grading根据同源基因的匹配程度进行打分, mg 指的是共线性区域中所允许最大空缺基因数。

wgdi -icl  ath.V2.conf

运行结束后，会得到ath.collinearity.txt，记录共线性区域。以 # 起始的行记录共线性区域的元信息，例如得分(score)，显著性(pvalue)，基因对数(N)等。

下面，我们就可以根据共线性结果来计算Ks。

同样滴，生成配置文件并编辑配置文件：

                       wgdi -ks \? >> ath.V3.conf

注：需要提供cds和pep文件，共线性分析输出文件(支持MCScanX的共线性结果）。WGDI会用muscle根据蛋白序列进行联配，然后使用pal2pal.pl 基于cds序列将蛋白联配转成密码子联配，最后用paml中的yn00计算ka和ks。

Ka和Ks在很多地方都用到，比如重测序分析的时候，用来挑选受选择的区域和基因。Ka和Ks分别指的是非同义替换位点数，和同义替换位点数。根据中性演化假说，基因的变化大多是中性突变，不会影响生物的生存，那么当一对同源基因分开的时间越早，不影响生存的碱基替换数就越多(Ks)，反之越多。

注：结果中共有6列，对应的是每个基因对的Ka和Ks值。

由于共线性和Ks值输出结果不方便后续使用和分析，所以将两者进行聚合，得到关于各个共线性区的汇总信息。

WGDI提供 -bi 参数帮助我们进行数据整合。同样也是生成配置文件，编辑配置文件。

wgdi -bi ? >> ath.V4.conf

然后运行：wgdi -bi  ath.V4.conf

生成的block.csv文件如下表所示：

输出共11列。

id 即共线性的结果的唯一标识；chr1,start1,end1 即参考基因组（点图的左边）的共线性范围；chr2,start2,end2 即参考基因组（点图的上边）的共线性范围

pvalue 即共线性结果评估（小于0.01的更合理些）；length 即共线性片段长度；ks_median 即共线性片段上所有基因对ks的中位数（主要用来评判ks分布的）；ks_average 即共线性片段上所有基因对ks的平均值；

homo...与multiple参数有关，即一共有homo+multiple列

block1,block2分别为共线性片段上基因order的位置。

ks共线性片段的上基因对的ks值

density1,density2 共线性片段的基因分布密集程度。值越小表示稀疏。

这个表格是我们后续分析的基础数据，比如，我们可以提取对角线上共线性区域的基因list。

进一步观察，可以发现这些基因的编号很多都是挨着的，说明这些都是基因簇（gene cluster）。

Ks点阵图

刚开始做的点图信息量很大，基本上涵盖了历史上发生的加倍事件所产生的共线性区。我们能大致的判断片段是否存在复制以及复制了多少次，至于这些片段是否来自于同一次加倍事件则不太好确定。借助于Ks信息 (Ks值可以反应一定尺度内的演化时间)，我们就可以较为容易地根据点阵图上共线性区域的颜色来区分多倍化事件。

和前面一样，也是创建配置文件，并编辑配置文件。

wgdi -bk \? >> ath.V5.conf

其中：blockinfo 是前面共线性分析和Ks分析的整合结果。 

pvalue: 共线性区的显著性， 对应blockinfo中的pvalue列

tandem: 是否过滤由串联基因所形成的共线性区，即点阵图中对角线部分

tandem_length: 如果过滤，那么评估tandem的标准就是两个区块的物理距离

block_length: 一个共线区的最小基因对的数量，对应blockinfo中的length列

然后运行：wgdi -bk ath.V5.conf

输出图片中的Ks的取值范围由参数area控制。图中的每个点都是各个共线性区内的基因对的Ks值。不难发现每个共线区内的Ks的颜色都差不多，意味着Ks值波动不大。因此，我们在判定Ks峰值的时候，采用共线性区的中位数就比其他方式要准确的多。在图中，我们大致能观察到2种颜色，绿色和蓝色，对应着两次比较近的加倍事件。

Ks频率分布图

除了用点阵图展示Ks外，我们还可以绘制Ks频率的分布情况。假如一个物种在历史上发生过多倍化，那么从那个时间点增加的基因经过相同时间的演化后，基因对之间的Ks值应该相差不多，即，归属于某个Ks区间的频率会明显高于其他区间。

一样的创建配置文件，修改配置文件。

wgdi -kp ? >> ath.V6.conf

最后运行：wgdi -kp ath.V6.conf

这一步除了输出Ks峰图（上图）外，还会输出一个根据输入文件( blockinfo )进行过滤后的输出文件(savefile)。过滤标准除了之前Ks点阵图提及的 tandem, pvalue, block_length 外，还多了三个参数, ks_area,multiple, homo.

pvalue: 共线性区的显著性， 对应blockinfo中的pvalue列，pvalue设置为0.05时会保留看着很不错的共线性片段，但是会导致古老片段的减少。

ks_area对应blockinfo中的ks列，该列记录了共线区所有基因对的ks值。ks_area=0,10 表示只保留ks值在0到10之间的基因对。 

multiple和blockinfo中的homo1,homo2,homo3,home4,homo5...列有关。一般默认为1, 表示选择homo1列用于后续的过滤。如果改成multiple=2,表示选择homo2

homo用于根据共线性中基因对的总体得分(取值范围为-1到1，值越高表明最佳匹配的基因对越多)对共线性区域进行过滤。当multiple=1, homo=-1,1时，表示根据homo1进行过滤，只保留取值范围在-1到1之间的共线性区。

从Ks频率分布图中，我们可以观察到2个peak，基本确定2个多倍化事件。

本文使用 文章同步助手 同步