同源染色体片段检测工具ColinearScan的使用

作者: Boer223 | 来源:发表于2019-01-17 11:35 被阅读154次

ColinearScan

写在前面

更好的阅读体验请移步这里。
确定染色体同源片段是基因组学研究的一个重要方面，有助于揭示基因组在历史上发生的多种多样的进化事件，如DNA复制、染色体重排、基因丢失等。研究发现，谷物之间、哺乳动物之间、分属不同种的酵母之间都存在大规模的染色体同源片段；物种内部也常发现由于大规模基因组复制而形成的同源片段；约80％的拟南芥基因组处于复制区，分析表明，拟南芥的进化过程中，至少发生了一次或三次多倍化事件。正是由于这种大规模的复制事件，大量染色体发生重排、丢失和插入等，从而使染色体片段间的关系变得面目全非，对同源染色体片段的判断造成了困难，因此便有了这样一种有效的检测工具——ColinearScan。

由于该软件的主要作者之一是本科时的学院院长，故在本科时的科研中使用较多，所以在此介绍一下基本用法。

简介

ColinearScan主要基于动态规划算法并采用基因共线性方法来检测同源染色体片段，该软件的主要特点是：

合理的统计推断
较强的适应性
计算的高效性
参数选择方面，尤其是相邻同源基因对距离的确定，依据基因组的特点做了合理的分析，对推断的同源区显著性也做了深入的统计学检验

之前的共线性的方法有一些缺陷，最大问题在于参数选择基于经验，没有深入合理的理论分析。例如相邻同源基因对之间的距离是一个重要参数，经验方法难以取定一个合适的值，而把不适当的值用于寻找同源区域，会使结果严重地偏离实际情况。

安装

安装环境

C++ 编译器和库

如果是从源码编译软件来安装，就需要与ANSI C++兼容的编译器，比如GNU GCC。目前，ColinearScan已经使用GCC 3.0（X86和AMD64平台）和GCC 4.0成功编译。

ColinearScan使用了C++标准库（STL），因此建议使用与GCC一致的STL。
PERL解释器和模块

PERL解释器，并将其加入环境变量或安装在 /usr/bin 目录下。

BioPerl模块，某些脚本可能需要该模块。

其它PERL模块，Getopt::Long和Pod::Usage，详细信息可搜索CPAN。
Make程序

如果是从源码编译安装，建议使用GNU make程序。

安装步骤

源码

下载源码压缩包。
解压缩。

gzip -dc ColinearScan-1.0.1.tar.gz | tar xv
cd ColinearScan

配置ColinearScan 的安装目录，下面以安装到/home/chip/tools/ColinearScan为例，如果没有改变安装位置，则默认安装到/usr/local目录。

./configure --prefix=/home/chip/tools/ColinearScan

编译并安装。

make
make install

将可执行程序路径加入环境变量。

echo 'export PATH=$PATH:/home/chip/tools/ColinearScan/bin' >> ~/.bashrc

预编译文件

官网提供了Linux X86和Linux AMD64两个平台的预编译安装文件，如果想下载，要先确定自己的系统是否适合。

下载适合自己的安装包。
解压缩安装包。

gzip -dc ColinearScan-1.0.1-i686-linux.tar.gz | tar xv

将可执行程序路径加入环境变量。

echo 'export PATH=$PATH:/home/chip/tools/ColinearScan/bin' >> ~/.bashrc

RPM包

RPM包是安装Linux软件最简单的方式。需要注意的是，官网提供的安装包是基于Fedora Core 4 X86平台的。

下载安装包。
安装RPM包需要root权限，因此需要先查看当前是否是root用户。

whoami

如果不是，需切换为root。

安装。

rpm -ivh ColinearScan-1.0.1-i686-fedora.rpm

运行步骤

下面的步骤展示了如何使用CoinearScan来检测拟南芥的2号和籼稻的5号染色体之间的共线性片段。

BLAST

BLAST用于获得2条染色体之间的同源基因对。拟南芥的2号和籼稻的5号染色体之间的BLAST比对结果文件已经准备好，即ath_chr2_indica_chr5.blast。

提取BLAST结果中的基因对

提取BLAST结果中符合要求的基因对，相关参数有两个，score和evalue，下面是提取score > 100的示例：

cat ath_chr2_indica_chr5.blast | get_pairs.pl --score 100 > ath_chr2_indica_chr5.pairs

屏蔽高度重复位点

由于连续的单基因复制而产生的高度重复位点会使得共线性区域检测变得十分困难。因此，可通过repeat_mask.pl程序来屏蔽那些高度重复位点。如果某些位点的基因对数量超过了设置的参数，则将其视为高度重复区域并将其从基因对文件中屏蔽。

cat ath_chr2_indica_chr5.pairs | repeat_mask.pl -n 5 > ath_chr2_indica_chr5.purged

估计最大间隙长度

最大间隙长度（mg）是检测共线性区域的重要参数。max_gap.pl程序可以使用基因对文件来估计mg值。需要注意的是，估计mg值还需要染色体长度文件（详细的文件格式介绍，请参阅文件格式部分）。

