因为自己的生信基础知识比较薄弱，最近在华农跟着王老师上了一些植物基因组课程，记录一下。这一次主要结合课堂内容、MCScanX官方文档以及自己的理解，演示基因组共线性工具MCScanX的用法。
1. 共线性分析
在樊龙江主编的《植物基因组学》中提到，植物起源于水生藻类，不同植物在基因水平上具有一定保守性（也就是具有一定的相同基因）。在一定的亲缘关系内，这种基因组水平上的保守性（基因组区块排列顺序保守性），也就是不同植物基因组间的共线性。
当我们拿到两个植物基因组数据，想要分析两个物种间是否存在基因进化历史、染色体结构变异、重要功能基因的插入缺失或者鉴定全基因组复制事件的时候，就可以利用一些基因组共线性分析工具进行直观的作图和分析。
做共线性分析之前需要区分两个描述基因组共线性的名词：
这里染色体打了个重点，因为所有的共线性分析都是基于染色体水平的基因组进行比较的，只有contig数据分析的话意义不大。现在的共线性研究以及开发的工具一般用的是collinearity（包括后面要说的MCScanx）。
共线性分析的原理主要都是分为三步：
- 获得基因在染色体上的位置（如基因组注释得到的gff文件）。
- 将基因的CDS/蛋白序列进行比对，得到高度相似的基因对（Anchoring）。
- 鉴定相同基因排列顺序的共线性区块（Chaining），不同软件算法不一样，这里不讨论算法只讨论如何应用。
其他就不过多介绍了，接下来主要讲讲MCScanX软件的用法。
2. 下载和编译MCScanX
按照官网克隆仓库，编译即可。MCScanX github官网
1
2
3
| git clone https://github.com/wyp1125/MCScanX.git
cd MCScanX
make
|
编译之后可以看到主文件的MCScanX、MCScanX_h和Duplicate_gene_classifier三个核心分析程序，以及downstream_analyses下游分析文件夹中的12个与下游分析有关的java和perl文件。
3. 运行MCScanX
3.1 数据预处理
MCScan支持的输入文件有两个：
- m8输出格式的BLASTP比对结果文件
- 记录染色体、基因名称以及起始和终止位点4个信息的gff文件
这里以棉花基因组和近缘物种可可基因组为例，两者都可以在NCBI上找到蛋白序列和注释的gff3文件。由于上机课中给的蛋白序列和gff文件都是处理好的，如果自己处理gff文件，总体思路是从gff文件的第3列提取’gene’关键词，从第9列分离基因名称信息，保留第1列染色体信息，保留第4列和第5列保留start和end信息即可。
这里输入的gff文件不是标准格式，简单写个python脚本处理：
1
2
3
4
5
6
7
8
9
10
11
| import os
import pandas as pd
file_path = './Theobroma_cacao.chr.gff3'
with open(file_path, "r") as file:
temp = pd.read_csv(file, sep = '\t',comment = '#', header = None)
temp = temp.drop(temp[temp[2] != 'gene'].index)
temp = temp[[0,8,3,4]]
temp[8] = temp[8].map(lambda x: x.split(';')[0].split('=')[1].split(':')[1])
temp.to_csv('Tcacao.gff3', header = None, index = None, sep = '\t')
|
当然，根据官网的表述，最好将染色体名称改为两个字母（物种缩写） + 数字（代表染色体编号）的形式。不改也没关系，两个基因组的染色体名字不要一样就好了，不改的话只有在下游分析对共线性区块分组的时候有点影响（添加的物种那列只显示两个字母）实质上不会有什么影响。
这步检查一下gene数量与pep文件的蛋白质条数是否一样，不一样的话可能是pep中有转录本蛋白序列，重新筛选完整的gff3文件，再用gffread提取蛋白序列，这里就不赘述了，如果处理有问题后续在更新。
本例中，棉花Gossypium herbaceum相关文件前缀为Gh；可可Theobroma cacao相关文件前缀为Tcacao。
1
2
3
4
5
6
7
8
9
10
| # 合并蛋白序列（方便做基因组组内和组间共线性分析）
cat Gh.pep Tcacao.pep > Gh_Tcacao.pep
# blast建库
makeblastdb -in /public/home/wlxie/biosoft/MCScanX/Data/ExerciseData/Gh_Tcacao.pep -dbtype prot -input_type fasta -out Gh_Tcacao
# 蛋白序列拆分40份（方便提交并行任务到集群）
perl fasta-splitter.pl --n-parts 40 Gh_Tcacao.pep
mkdir Gh_Tcacao
mv Gh_Tcacao.part* Gh_Tcacao
|
用到的蛋白序列拆分perl脚本来自于Kirill Kryukov，具体在哪个仓库找不到了……这里提供下载，拆分后文件夹名称如下：
上面拆分的40个蛋白序列名称数字部分是等宽的，不方便提交并行任务到集群，这里简单修改下文件名称。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| import os
path = '/public/home/wlxie/biosoft/MCScanX/Data/ExerciseData/Gh_Tcacao/'
fileList = os.listdir(path)
prefix = 'Gh_Tcacao_'
suffix = '.pep'
m = 1
for file in fileList:
old_path = path + os.sep + file
if os.path.isdir(old_path):
continue
new_path = path + os.sep + prefix + str(m) + suffix
os.rename(old_path, new_path)
m += 1
|
修改之后提交40个blastp并行任务，每个任务4核，很快就可以跑完。
1
2
3
4
5
6
7
| #!/bin/bash
#SBATCH --array=1-40
#SBATCH --cpus-per-task=4
echo start on $(date)
srun blastp -query Gh_Tcacao${SLURM_ARRAY_TASK_ID}.pep -out Gh_Tcacao_${SLURM_ARRAY_TASK_ID}.blast -db Gh_Tcacao -outfmt 6 -num_threads 4 -num_alignments 5 -evalue 1e-10
echo end on $(date)
|
最后合并blastp结果。
1
2
| cat Gh_Tcacao_*.blast > Gh_Tcacao.blast
rm Gh_Tcacao_*.blast
|
3.2 运行MCScanX和结果解读
MCScanX有三个主命令：
- MCScanX共线性分析
- MCScanX_h也是共线性分析，输入不是BLASTP文件，而是第三方检测的以制表符分隔的成对同源关系文件
- Duplicate_gene_classifier使用MCScanX的算法鉴定singleton（单基因）和重复基因
这里用MCScanX，前面处理好了的blastp结果文件和gff文件放在同一个文件夹，输入的文件相对路径要到文件名的前缀部分
1
| ./MCScanX Data/ExerciseData/Gh_Tcacao
|
结果文件如下：
.collinearity 是两个基因组共线性结果文件（默认分析collinearity），包含了本次运行的参数、计算出的共线性区块等。还可以根据这个文件用grep的方法提取两个物种之间的同源基因（提取第二列第三列物种名字不一样的行，去重复，计数），这里也不赘述了。
.html 这个文件夹每个文件对应一条染色体共线性分析结果，包括基因座的共线性区块数量、基因名称（串联重复基因会标红）和具体比对上哪些共线性区块：
.tandem文件显示基因组内所有串联重复基因对：
3.3 下游分析和作图
官方在downstream_analyses文件夹中提供了12个下游分析的perl和java脚本。我个人将这个文件下的分析工具分为两类：
- 功能相关（此部分与数据处理相关）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| # Detect_syntenic_tandem_arrays 检测串联重复基因
../../../downstream_analyses/detect_collinear_tandem_arrays -g Gh_Tcacao.gff -b Gh_Tcacao.blast -c Gh_Tcacao.collinearity -o tandem_arrays
# Dissect_multiple_alignment 将共线性区块分为物种内和物种间共线性区块
../../../downstream_analyses/dissect_multiple_alignment -g Gh_Tcacao.gff -c Gh_Tcacao.collinearity -o dissect_multiple_alignment
# add_ka_and_ks_to_collinearity.pl 计算ka和ks值（在collinearity结果后面增加两列，执行时间较长）
cat Gh.cds Tcacao.cds > Gh_Tcacao.cds # 需要cds序列，且与collinearity文件中的序列名要一致，否则算出来全部都是-2
perl ../../../downstream_analyses/add_ka_and_ks_to_collinearity.pl -i Gh_Tcacao.collinearity -d Gh_Tcacao.cds -o add_kaks_to_synteny
# group_collinear_genes.pl 基因分组，构建基因家族（结果文件需要去掉第一行无效信息）
perl ../../../downstream_analyses/group_collinear_genes.pl -i Gh_Tcacao.collinearity -o group_collinear_genes
# 以及以下几种功能就不一一试了
# detect_collinearity_within_gene_families.pl 检测基因家族的共线性基因对，需要用到前面构建的基因家族文件
# origin_enrichment_analysis.pl 根据Duplicate_gene_classifier的结果识别输入基因家族的重复基因起源的潜在富集？不太懂什么意思
|
- 作图相关（此部分均有对应的配置文件，在downstream_analyses文件夹操作，其他文件路径java会发生未知错误）
1
2
3
4
5
6
7
| # dot_plotter 两组染色体所有共线性区块做点图
java dot_plotter -g ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.gff -s ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.collinearity -c dot.ctl -o dot.png
# dot.ctl配置文件内容如下
800 //x轴维度（像素大小，下同，不再赘述）
800 //y轴维度
Ghir_A01,Ghir_A02,Ghir_A03,Ghir_A04,Ghir_A05,Ghir_A06,Ghir_A07,Ghir_A08,Ghir_A09,Ghir_A10,Ghir_A11,Ghir_A12,Ghir_A13 //x轴染色体名称
Chromosome_1,Chromosome_2,Chromosome_3,Chromosome_4,Chromosome_5,Chromosome_6,Chromosome_7,Chromosome_8,Chromosome_9 //y轴染色体名称
|
1
2
3
4
5
6
7
| # dual_synteny_plotter
java dual_synteny_plotter -g ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.gff -s ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.collinearity -c dual_synteny.ctl -o dual_synteny.png
# dual_synteny.ctl
600 //图宽
800 //图高
Ghir_A01,Ghir_A02,Ghir_A03,Ghir_A04,Ghir_A05,Ghir_A06,Ghir_A07,Ghir_A08,Ghir_A09,Ghir_A10,Ghir_A11,Ghir_A12,Ghir_A13 //位于左边的染色体名称
Chromosome_1,Chromosome_2,Chromosome_3,Chromosome_4,Chromosome_5,Chromosome_6,Chromosome_7,Chromosome_8,Chromosome_9 //位于右边的染色体名称
|
1
2
3
4
5
| # circle_plotter
java circle_plotter -g ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.gff -s ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.collinearity -c circle.ctl -o circle.png
# circle.ctl
800 //图宽和高
Ghir_A01,Ghir_A02,Ghir_A03,Ghir_A04,Ghir_A05,Ghir_A06,Ghir_A07,Ghir_A08,Ghir_A09,Ghir_A10,Ghir_A11,Ghir_A12,Ghir_A13,Chromosome_1,Chromosome_2,Chromosome_3,Chromosome_4,Chromosome_5,Chromosome_6,Chromosome_7,Chromosome_8,Chromosome_9 //环状图的染色体名称
|
1
2
3
4
5
6
7
| # bar_plotter
java bar_plotter -g ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.gff -s ../Data/ExerciseData/Gh_Tcacao/Gh_Tcacao.collinearity -c bar.ctl -o bar.png
# bar.ctl
800 //x轴维度
800 //y轴维度
Ghir_A01,Ghir_A02,Ghir_A03,Ghir_A04,Ghir_A05,Ghir_A06,Ghir_A07,Ghir_A08,Ghir_A09,Ghir_A10,Ghir_A11,Ghir_A12,Ghir_A13 //参考染色体
Chromosome_1,Chromosome_2,Chromosome_3,Chromosome_4,Chromosome_5,Chromosome_6,Chromosome_7,Chromosome_8,Chromosome_9 //目标染色体
|
1
2
3
| # 还有两个作图的工具也不一一试了，以后有需要再做
# family_circle_plotter 基因家族同源基因圆形图，红线连接一个基因家族所有同源基因
# family_tree_plotter 基因家族树，共线性基因对用红线连接，串联重复基因用蓝线连接
|
可以看到这个软件绘图还是非常简单方便的，不过图片质量不是很好，作图的像素点大小都是由自己定义的，因此排版缩放也很容易失真。还是比较推荐JCVI一类的软件，引入多种过滤参数功能更强大，且做的图是矢量图，比较好看……而且有docker镜像，不怕安装依赖的问题：
1
| singularity -d build jcvi.sif docker://tanghaibao/jcvi
|
以后有空再更新吧~
2023.5.8 更新
上面例子仅仅只是作图，没有阐明具体的生物学问题。这里通过以上的数据补充几个思考：
- 棉花和可可的WGD（全基因组加倍）事件分别在多少年前？
- 在多少年前棉花和可可发生了物种分化？
在解决上面两个问题前，需要知道一个基本概念——中性进化论。
在中性进化理论中，分子水平的变异是中性的，不受自然选择的影响。也就是说对于一个基因序列而言，每个位点上的演化（也就是发生突变）速率都是恒定的，如果一个基因位点发生同义突变，氨基酸序列并未发生变化，则这种突变不会影响物种的适应性。
同义突变率ds（Ks）指平均每个同义位点上发生同义置换的数目，在物种进化中代表了进化过程的背景碱基替换率，两个物种或者同个物种之间的Ks值可以通过上面的MCScanX的下游分析程序add_ka_and_ks_to_collinearity.pl计算得出。
如果一个物种发生了全基因组加倍事件，则会产生大量的旁系同源基因，反映在Ks值上就是有大量的Ks值接近的同源基因对产生，在Ks统计图中会出现一个峰（peak）；如果两个物种之间发生了物种分化，同样产生大量的直系同源基因（物种形成的伴随事件），同样是Ks值接近的同源基因产生，并且在Ks统计图出现峰值。因此，我们可以根据同一个物种内以及不同物种间的同源基因Ks值来反推物种的WGD事件和物种分歧事件，Ks峰值处发生的事件即为最近一次的物种WGD事件或物种分歧事件。
这些事件发生的时间点，如果有已知的化石证据则最准确（根据放射性同位素衰变），如果没有就需要根据分子钟理论计算对应的时间，我们用最基础的公式T=Ks/2r。Ks，即同义突变率，平均每个位点的突变次数；r是核酸突变速度，也就是这个分类的物种每个位点每年的突变概率。
解释一下这个公式怎么来的，在两个物种分化一定的时间T后，两者都以相似的速度r累积突变（分化后是近缘物种，核酸突变速度类似），则两个物种之间核酸替换率K=T*（r+r）。实际上为了避免进化选择对突变速率的影响，这里一般用同义突变率替代核酸替换率。r值是怎么计算的呢？同样要依赖化石证据，根据两个物种共同祖先的化石时间反推r值，假设同一类物种核酸突变率类似，则这个r值还可以进一步用于别的近缘物种。
以上理论依据建立在同一类物种核酸突变率类似的假设中，实际不一定如我们所愿，且化石依据也会存在一定误差。这个r值也只是估计个大概，实际上只要在10的-9次方数量级，算出来的时间能自圆其说就行。
1
2
3
4
5
| # 利用MCScanX计算物种内部的共线性结果
./MCScanX -b 1 Data/ExerciseData/Gh_Tcacao_1
# 利用MCScanX计算两个物种之间的共线性结果
./MCScanX -b 2 Data/ExerciseData/Gh_Tcacao_2
|
分别整理棉花与棉花、可可与可可、棉花与可可的基因组共线性结果，只选取最后一列ks值，用R做Ks密度分布图，可以参考现有的教程：
Ks密度曲线分布图绘图 - 简书 (jianshu.com)
Finding Peak Values For a Density Distribution (ianmadd.github.io)
找到Ks密度分布图的峰值，选一个参考文献里合适的r值，就可以计算上面两个问题了。
顺便补充一下，下面这篇教程里有一个python脚本可以提取最长转录本，因为原作者给的是图片懒得敲下来试了，看了下代码逻辑是处理pep文件，根据序列名比较同一个基因的不同转录本长度，选取最长的那个。逻辑没有问题，因为不是刚需，有需要自己再去复现，就不重复造轮子了。
WGD（全基因组复制）分析——Ka/Ks及4Dtv值计算 - 简书 (jianshu.com)
