前段时间因为课题需求，需要获取两个物种间的所有一对一直系同源基因对（One-to-One ortholog, 为避免文章内容冗余后面统称 OTO）。由于两个目标物种在 Ensembl 上都有，所以可以直接通过 BioMart 下载其基因的同源关系并提取出来。

但经过筛选和过滤后，发现 BioMart 中能确定的 OTO 数量极低，不太能满足分析需求。考虑到我所分析的目标基因类型在不同物种间具有高顺序保守性和序列保守性，因此想看看能不能通过 JCVI 的共线性结果来获取更多的 OTO 作为补充，发现效果还挺不错，固有此文。

JCVI 介绍

JCVI 是一个用于比较基因组学分析的多功能工具包，共线性分析只是它其中的一个子功能。

JCVI 的安装，最简单的方式是通过 pip：

1
pip install jcvi

也可通过 conda 或 mamba 安装：

1
conda install bioconda::jcvi

此外也需要安装用于比对的 LAST：

1
conda install bioconda::last

关于 JCVI 的运行，更多细节可见：JCVI tutorial

请注意，JCVI 的 wiki 写的并不全面，很多参数需要灵活调整以适配自身需求。本文所用到的文件（如注释等）均来自 Ensembl 数据库，如果使用其他来源的文件，请灵活调整参数（或调整源文件格式）以达成相同目的。

JCVI 寻找 OTO 示例

本文使用到的文件下载（人和黑猩猩 Ensembl V110）：

1
2
3
4
wget ftp://ftp.ensembl.org/pub/release-110/fasta/pan_troglodytes/cds/Pan_troglodytes.Pan_tro_3.0.cds.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-110/gff3/pan_troglodytes/Pan_troglodytes.Pan_tro_3.0.110.gff3.gz
wget ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-110/gff3/homo_sapiens/Homo_sapiens.GRCh38.110.gff3.gz

下文中的命令请根据自身情况修改相关输入和输出文件名，同时检查是否需要修改参数。

①、首先使用 jcvi.formats.gff 将其 gff 文件中特定的注释类型转为 bed 文件格式：

1
2
3
4
5
6
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id --primary_only Pan_troglodytes.Pan_tro_3.0.110.gff3.gz -o chimpanzee.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id --primary_only Homo_sapiens.GRCh38.110.gff3.gz -o human.bed
python -m jcvi.formats.bed uniq chimpanzee.bed
python -m jcvi.formats.bed uniq human.bed
mv chimpanzee.uniq.bed chimpanzee.bed
mv human.uniq.bed human.bed

参数说明：

  • --type：指定想要提取的注释类型，可以通过逗号分隔同时提取多种信息。这里指定提取蛋白编码基因的 mRNA 序列区间范围。
  • --key：指定提取时使用的标志信息。由于 Ensembl CDS 文件中序列以 transcipt id 开头，因此此处使用 transcript_id 以和 CDS 文件对应。
1
2
3
4
5
6
7
8
# Ensembl CDS 文件示例
>ENST00000633790.1 cds scaffold:GRCh38:HSCHR7_2_CTG6:503684:504150:1 gene:ENSG00000282748.1 gene_biotype:TR_V_gene transcript_biotype:TR_V_gene gene_symbol:TRBV5-7 description:T cell receptor beta variable 5-7 (non-functional) [Source:HGNC Symbol;Acc:HGNC:12224]
ATGGGCCCCGGGCTCCTCTGCTGGGTGCTGCTTTGTCCCCTAGGAGAAGGCCCAGTGGAC
GCTGGAGTCACCCAAAGTCCCACACACCTGATCAAAACGAGAGGACAGCACGTGACTCTG
AGATGCTCTCCTATCTCTGGGCACACCAGTGTGTCCTCGTACCAACAGGCCCTGGGTCAG
GGGCCCCAGTTTATCTTTCAGTATTATGAGAAAGAAGAGAGAGGAAGAGGAAACTTCCCT
GATCAATTCTCAGGTCACCAGTTCCCTAACTATAGCTCTGAGCTGAATGTGAACGCCTTG
TTGCTAGGGGACTCGGCCCTCTATCTCTGTGCCAGCAGCTTGG
  • python -m jcvi.formats.bed uniq：此步将提取的注释中存在重叠的部分删除，最终保留更长的那一个转录本。
  • --primary：此步仅保留每个 transcript id 下最长的 mRNA。该步骤和 uniq 步骤的目的都是为了尽可能去除冗余性以提高后续的共线性得分。

