基因组共线性工具MCScanX使用说明

基因组共线性工具MCScanX使用说明

简介

MCScanX工具集对MCScan算法进行了调整，用于检测共线性和同线性区域，还增加了可视化和下游分析。。MCscanX有三个核心工具,以及12个下游分析工具

实际使用

第一步: BLASTP比对

根据"MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity"，在使用BLASTP时需要注意如下事项：

• 如果一个基因有多个转录本，只选择注释中的第一条
• 寻找同源序列时，每个基因都要和自身以及其他基因组比较
• BLASTP的输出格式为Tabular(-outfmt 6)
• 保留前5条最佳比对序列(num_alignments 5)
• E值的阈值设置为1e-5(-evalue 1e-5)

注：上面的选项基于NCBI BLAST+

以拟南芥和水稻为例，拟南芥和水稻的基因组和注释文件都可以在EnsemblPlants上

# arabidopsis thaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.41.gff3.gz
# rice
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/oryza_sativa/pep/Oryza_sativa.IRGSP-1.0.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.41.gff3.gz

在建立BLAST的索引之前，我们需要对序列进行预处理，仅保留每个基因中的一个转录本。由于不同物种的基因命名方法不同，所以得先用less查看GFF文件，对其有一个总体的了解。

例如拟南芥的编号都是ATXGXXXXX.Y, 并且后面".Y"都是从1开始， 而水稻的是OsXXgXXXXXX表示基因，OsXXtXXXXXX-0Y表示转录本, 而"0Y"，一些转录本是只有00，而另一些则是从01开始，有02，03等。我写了一个脚本用于提取每个基因对应的最长转录本以及它的位置

我在提取的同时删除了一些不需要的字符串

zcat Arabidopsis_thaliana.TAIR10.41.gff3.gz| python get_the_longest_transcripts.py | sed -e 's/gene://' -e 's/transcript://'> at_lst_gene.txt
zcat Oryza_sativa.IRGSP-1.0.41.gff3.gz | python get_the_longest_transcripts.py | sed -e 's/gene://' -e 's/transcript://' > os_lst_gene.txt

同时考虑到我们也不一定需要叶绿体和线粒体的信息，所以还可以进一步删掉这些数据

sed -i -e '/Mt/d' -e '/Pt/d' at_lst_gene.txt
sed -i -e '/Mt/d' -e '/Pt/d' os_lst_gene.txt

使用seqkit grep提取蛋白质序列.

seqkit grep -j 20 -f <(cut  -f 2 at_lst_gene.txt ) Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | sed 's/\..*//' > at.fa
seqkit grep -j 20 -f <(cut  -f 2 os_lst_gene.txt) Oryza_sativa.IRGSP-1.0.pep.all.fa.gz | sed 's/-.*//'> os.fa

随后用makeblastdb建立索引

makeblastdb -in at.fa -dbtype prot -out index/at -parse_seqids

最后进行蛋白比对

blastp -query os.fa -db index/at -out at_rice.blast -evalue 1e-5 -num_threads 100 -outfmt 6 -num_alignments 5

注: 如果要做基因组组内和组间的共线性，那么就要将这两个基因组先进行合并, cat at.fa rice.fa > all.fa, 然后用all.fa建索引，用all.fa进行比对

第二步：构建gff

MCscanX要求的gff文件和标准的gff文件不一样，它只有四列, 其中"sp#"的sp意味着你要用2个字母代表物种，#则表示是哪条染色体。而"gene"则要是你蛋白序列的基因名

sp# gene    starting_position  ending_position

基于我的get_the_longest_transcripts.py脚本，很容易用awk进行创建该信息

awk 'BEGIN{OFS="\t"} {print "at"$3,$1,$4,$5}' at_lst_gene.txt > at_os.gff
awk 'BEGIN{OFS="\t"} {print "at"$3,$1,$4,$5}' os_lst_gene.txt >> at_os.gff

第三步: MCScanX寻找共线性区块

在经过上面两步后，后续的MCScanX基本上是瞬间完成

MCScanX ./at_rice

如果你在输出信息中发现 0 matches imported (xxxxx discarded), 那么一定是你的GFF文件里的基因名和blast结果的基因名不对应导致

输出文件分为两个：

第一个是at_rice.collinearity， 记录着共线性区块(collinear blocks)

共线性文本文件

第二个则是以网页的形式展示共线性情况。

网页展示

说明：第一列是每个基因座位的重复深度，第二列是参考基因（串联重复以红色标注，仅在序列既用于建立索引又用于搜索时才出现），之后的每一列都是匹配的共线性区块

下游分析

MCSCANX提供了一些列Java程序和Perl脚本用于下游分析，建议在Perl脚本的开头加入#!/usr/bin/env perl，此外部分软件需要安装bioPerl。

• group_collinear_genes.pl
• add_ka_and_ks_to_collinearity.pl
• detect_collinearity_within_gene_families.pl
• origin_enrichment_analysis.pl
• detect_collinear_tandem_arrays
• dissect_multiple_alignment

MCSCANX自带了一些结果可视化相关的工具(如下），可以快速做图了解大致情况，但是画的图不好看，后续会尝试自己用R语言实现作图。

• dot_plotter.java
• dual_synteny_plotter.java
• circle_plotter.java
• bar_plotter.java
• family_circle_plotter.java
• family_tree_plotter.java

软件安装

一般而言，按照下面的方法就能顺利安装。显然我的过程并没有那么顺利，所以我得多说几句

wget http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip
unzip MCScanX.zip
cd MCScanx
make
# 在.bashrc添加下有分析java类的路径, /path/to 一定要替换成自己的实际路径
CLASSPATH=/path/to/MCScanx/downstream_analyses

报错1: "msa.cc:289:22: error: ‘chdir’ was not declared in this scope"
解决方案: 打开msa.cc，在顶部加上#include <unistd.h>

报错2: "dissect_multiple_alignment.cc:252:44: error: ‘getopt’ was not declared in this scope"
解决方案: 打开"dissect_multiple_alignment.cc"，在顶部加上#include <getopt.h>

报错3: "detect_collinear_tandem_arrays.cc:286:46: error: ‘getopt’ was not declared in this scope"
解决方案： 打开"detect_collinear_tandem_arrays.cc"，在顶部加上#include <getopt.h>

报错4: "make[1]: javac: Command not found"
解决方案: 在https://www.oracle.com/technetwork/java/javase/downloads/index.html下载JDK，安装Java环境
经历了以上挫折后终于编译成功，下一步就是测试软件能否顺利运行。在MCScanx文件夹下有一个data目录, 里面准备了一些测试数据

./MCScanx data/at
duplicate_gene_classifier data/at
detect_collinear_tandem_arrays -g data/at.gff -b data/at.blast -c data/at.collinearity -o data/at.syn

脚本：get_the_longest_transcripts.py, 输出文件分为: 基因ID, 最长的转录本ID，染色体编号，start， end， strand