max_gap.pl --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged

检测共线性片段

当有了mg值后，就可以用它来检测染色体之间的共线性区域。

block_scan.pl --mg 321000 --mg 507000 --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged

为了提高效率，上述流程也可以写为一个shell脚本，代码如下：

#!/bin/sh

do_error()
{
 echo "Error occured when running $1"
 exit 1
}

echo "Start to run the working example..."
echo

echo "* STEP1 Extract pairs from BLAST results"
echo "  We should parse BLAST results and extract pairs of anchors (genes in this example) satisfying our rule (score >= 100)."
echo
echo "  > cat ath_chr2_indica_chr5.blast | get_pairs.pl --score 100 > ath_chr2_indica_chr5.pairs"
echo
cat ath_chr2_indica_chr5.blast | get_pairs.pl --score 100 > ath_chr2_indica_chr5.pairs || do_error get_pairs.pl
echo

echo "* STEP2 Mask highly repeated anchor"
echo "  Highly repeated anchors which are mostly generated by continuous single gene duplication events make those colinear segements vague to be detected. We mask them off using a very simple algorithm."
echo
echo "  > cat ath_chr2_indica_chr5.pairs | repeat_mask.pl -n 5 > ath_chr2_indica_chr5.purged"
echo
cat ath_chr2_indica_chr5.pairs | repeat_mask.pl -n 5 > ath_chr2_indica_chr5.purged || do_error repeat_mask.pl
echo

echo "* STEP3 Estimate maximum gap length"
echo "  Use pair files with repeats masked to estimate mg values which will be used to detected colinear blocks."
echo
echo "  > max_gap.pl --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged"
echo
max_gap.pl --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged || do_error max_gap.pl
echo

echo "* SETP4 Detect blocks from pair file(s)"
echo "  Everything's ready do scan at last."
echo
echo "  > block_scan.pl --mg 321000 --mg 507000 --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged"
echo
block_scan.pl --mg 321000 --mg 507000 --lenfile ath_chrs.lens --lenfile indica_chrs.lens --suffix purged || do_error block_scan.pl
echo

echo "Now ath_chr2_indica_chr5.blocks contains predicted colinear blocks."

文件格式

基因对文件

基因对文件包含在2条染色体之间或1条染色体内发现的基因对。文件的每一行是一对基因或标记位点的记录。

#anchor orientation position anchor orientation position
OsIBCD007232 1 26277701 OsIBCD015822 -1 32139191

染色体长度文件

染色体长度文件包含物种的染色体长度。每行对应于物种的一条染色体，包括2个部分。首先是染色体号，第二个是长度（以bp或cM表示）。

5       31202585

共线性区域文件

在该文件中包括预测的共线性区域。每个区域都有一个描述行（第一行），显示用于搜索共线性区域的mg值。

+++++++++++++++++ MAXIMUM GAP LENGTH 500000 500000

每个共线性区域包含多行，包括区域的首行，中间几行，以及p值的统计行。

the 1th path length 15
OsIBCD029268 2.01061e+07 LOC_Os10 1.81952e+07 1
OsIBCD029264 2.00692e+07 LOC_Os10 1.83001e+07 -1
OsIBCD029250 1.99637e+07 LOC_Os10 1.85672e+07 1
OsIBCD029222 1.97822e+07 LOC_Os10 1.89026e+07 -1
OsIBCD029198 1.95817e+07 LOC_Os10 1.92271e+07 -1
OsIBCD029186 1.95031e+07 LOC_Os10 1.95171e+07 -1
OsIBCD029169 1.93285e+07 LOC_Os10 1.96994e+07 -1
OsIBCD029166 1.93003e+07 LOC_Os10 1.97443e+07 -1
OsIBCD044787 1.92408e+07 LOC_Os10 1.97995e+07 -1
OsIBCD029122 1.88742e+07 LOC_Os10 2.01553e+07 1
OsIBCD029076 1.84667e+07 LOC_Os10 2.0168e+07 1
OsIBCD029060 1.83421e+07 LOC_Os10 2.03553e+07 1
OsIBCD028992 1.78464e+07 LOC_Os10 2.04401e+07 1
OsIBCD028983 1.77888e+07 LOC_Os10 2.06978e+07 1
OsIBCD044725 1.73026e+07 LOC_Os10 2.07117e+07 -1
>LOCALE p-value : 7.78723e-15

其它相关软件

MCScanX

MCScan是一个通过扫描多个基因组或亚基因组以识别染色体的同源区域的工具，并可以使用基因位点来比对这些区域。MCScanX工具包对原有的MCScan工具包做了相应的调整，以用于检测共线性区域，并整合了15个程序来对结果进行可视化和进一步的分析。

image

SynFind

SynFind是在一个基因组中检测给定基因组的共线性区域。在分析时，SynFind可以检测待搜索基因组中所有基因相关的共线性区域，并且可以下载完整的共线性基因集，同时可以生成同线性深度表以获得查询基因组和每个目标基因组之间的多倍化水平。