②、对于 CDS，使用 jcvi.formats.fasta 将其 ID 处理为不带版本尾缀的格式：

1
2
python -m jcvi.formats.fasta format --sep=. Pan_troglodytes.Pan_tro_3.0.cds.all.fa.gz chimpanzee.cds
python -m jcvi.formats.fasta format --sep=. Homo_sapiens.GRCh38.cds.all.fa.gz human.cds
1
2
3
4
# 处理前
>ENSPTRT00000098376.1
# 处理后
>ENSPTRT00000098376

进行该步骤的原因主要是 --key 提取的 bed 文件中 transcript_id 里并不携带 version，因此这里通过 --sep=. 去除 CDS ID 中的 version 尾缀使 bed 文件和 cds 文件能够相吻合（在 Ensembl 同一版本的注释信息中不会出现一个 transcript id 同时具有多个 version 的情况）。

③、进行成对共线性搜索：

1
python -m jcvi.compara.catalog ortholog human chimpanzee --cscore=.99 --no_strip_names

该步骤通过 LAST 进行比对，并通过 --cscore=.99 过滤 hit 以保留 reciprocal best hit (RBH)。

请注意 jcvi.compara.catalog ortholog 后紧接两个需要进行分析的物种文件名，此处应与前文输出的文件前缀对应，例如 human 对应 human.bed 及 human.cds。分析后输出的 .anchors 中即包含保留下的高质量共线性块。

1
2
3
4
5
6
7
8
# anchors 文件内容
###
ENST00000616016 ENSPTRT00000090571      2030
ENST00000338591 ENSPTRT00000103781      3000
ENST00000428771 ENSPTRT00000073784      871
ENST00000624697 ENSPTRT00000076196      742
ENST00000453464 ENSPTRT00000098436      1010
ENST00000379339 ENSPTRT00000106256      459

④、筛选 OTO：

.anchors 文件中绝大多数都为 one-to-one 的对应关系，但也依然存在少数为多对一或一对多，因此需要进行进一步的过滤，该步骤可以通过 Python 实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import pandas as pd
import sys

jcvi_result_df = pd.read_csv(sys.argv[1], sep="\t", comment="#", header=None)
jcvi_result_df.columns = ["spe1", "spe2", "score"]

spe1_counts = jcvi_result_df["spe1"].value_counts()
spe1_only_once = spe1_counts[spe1_counts == 1].index.to_list()

spe2_counts = jcvi_result_df["spe2"].value_counts()
spe2_only_once = spe2_counts[spe2_counts == 1].index.to_list()

jcvi_result_filter_df = jcvi_result_df.loc[
    (jcvi_result_df["spe1"].isin(spe1_only_once)) &
    (jcvi_result_df["spe2"].isin(spe2_only_once))
]

jcvi_result_filter_df[["spe1", "spe2"]].to_csv(f"{sys.argv[1]}.one2one", sep="\t", index=False, header=False)

将上述代码保存到一个名为 jcvi_filter.py 的文件中，运行以下命令：

1
python jcvi_filter.py human.chimpanzee.anchors

随后将生成一个名为 human.chimpanzee.anchors.one2one 的文件，其中每一行都为一个 human 和 chimpanzee 间的共线性一对一直系同源基因。

其他需要注意的地方

由于 JCVI 除了考虑比对情况以外，也考虑了这些基因在基因组位置上的线性关系，因此使用 JCVI 有利于回收在两种物种基因组上位置顺序保守的直系同源基因对，但相对的也掩盖了一些在不同物种间位置不一致的直系同源关系。如果你并不关心基因位置信息，可以考虑直接使用 RBH 方法进行确定。