用GENESPACE建立基于orthologues的染色体共线性关系图

发表于 分类于
GENESPACE使用记录。

为什么选择GENESPACE

目前多数建立染色体共线性分析的软件都是直接基于核酸序列比对结果的，这种分析在用于亲缘关系非常近的物种时是适用的，但如果想要分析属或以上分类关系的物种时往往会因为核酸差异太大而得不到足够连续的synteny block，可视化效果非常差。

GENESPACE的分析利用的是orthologues而不是核酸序列本身的比对相似性，因此可以在亲缘关系较远的物种之间找到较为连续的保守区块，从而得到不错的可视化效果。下面是官网教程给出的几种亲缘关系很远的物种染色体共线性分析结果：

Riparian plotting with GENESPACE (htmlpreview.github.io)

 

安装与运行（以及过程中的debug记录）

官网给出的其中一种方法是通过conda单独创建一个环境，将包括orthofinder在内的所有dependency全部装在这里面。但我在安装的时候遇到了很多很多问题，因为官方维护的问题所以在bioconda和conda-forge能够下载到的dependency版本存在大量的不兼容，最后我折腾出了一个折中的方案，即：

①建立一个orthofinder专用的conda环境；

②建立一个装有R 4.2.3的conda环境，将GENESPACE这个R包需要的其他R包装在这个环境中；

③在实际使用的时候激活orthofinder的环境，并直接在那里跨环境调用R 4.2.3；

④在conda以外手动安装MCScanX。

感觉相当于曲线救国了，整个安装调试下来差不多折腾了两个月。不过这样做的好处是不会轻易搞乱早已调试好的环境里的dependency，以防之后出错。

调试与初始化

在R环境下：

1
2
3
4
5
library(GENESPACE)

wd <- "/GENESPACE/workingDirectory"
path2mcscanx <- "/home/user/software/MCScanX/MCScanX-master"
gpar <- init_genespace(wd = wd, path2mcscanx = path2mcscanx, nCores = 10)

此时得到了如下的报错：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Checking dependencies ...
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/tree.py:367: SyntaxWarning: invalid escape sequence '\-'
  """
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/tree.py:1422: SyntaxWarning: invalid escape sequence '\-'
  """
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py:54: SyntaxWarning: invalid escape sequence '\['
  _ILEGAL_NEWICK_CHARS = ":;(),\[\]\t\n\r="
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py:57: SyntaxWarning: invalid escape sequence '\['
  _NHX_RE = "\[&&NHX:[^\]]*\]"
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py:58: SyntaxWarning: invalid escape sequence '\d'
  _FLOAT_RE = "[+-]?\d+\.?\d*(?:[eE][-+]\d+)?"
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py:60: SyntaxWarning: invalid escape sequence '\['
  _NAME_RE = "[^():,;\[\]]+"
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py:337: SyntaxWarning: invalid escape sequence '\s'
  MATCH = '%s\s*%s\s*(%s)?' % (FIRST_MATCH, SECOND_MATCH, _NHX_RE)
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/probroot.py:10: SyntaxWarning: invalid escape sequence '\i'
  """
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/probroot.py:201: SyntaxWarning: invalid escape sequence '\l'
  """
