写在前面
Emmm，前面推送了几期好基友-比较基因组小王子的 WGDI 软件使用和比较基因组图表解读。后台常常受到消息，大体意思“快更新！”。原本，我和小王子商量的是一周一更。但这两天小王子就要结婚了，接下一小段时间估计也忙事情。索性，今天我们把WGDI目前已成稿的内容一次推出，以飨读者。
开展批量化的 ks 计算
有比较基因组分析经验，更或者多少接触过相关领域文章的朋友应对ks有了解。ks值的大小可以反映一定时间尺度的演化时间。wgdi中ka、ks批量计算，通过paml中的yn00实现（相对于其他实现，往往更容易得到认可）。
当然，这需要安装
pal2nal.pl和paml，以及需要 在wgdi安装包目录（site-package/wgdi）的conf.ini 中配置路径。我们在前面的《安装 | 比较基因组系列之一 - WGDI 软件安装与配置》一文中已详细介绍安装和配置方式。
一切准备就绪。wgdi的操作仍然简单。
首先是获取参数
wgdi -ks ? >> total.conf
随后编辑配置文件，调整参数即可
[ks]
cds_file = Grape.cds.fa
pep_file = Grape.pep.fa
align_software = muscle
pairs_file = collinearity file
ks_file = Grape.ks
其中，
- align_software 支持 muscle 和 mafft
- pairs_file 自动支持 colinearscan 和 mcscanx 的输出结果，当然也可以是 tab 分隔的 基因对文本文件
- ks_file 即输出文件
随后，直接运行即可
wgdi -ks total.conf
计算需要一点时间，毕竟是大量基因对。大体结果如下。
注意到，
wgdi同时给出了ng86和yn00的计算结果，用户可以自己选择。如果基因对算不出来ks值，那么默认为 -1。使用
wgdi好处之一是支持断点续跑。只要用户保持输出文件名字不变化，那么即使运行中断，下一次运行，wgdi会从断点开始计算，节省时间。针对来源不同的输入基因对文件，相应计算结果都会输出至用户指定文件。
精细的共线性区块分析 - BlockInfo
共线性和ks值计算结果信息量极大，事实上，并不容易展示。此外，严谨的比较基因组分析项目上，我们常常需要对共线性结果进行筛选。因此，将相关信息汇总成表，会便于删选。
wgdi中直接提供这一整合信息的功能。
首先获取功能参数
wgdi -bi ? >> total.conf
调整对应参数
vim total.conf
[blockinfo]
blast = Grape.blastp
gff1 = Grape.gff
gff2 = Grape.gff
lens1 = Grape.len
lens2 = Grape.len
collinearity = collinearity file
score = 100
evalue = 1e-5
repeat_number = 20
position = order
ks = Grape.ks
ks_col = ks_NG86
savefile = block_information.csv
运行即可
wgdi -bi total.conf
输出结果解读，直接在本地使用 Excel 打开即可
结果文件中内容较多，但我们不需要被吓到。因为更全面的信息，可以帮助我们做更准确的判断。以下，我们逐列说明具体含义：
-
id即共线性的结果的唯一标识 -
chr1,start1,end1即参考基因组（点图的左边）的共线性范围 -
chr2,start2,end2即参考基因组（点图的上边）的共线性范围 -
pvalue即共线性结果评估，常常认为小于0.01的更合理些 -
length即共线性片段长度 -
ks_median即共线性片段上所有基因对ks的中位数（主要用来评判ks分布的） -
ks_average即共线性片段上所有基因对ks的平均值 -
homo1,homo2,homo3,homo4,homo5与multiple参数有关，即一共有homo+multiple列
主要规则是：基因对如果在点图中为红色，赋值为1，蓝色赋值为0，灰色赋值为-1（颜色参考前述wgdi点图 相关推文）。共线性片段上所有基因对赋值后求平均值，这样就赋予该共线性一个-1,1的值。如果共线性的点大部分为红色，那么该值接近于1。可以作为共线性片段的筛选。 -
block1,block2分别为共线性片段上基因order的位置。 -
ks共线性片段的上基因对的ks值 -
density1,density2共线性片段的基因分布密集程度。值越小表示稀疏
神之一手 - 使用 ks 对 dotplot 着色
纯粹的点图，我们只能模糊地看到片段是否复制更或者复制了多少次。如果将 ks 和 共线性 在点图中展示，那么就可以清晰看到共线性上ks值的变化（Ks值，事实上又对应了一定尺度上的演化时间）。
老规矩，第一步，获取参数配置文件
wgdi -bk ? >> total.conf
修改配置文件即可
vim total.conf
[blockks]
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
blockinfo = block_information.csv
pvalue = 0.01
tandem = true
tandem_length = 20
markersize = 1
area = 0,2
block_length = 5
figsize = 8,8
savefig = Grape.ks.dotplot.pdf
开始绘制
wgdi -bk total.conf
由于
repeat_number 的影响，近期加倍的片段影响较小，如果对古老多倍化的某些共线性感兴趣，强烈建议修改lens文件(去除一些染色体，保留感兴趣的染色体)，针对性的计算dotplot和ks dotplot。从
ks dotplot可以清楚看到，共线性上ks值的变化。你会发现，近期加倍产生的共线性片段，ks的波动一般很小（有时候，观测到波动，那么就说明是染色体结构变异区）。基于这个现象，我们在判定ks 峰值的时候，采用共线性的中位数要相比其他方式要准确地多。另外，有幸发现神奇的现象，如水稻11,12号就属于中奖了。
正确地计算 Ks 峰值
在判断ks峰值的时候，目前仍有不少分析报道直接采用所有旁系同源基因ks值进行计算。但这存在不少问题。比如串联重复基因积累，会影响Ks峰值分布。基于上面的ks dotplot，我们可以直接选取近期加倍产生的共线性片段来计算ks峰值。
首先，获取配置文件参数
wgdi -kp ? >> total.conf
编辑参数文件
vim total.conf
[kspeaks]
blockinfo = block_information.csv
pvalue = 0.05 # 可以用来简单筛选
tandem_length = 20 # 过滤串联重复
tandem = true # 同ks dotplot
block_length = 5 # 如调整为 10
ks_area = 0,10 # ks值的取值范围
multiple = 1
homo = -1,1 # 在`blockinfo`介绍过，可以快速筛选不同的共线性片段，如 0-1 或 0.3-1
fontsize = 9
area = 0,3 # 横轴的范围
figsize = 10,6.18
savefig = Grape.ks_median.distri.pdf
savefile = Grape.ks_median.distri.csv # 共线性中位数的数据集合，用来绘制ks峰的分布图，这个文件可以重新回到上一步 点图 绘制，确定是否提取到感兴趣的共线性区块
得到这个图之后，我们也可以通过点图，重新确认是否提取准确。不断修改参数，更或者手动筛选，直到准确为止。
[blockks]
lens1 = Grape.len
lens2 = Grape.len
genome1_name = Grape
genome2_name = Grape
blockinfo = Grape.ks_median.distri.csv
pvalue = 0.01
tandem = true
tandem_length = 20
markersize = 1
area = 0,2
block_length = 5
figsize = 8,8
savefig = Grape.ks.dotplot.filter.pdf
举个例子，简单修改参数，就可以得到一个稳定的分布，如
block_length设定为10，homo设定为0,1，area范围是0,2。因此，葡萄近期加倍的峰值应该为1.25。多次调整参数，可以测试峰值是否稳定。
注意，Ks值峰值接近于0的时候，要小心，因为共线性片段中有很多等于0的值。这些等于0的值如何处理，如何判定分化时间，我个人没确凿的把握，建议重视。
写在最后
比较基因组分析，事实上是一个非常精细的工作。软件的使用固然很重要，但具体的分析思维更重要。wgdi是我个人基于近几年在比较基因组工作的实战经验开发出来的一款瑞士军刀式软件。几乎覆盖了目前绝大多数基础的比较基因组分析所需功能。当然，更多的功能，如祖先染色体重构的，目前还在进一步测试，相信很快会与大家见面。
总的来说，希望通过这一系列教程，能让大家更为全面地认识wgdi，做出原则上无误的比较基因组数据分析。
网友评论