/home/user/miniconda3/envs/orthofinder/bin/scripts_of/probroot.py:267: SyntaxWarning: invalid escape sequence '\l'
  """
        Found valid path to OrthoFinder v2.55: `orthofinder`
        Found valid path to DIAMOND2 v2.19: `diamond`
        Found valid MCScanX_h executable: `/home/user/software/MCScanX/MCScanX-master/MCScanX_h`

个人感觉这里不是软件本身的问题，而是因为我在安装的时候跨了太多个环境，中途更换的dependency版本可能导致python版本的不兼容，在调用tree.py、newick.py和probroot.py这三个脚本时检测的raw strings字符串包含了“\”这种正则符号，因此无法正常运行。最后的解决方式就是挨个更改这三个python脚本，在出现报错那几行的“”符号前加上一个“r”，代表这后面引号框起来的都是纯字符串而不是正则符，以防被系统检测为syntax error。具体例子：

1
2
3
4
5
6
7
vim /home/user/miniconda3/envs/orthofinder/bin/scripts_of/newick.py

#原先报错的第54行是：
_ILEGAL_NEWICK_CHARS = ":;(),\[\]\t\n\r="

#把它更改成：
_ILEGAL_NEWICK_CHARS = r":;(),\[\]\t\n\r="

这样就可以解决报错问题了。

输入文件准备（最麻烦的一步）

GENESPACE本身是有一个parse_annotations()的功能的。但是软件预设能够检测的注释格式很少，在教程里只有“ncbi”和“phytozome”两个格式。我用到的几个基因组都是从ensembl下载的，注释格式和GENESPACE能检测的不一致，所以需要自己手动转换格式。在转换过程中要格外注意检查每一步的中间文件，确保得到的文件能被GENESPACE正确识别。

bed格式的注释文件

GENESPACE本身需要的bed文件格式如下：

1
2
chr1  1500  2000  gene1
chr1  2500  3200  gene2

而我本身的注释文件是ensembl的gff3格式，因此这里要用BEDOPS软件包中的gff2bed工具转换格式，并去除分类为“gene”以外的其他类别，比如“mRNA”“CDS”和“exon”之类的。

1
2
3
4
5
6
7
# Generating bed file from gff3 annotations
gff2bed < species.gff3 \
> species_raw.bed

# Reformatting bed files
awk -F "\t" '$8 == "gene" {print $1, $2, $3, $4}' species_raw.bed \
> species.bed

fasta格式的pep多肽序列

这个文件需要从下载参考基因组的地方找，有些基因组组装完成度没达到这个级别所以很可能只有核酸序列没有蛋白组序列，为了保持基因组质量的一致性，我在这里选的基因组都是本身就自带pep.fa的那些。

fasta文件有两个需要手动更正的格式问题，第一是header太长太冗余，我们需要的只有“gene:”后面接的那个编号而已；第二是在序列中包含了大量的换行，需要把它们全都删掉，变成如下的格式：

1
2
3
4
>gene1
MGASGRGAGEQQSPSSTGMLVMSECKGRDRSPSSSMMLVMSECKGRDRSPSSSMMLVMSECKGRDRSPSSSMMLVMSECKGRDRSPSSSMMLVMSECKGRDRSPSSSMMLVMSECKGRDRSPSSSM
>gene2
MLVMSECKGRDRSPSSSMRDRSPSSSMMLVMSECKGRDRSPSSRDRSPSSSMMLVMSECKGRDRSPSSRDRSPSSSMMLVMSECKGRDRSPSSRDRSPSSSMMLVMSECKGRDRSPSSRDRSPSSSMMLVMSECKGRDRSPSS

首先用seqtk软件去除fasta中的换行，之后再用awk对每个蛋白序列的header进行格式转换：

1
2
3
4
5
6
7
# Removing multiple lines in original fasta
seqtk seq -l 0 species.pep.fa \
> species.pep_raw.fa

# Parsing fasta headers
awk '$0 ~ "^>" {match($0, /gene:([0-9A-Z_\.]+)/, gene); print ">"gene[1]; next;}1' species.pep_raw.fa \
> species.fa

▲这里要注意一点，每个物种编辑好的bed和fa文件命名必须是相同的，比如species.bed对应的就只能是species.fa，否则GENESPACE也会检测不到。

工作路径与输入文件存放位置

建立一个空文件夹（比如叫“/GENESPACE/workingDirectory”），将上一步编辑好的几个bed和fa文件分别移动至里面相应的位置：

1
2
3
4
5
6
7
/workingDirectory
└─ peptide
    └─ species1.fa
    └─ species2.fa
└─ bed
    └─ species1.bed
    └─ species2.bed

运行GENESPACE.sh

1
2
3
source activate orthofinder

/home/user/miniconda3/envs/R_base_4/lib/R/bin/R --slave --vanilla < GENESPACE.R > GENESPACE.Rout

调用的R脚本GENESPACE.R

1
2
3
4
5
6
7
library(GENESPACE)

wd <- "/GENESPACE/workingDirectory"
path2mcscanx <- "/home/user/software/MCScanX/MCScanX-master"
gpar <- init_genespace(wd = wd, path2mcscanx = path2mcscanx, nCores = 10)

out <- run_genespace(gpar)

用10个cpu的大概运行时间是一个小时左右，软件会自动根据不同基因组的摆放顺序分别作图，效果还不错，后期也可以手动调整作图的参数。