Nature Communications: 使用自然语言处理解密微生物基因功能

Nature Communications: 使用自然语言处理解密微生物基因功能

1. 摘要

通讯作者为：以色列特拉维夫大学的David (Dudu) Burstein

在测序数据量不断增加的时代，揭示未表征基因的功能是一项基本挑战。在这里，我们提出了一个使用自然语言处理（NLP）中采用的深度学习方法来应对这一挑战的概念。我们重新调整了NLP算法的用途，以基于其基因组背景下超过3.6亿个微生物基因的生物语料库为基础，对**“基因语义(gene semantics)”进行建模。我们使用语言模型预测了56, 617个基因的功能类别，发现在1369个与最近发现的防御系统相关的基因中，98%的基因被正确推断。然后，我们系统地评估了不同功能类别的“发现潜力(discovery potential)”**，确定了那些具有最多基因尚未表征的功能类别。最后，我们展示了我们的方法发现与微生物相互作用和防御相关的系统的能力。我们的研究结果强调，将微生物基因组学和语言模型相结合是揭示基因工程的一个很有前途的途径。

2. 引言

在后基因组时代，大量的遗传数据正在迅速积累。特别是，宏基因组学，即直接从其生态系统中对微生物群落进行DNA测序，提供了未开发的数据，包括从未在实验室环境中培养过的大量微生物^[1，2]。人们对这些微生物编码的相当一部分基因的功能知之甚少。破译未表征基因的功能是当今微生物基因组学的一大挑战。这些基因作为基因组操作工具、抗菌药物、递送系统等，可能对生物技术和医学具有巨大价值^[3-5]。
实验和计算研究表明，基因组背景，即位于给定基因附近的一组基因，承载着有关基因功能的重要信息^[6-10]。这种现象在原核生物中很突出，在原核生体中，共功能基因通常在基因组中组织成簇。这种基因组基因座的一个突出例子是CRISPR-Cas系统，它编码一系列赋予对外来遗传元素抗性的基因。虽然不同系统类型的cas基因含量不同，但CRISPR-cas基因座内cas基因亚群的共存是该系统的一个强大的基因组特征^[10-13]。

利用上下文推断意义是自然语言处理领域的一个关键概念。许多应用于英语等自然语言的模型，使用句子中单词的上下文来学习其语义^[14，15]。现代NLP方法在大型文本语料库（如维基百科、新闻文章和其他特定领域的数据源）上训练深度学习算法，为单词提供有意义的数字表示，从而破译其含义和语义关系。这些数字表示被称为“嵌入embeddings”，用于各种下游应用，从主题文本分类到模拟对话的聊天机器人。最近，基于NLP的方法已被应用于对“蛋白质语言”进行建模，即基于氨基酸在属于特定蛋白质家族的序列语料库中的上下文来预测氨基酸的特性。此类应用已被用于模拟各种蛋白质特征^[16-19]，发现抗微生物肽^[20]，甚至预测导致病毒逃逸的抗原^[21]。一种不同的应用旨在使用Pfam结构域而不是氨基酸作为语言模型的输入对生物合成基因簇进行分类^[22]。其他研究将NLP算法应用于DNA k-mers进行分类^[23]，预测增强子-启动子相互作用^[24]和染色质可及性^[25]。

在这里，我们在更高层次的表示上使用了NLP，试图创建一个通用的“基因语义”模型。在我们的模型中，基因家族是由“基因组句子”组成的“单词”。为了生成这些句子，我们重新注释并分析了公开可用的基因组和宏基因组数据集，由超过2.5万亿碱基对的组装序列数据组成。我们将基因数据转换为语料库，通过将基因聚类到家族中添加了一层抽象。我们根据基因家族的基因组背景对其进行建模，以研究其“语义”，并预测数万个未表征基因的功能。我们验证了我们的方法，证明它可以正确地恢复最近发现的系统。最后，我们评估了哪些功能类别具有最高的“发现潜力”，并强调了我们的方法揭示的先前未表征系统的三个例子。

本节提到的重要论文:

16. Elnaggar, A. et al. ProtTrans: towards cracking the language of lifes code through self-supervised deep learning and high performance computing. IEEE Trans. Pattern Anal. Mach. Intell. 14, 7112–7127 (2021).
17. Asgari, E. & Mofrad, M. R. K. Continuous distributed representation of biological sequences for deep proteomics and genomics. PLoS ONE 10, e0141287 (2015).
18. Rives, A. et al. Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences.
    Proc. Natl Acad. Sci. USA 118, e2016239118 (2021).
19. Bepler, T. & Berger, B. Learning the protein language: evolution, structure, and function. Cell Syst. 12, 654–669.e3 (2021).
20. Ma, Y. et al. Identification of antimicrobial peptides from the human gut microbiome using deep learning. Nat. Biotechnol. 40, 921–931.
    https://doi.org/10.1038/s41587-022-01226-0 (2022).
21. Hie, B., Zhong, E. D., Berger, B. & Bryson, B. Learning the language of viral evolution and escape. Science 371, 284–288 (2021).
22. Hannigan, G. D. et al. A deep learning genome-mining strategy for biosynthetic gene cluster prediction. Nucleic Acids Res. 47, e110 (2019).

3. 方法

3.1 数据集汇编和初始基因注释

我们下载了除**mettazoa(后生动物, 除原生动物外所有其他动物的总称)、Fungi和Viridiplantae(植物界)**外的所有基因组，以及NCBI WGS²⁶和EBI Mgnify²⁷中公开的所有宏基因组组装。

总体而言，该数据集包括596,338个基因组和22,923个宏基因组。这些组合代表了巨大的系统发育多样性和各种生态系统，但它们高度偏向于与人类相关的微生物(参见补充图1基因组的分类分布和补充表2所使用样本的细分)。**为了去除基因水平上上下文数据较少的contigs，我们过滤掉了小于10 kbp的contigs。**正如预期的那样，这对宏基因组序列的影响比基因组序列的影响要大得多:只有26.4%的宏基因组序列数据通过了这个阈值(在2.01 Tbp的初始数据集中保留了530 Gbp)，而基因组序列的这一比例为95.1%(在2.25 Tbp的初始数据集中保留了2.14 Tbp)。使用prodigal版本3.0.0⁵⁸进行基因预测，使用Prokka管道⁵⁹(版本1.14.6)进行初始基因注释。为了去除重复的基因组并减少微生物组成的偏差，我们只使用了与Uniprot非冗余蛋白组对应的基因组⁶⁰。该集合由在物种水平上聚集的基因组组成，并经过手动和自动管理以消除冗余⁶¹。对于宏基因组，我们使用BBMAP dedupe utility⁶²(38.69版本)来去除高度相似的contigs。在这些筛选之后，数据集包括在11,119,550个contigs上编码的394,374,454个基因(1100万条contigs上的3.9亿个基因)。

基于KO⁶³进行功能标注。截至2021年5月14日，共从KEGG数据库中下载了24307个KOs中的17,107,806个蛋白。使用参数为-s 7.5 -c 0.5的mmseqs2 cluster⁶⁴进一步对与每个KO相关的蛋白进行亚聚类。每个含有5个以上KEGG蛋白的子簇与MAFFT⁶⁵进行比对，并使用HMMer suite³⁷ version 3.3.2构建HMM图谱。总的来说，KO HMM数据库包括63234个HMM，代表23199个KO。

3.2 语料库生成(Corpus generation)

为了生成基因语料库，使用上面描述的KO HMM数据库扫描由非冗余组编码的蛋白质。与KO HMM显著匹配的蛋白(使用hmmsearch, E-value阈值为10−6)，以KXXXXX的形式根据得分最高的KO子簇进行分类。其中，KXXXXX为KEGG的KO, YY为子簇索引。未映射到KOs的蛋白质被认为是**“假设的蛋白质(hypothetical proteins)”**，并根据氨基酸序列相似性迭代聚类如下。

首先，使用参数为-s 0.80 -c 0.80的CD-HIT⁶⁶对高度相似的假设蛋白进行聚类。其次，使用mmseqs2⁶⁴将CD-HIT代表聚类到基因家族中，其参数与用于对KO进行子聚类的参数相同(-s 7.5 -c 0.5)。每一簇“假设的蛋白质”被分配了一个标识符hypo.clst.ZZZZ，其中ZZZZ是唯一的簇索引。通过测试不同的阈值，以每个家族至少24次出现的阈值过滤掉罕见的基因家族(token)(见补充表3)。出现频率大于10⁻³的无所不在的token也被从语料库中过滤掉，从而得到由563,589个唯一“单词”代表的360,039,110个基因的语料库大小。

我们基于KO数据库进行基因家族注释，因为KO数据库非常广泛，并且包含每个基因家族功能的详细元数据和系统信息。然而，并不是所有的注释蛋白都能被指定KO。为了评估有功能注释但没有KO分配的基因数量，我们使用DIAMOND version 2.0.11⁶⁷在NCBI的非冗余蛋白数据库(NR)中搜索了每个假设集群的代表。如果蛋白质没有在E-value<10⁻⁴下没有命中(hits)，则认为它们未被注释(“假设”)。如果蛋白质具有显著的同源性(E-value<10⁻¹⁰)，则被认为具有可靠的注释(不包括假设蛋白质、假定功能和未知功能域的注释)。E-value在10−4和10−10之间的蛋白质被认为是数据库命中的蛋白质，但该注释被认为不可靠，不能用于下游分析。按照注释程序，每个基因组或宏基因组组被视为一个句子。分配给每个ORF的标识符(基于KO子集群或假设集群)用作句子中的单词(令牌)。

3.3 基因空间生成(gene space generation)

我们的目标是以无监督的方式学习语料库中每个基因家族的数字表示，这里表示为**“嵌入(embedding)”。更正式地说，让我们定义K为所有KEGG标识符的集合，H为所有假设的基因家族标识符的集合。给定一组记号$G=(g_1，\dots ，g_n)$，我们定义了一个函数f，使得$f:g_i\to R^k,k\in N$。基于基因组上下文可以用来暗示基因功能的假设，我们使用了word2vec^29,68中提出的结构和Skip-gram模型**。这个架构是基于单层神经网络,旨在学习嵌入一个简单的任务:预测给定单词的相邻的词,也就是说,最大化对数概率$p(g_{i+j}\mid g_i)$, 其中 $w$是中心词$g_i$周围选择的窗口大小：$1\le i\le \mid K\cup H|,-w\le j\le w\in N$:
$$\frac{1}{\mid K\cup h\mid }\sum _{i=1}^{\mid K\cup h\mid }\sum _{w\le j\le w}logp(g_{i+j}\mid g_i)$$
在实践中，为了加速，我们使用Mikolov等人提出的负采样目标，而不是使用softmax函数来评估$p(g_{i+j}\mid g_i)$。每个$g_i$的最终嵌入可以得到:
$$f(g_i)=W^T\cdot v_{g_i}$$
其中$\mid K\cup H\mid =n$, $W^{n\times k}$是网络第一层的权值矩阵，${\nu }_{g_i}^{n\times 1}$是令牌gi的单热编码向量。我们用窗口大小$w=5$和$k=300$ 维的向量表示来描述每个标记的“基因嵌入空间”。

NLP领域中的token和tokenization到底指的是什么？

• token(符号):包括单词和标点
• tokenization(分词)：我是中国人->[‘我’, ‘是’, ‘中国人’]

3.4 基于嵌入的分类模型

我们利用KO令牌的预训练嵌入作为分类任务中的唯一特征(参见下面应用的交叉验证的详细信息)。我们使用了一个简单的四层深度神经网络，权重为$W_i$，偏置为$b_i$, $i\in [0,1,2,3]$，隐藏层H1,H2,H3的大小分别为256,128和64。我们定义了输入嵌入矩阵$X\in R^{\mid K\mid \times K}$和输出层$Y\in [0,1{]}^{\mid K\mid \times c}$，其中c是我们分类中的类别数量，它为每个预测类别提供了标签分数(见下文选择类别的标准)。

预测的类别定义如下:
$$argmaxsoftmax(H_3^T\ast W_3+b_3)$$
为了避免过拟合，我们为所有隐藏层添加了0.2的dropout。我们使用ReLU作为每个隐藏层和输出层的softmax的激活函数。使用Adam优化器、分类交叉熵损失和20次epoch对网络进行训练。在其他参数固定的情况下，根据损失函数的收敛性选择epoch的个数。我们测试了六种不同的架构:一个或两个隐藏层256、128或64个神经元，保持第n - 1层的大小大于第n层。

我们还测试了SVM³²、随机森林³³和梯度增强分类器³⁴。每个模型的参数都使用嵌套的五倍交叉验证网格搜索对用于模型评估的训练数据进行优化。我们总共训练了336个模型来探索多种参数组合(补充表5)，并根据F1得分选择最佳参数组合。各分类器类型中表现最好的模型为:RBF核、正则化为1 (C = 1)的SVM;随机森林有1000个估计器，最大深度为50;梯度增强分类器有800个估计器，最大深度为6，学习率为0.05。各种模型的训练时间对比见补充图4和补充数据集4。

选择深度神经网络是因为它的性能和效率。为了更好地对其预测进行优先级排序，我们通过加权分类分数和我们的模型在AUPR方面对相关函数的准确性来分配预测分数。作为信度阈值，我们使用t = 0.9作为根据AUPR加权的分数，t’= 0.99为非加权分数。

3.5 功能分类采样

KEGG数据库包括对每个KO标识符的多级功能标注。我们使用每个KO的层次结构的第三级功能类别。对于与多个功能类别相关的KO，根据丰度、与微生物功能的相关性和详细程度选择单一类别(KO及其类别分配表参见补充数据集8)。所有类别的集合 C = { c 1 , c 2 ， … c n } C = \{c_1,c_2，…c_n\} C={c1​,c2​，…cn​}被简化为一个更小的集合 C ∗ = { c 1 ∗ , … , c m ∗ } C^* = \{c^*_1,…,c^*_m\} C∗={c1∗​,…,cm∗​}，使得$c_i^{\ast }\in C$ , $c_i$在语料库中有80多个出现(第90百分位)，其大部分KO属于一个或两个功能类别，其基因具有一定的语境组织。

所选择的类别是氨基糖和核苷酸糖代谢(312个基因)、苯甲酸酯降解(168个基因)、能量代谢(209个基因)、氧化磷酸化(409个基因)、卟啉和叶绿素代谢(274个基因)、原核防御(871个基因)、核糖体(489个基因)、分泌系统(1617个基因)和双组分系统(1174个基因);每个KO的功能类别在补充数据集8中指定。为了解释任何$c_i\in C\setminus C^{\ast }$，我们创建了一个标记为"其他"的类别，从中不包括在$C^{\ast }$的每个类别中随机抽取20个KO，这些KO出现次数超过50次。

3.6 对远程同源方法进行基准测试

我们将该方法的功能分配与三种广泛使用的基于序列的远程同源搜索方法(HHblits, HMMer和PSI-BLAST)进行了比较。对于我们的模型，我们应用了下面描述的**“leave-one-KO-out cross-validation”**过程。对于基于序列的方法，我们针对数据集中未分配与查询相同KO的所有单词搜索每个单词，并根据最佳命中确定预测的功能类别。采用HHsuite3³⁸构建HHsuite数据库，采用HHblits 3.3.0版本，2次迭代，1次迭代E-value阈值为10^-3。为了运行具有两次迭代和参数-N 2 -E 0.001—incE 0.001的jackhmmer，使用了H M M e r³⁷版本3.3.2。对于PSI-BLAST，使用BLAST 2.7.1版本构建BLAST数据库，使用PSI-BLAST³⁶进行2次迭代搜索。对E-value < 10^-3 (value- thrresh参数)的结果进行迭代过滤。如果没有找到匹配项，则预测类别为“未命中”。超参数选择基于对每种方法进行两次、三次或四次迭代测试，e值包含阈值为10^-3、10^-4、或10^-6(见补充图、补充表6和补充数据集4)。

3.7 交叉验证

为了在测试分类器时考虑进化相关性并减少保留集中可能存在的偏差，我们对语料库进行了两个交叉验证程序:(i)朴素的五倍交叉验证和(ii)“leave-one-taxonomy-group-out cross-validation”。在朴素方法中，整个语料库数据被分成五个独立的集，其中20%的组态被忽略。在80%的句子上使用word2vec计算嵌入。在测试分类器时，针对KEGG数据库搜索holdout集中的每个基因以提取相关的预训练嵌入，并将其功能类别与模型的预测进行比较。

对于分类感知交叉验证(taxonomy-aware cross-validation)，我们仅使用基因组数据，其中包括56,736个基因组(在补充数据集2中列出)。为了确保可靠的分类定位，我们排除了宏基因组数据。我们选择了五个最具代表性的分类群:γ变形菌门和α变形菌门，放线菌门，厚壁菌门和拟杆菌门。使用这五组，我们进行了交叉验证(如上所述的初级交叉验证)。

为了将远程同源搜索算法与我们的方法进行比较，我们使用“leave-one-KO-out”交叉验证测试了它们为每个KO正确分配功能的能力。我们在上面描述的交叉验证中使用了同样的7043个单词来进行性能评估，这样对于数据集中所有n = 2915个唯一KO，每个fold包含n - 1个KOs作为训练集,属于其他KO的words作为holdout集。 我们在训练和测试数据集中考虑了属于所选KO的所有词(子簇)。

对于交叉验证程序的性能评估，我们获得了每个功能类别的F1分数，考虑所有类别的加权F1分数，以及每个类别的AUPR。每个类别的AUPR是给定功能类别中所有折叠的微观平均值(micro-average)。

3.8 稀疏分析

为了根据我们的预测评估发现潜力，我们对每个功能类别进行了单独的稀疏分析。

我们定义 P c = { p 1 ， … , p k } P_c = \{p_1，…, p_k\} Pc​={p1​，…,pk​}为功能类别$c$中的预测单词(假设基因家庭)集合。并且， 我们还定义了词集$G_c^p$和词数$W_c^p$,相应地,功能类别$c$的预测词$p$。这个量是指基因的数量为代表的语料库,属于词p所代表的家族。最后, 我们定义 G c = { G c p ∣ p ∈ P c } G_c = \{G^p_ c∣p \in P_c\} Gc​={Gcp​∣p∈Pc​}作为c类预测的所有基因的集合。我们在$n\in [10^3,10^6]$为子样本量，因此我们从$G_c$中均匀抽取n个基因，得到该子样本所代表的独特基因家族数量(记为$P_c^n$)。通过将每个大小为$n$的子样本与其对应的采样基因家族总数$P_c^n$对齐来计算稀疏曲线。置信区间通过10,000个独立样本的自举获得。

3.9 预测的系统候选选择

为了寻找组成未知系统的基因，我们搜索了具有相同预测功能类别的共现基因家族，这些基因家族在语料库中非常频繁，并且倾向于同时出现在多个contigs中。具体来说，我们选择了在语料库中出现次数超过1000次的基因家族，构建了一个二元矩阵${\mathcal{I}}^{n\times m}$, 其中n个contigs, m个selected words，如果contigs有单词j，使得${\mathcal{I}}_{ij}=1$，，则${\mathcal{I}}_{ij}=0$，否则${\mathcal{I}}_{ij}=0$。该矩阵用于编制相关矩阵，从中提取高度相关的集群作为系统候选。使用HHpred³⁸ version 57c8707149031cc9f8edceba362c71a3762bd bf8，使用默认参数对PDB mmCIF70 12Oct、PfamA v35、NCBI conserveddomains (CD) v3.18和TIGRFAMs v15.0数据库进行域和远程同源性探索。使用BLAST69在NCBI的NR数据库中验证所有候选基因，并针对PADLOC⁵³和DefenseFinder⁵⁴测试防御系统基因家族的新颖性。

4. 结果

4.1 基因组语料库汇编

我们的语料库是从NCBI’s²⁶和EBI’s²⁷数据库中公开的所有组装好的宏基因组和基因组(不包括绿色植物、真菌和动物)中编译而成。在去除冗余和短contigs后，该数据集包含1100万个contigs，编码约3.6亿个基因。基于KEGG同源群进行基因家族注释²⁸(图1a)，导致我们数据集中74%的基因被注释。其余基因缺乏明确的KEGG注释，根据序列相似性将其聚类到基因家族中**(图1b)**。每个有足够代表性的基因家族(≥24个基因，见方法)，无论是有注释的还是没有注释的，都被认为是我们基因组语料库中的一个“单词”，从而产生563,589个单词的“基因组词汇”。

值得注意的是，尽管标注的基因在语料库中占大多数，但在将它们聚类后，它们仅占语料库中基因家族或独特“词”的7.8%。**这意味着具有核心微生物功能的已被充分表征的基因聚集在相对较少的、非常大的家族中，而语料库中的大多数遗传多样性(92.2%)没有得到很好的注释。**为了尝试对未被分配KEGG同源的基因进行注释，我们搜索了NCBI的非冗余蛋白数据库(NR)，并恢复了14%的KEGG未注释的基因家族的信息注释(在排除诸如“假设的蛋白质”或“未知功能域”等描述后)。总体而言，约80%的基因家族在两个数据库中都没有信息注释(图1b和补充表3)。

**图1 |模型工作流程。a数据处理和注释。**从公共数据库中下载组装的contigs，并进行基因识别和注释。加注释和未加注释的基因都聚集在基因家族中。**b语料库中标注和假设基因的分布。**左边的条形图代表了语料库中使用的所有约3.6亿个基因。右边的条形图表示约560,000个独特的基因家族(即语料库“词汇”)。**c.英语语料库与基因组语料库的比较。**基因组语料库中的“句子”是由基因家族标识符组成的“词”。**d.嵌入生成和功能预测管道。**嵌入(数字向量表示)由word2vec算法生成，并作为基因功能分类的深度神经网络的输入。

4.2 基因注释嵌入空间

由于大量未注释的基因引起了我们的兴趣，我们试图使用NLP方法更好地理解它们的功能。这种方法依赖于神经网络算法，该算法用于根据文本上下文将单词编码为数字向量。这些向量，或“词嵌入”，旨在封装词的语义，假设一个词是一个词，它的意思是一个词，它的意思是一个词。在这里，我们将学习单词的分布式表示的问题转移到学习基因功能的相同表示(图1c)。也就是说，我们将基因组区域分析为由基因组成的句子而不是单词。在实践中，我们在整个基因语料库上训练了一个简单的无监督神经网络模型word2vec²⁹，并使用该模型学习到的嵌入来创建一个“基因注释空间”(图2a)。**该模型捕获了整个基因组语料库中的基因共现关系，这样具有相似上下文的基因将在基因嵌入空间中相邻。**此外，表示基因对的载体之间类似的代数关系可能意味着它们之间类似的功能相互作用(例如，传感器和调节器对之间类似的距离和方向，见补充说明和补充表1)。

我们首先检查了功能注释基因的嵌入。令人欣慰的是，具有相似功能的基因在基因注释空间中倾向于聚集在一起。例如，几乎所有的cas (crispr相关)基因在嵌入空间中聚集在一起，并且位于其他原核防御系统的附近。虽然CRISPR-Cas簇位于嵌入空间附近，但与其他原核生物防御簇明显可分离**(图2b)**，符合目前对防御岛的理解¹⁰。我们还发现了代表功能群的簇，如分泌系统，它们根据系统类型明显聚类，甚至在一定程度上保留了进化关系: 例如，IV型分泌系统具有与之共享祖先的共轭系统相交的簇³⁰(图2c)。总的来说，我们发现已知功能的基因倾向于在我们创造的高维空间中聚集在一起。

虽然具有相似注释的基因通常在基因空间中聚在一起，但也有一些有趣的偏差。我们发现具有相同注释的基因在嵌入空间中彼此相距很远，这意味着这些基因家族可能在不同的基因组环境中运作。例如，与大多数cas基因不同，并非所有cas4“词”都位于CRISPR-Cas区域内。它们主要分布在基因嵌入空间的四个不同的簇中**(补充图6)**。最大的簇位于CRISPR-Cas区域内，另外两个簇靠近限制性修饰基因，另一个小簇位于我们的嵌入空间中多个未知基因旁边。最近的研究结果表明，cas4家族的基因同时具有crispr相关和非crispr相关的功能，这可以解释这一点³¹。这使得我们的方法能够捕获具有相同注释的基因所执行的不同功能。

**图2 |基因注释嵌入所跨越空间的二维表示。a.全球基因嵌入空间，包括所有563,589个基因家族。**空间中的每个点代表一个基因家族，其中浅红色点代表已注释的基因家族，灰色点代表未注释的基因家族。橙色圆圈标记了包含最多CRISPR-Cas基因的区域(与其他防御基因一起在图b中放大)。**b. 用圆圈标记的防御系统集群区域:cas基因为浅红色，浅蓝色点代表已知的原核防御基因。**红圈集中于上部CRISPR簇，富含I型CRISPR- cas系统基因。**c.按系统类型用颜色编码的包含大多数注释的分泌系统簇的区域。**源数据作为源数据文件提供。

4.3 基于嵌入的功能分类

接下来，我们希望使用注释基因的嵌入来训练基因功能预测模型。我们使用KEGG的功能层次结构用十个功能类别中的一个来标记每个注释的基因(如图3所示)。以基因嵌入为输入，我们训练了四个分类器:支持向量机³²(SVM)、随机森林³³、XGBoos t³⁴和深度神经网络³⁵(DNN)来将基因分配到一个功能类别。我们使用基于分类学的交叉验证来评估分类器的性能。我们将基因组数据集划分为5个主要的分类类群，并迭代测试了基于嵌入和基于其他4个类群训练的classifier准确预测一个主要分类类群中基因功能类别的能力(见方法)。这个过程评估模型在看不见的、进化距离较远的基因组上的表现，我们称之为**“留下一个分类群 leave-one-taxonomic-group-out”交叉验证。在分类性能和速度方面，表现最好的算法是DNN模型(图3a和补充图4)。对于精度召回率曲线(AUPR)下面积为0.56-0.97**的大多数类别，分类器的每个功能类别的性能都很高(图3b和补充图2)。这些结果表明，即使在大进化距离上，基于上下文的功能推断也是有效的。由于分类完全基于基因组背景，没有考虑基因序列或任何其他外部信息。

预测结果表明，嵌入所捕获的功能信息在所有类别中并不统一(图3和补充表4)。例如，嵌入对分泌系统、原核防御系统和氧化磷酸化基因的分类具有很高的信息性。这些功能通常具有很强的环境约束，由多个基因组成，其基因很少属于其他功能途径。相比之下，基因嵌入对功能类别(如氨基糖代谢)的信息较少。通常，缺乏上下文特征和与多种不同途径相关的基因的系统会阻碍分类器的目标(问题：哪些基因缺乏上下文特征或者与多种不同途径相关？占比有多大?)。

鉴于我们的方法具有很高的预测性能，我们希望将它们与基于序列的远程同源搜索算法进行比较。为此，我们将我们的方法与三种成熟的搜索方法(PSI-BLAST³⁶、HMMer³⁷和HHblits³⁸)在KEGG直系同源数据集上进行了基准测试。为了对10个主要功能类别的基因进行交叉验证，我们在KEGG的直系家族(KOs)水平上对属于10个主要功能类别的基因进行交叉验证: 训练了2915个独立的模型，每个模型都遗漏了一个直系家族。在省略“自我命中”(对查询KO的命中)之后，将预测结果与基于序列的搜索结果进行比较。总体而言，在不同的功能类别中，我们的NLP方法通常优于或可与最敏感的基于同源的工具HHblits相比较，尽管仅依赖于嵌入空间中的基因坐标(图3c和补充数据集4)。基于嵌入的分类平均比性能最好的基于同源的方法敏感1.4倍，因为27-44%的测试基因在序列搜索算法中没有显著的命中。我们注意到，与我们的方法相比，某些功能类别的HHblits的精度高出4-6%，假阳性略少，这在依赖序列同源性时是可以预期的。

**图3 |基于嵌入的函数预测性能评估与基准测试。a. 分类器比较。**使用leave-one-taxonomic-group-out交叉验证获得每个类别的F1分数。每个点表示n = 5次折叠的平均值，误差柱为±1 SD。DNN深度神经网络，RF随机森林，SVM支持向量机，XGB XGBoost。b. 使用leave-one-taxonomic-group-out交叉验证计算每个功能类别的Precision-recall曲线(另见补充图2)。“整体”组是指所有类别组合的微平均值(即，将所有类别的预测汇总以计算平均值)。曲线下面积的数值在图中的圆括号中表示为每个功能类别。每条类别线表示交叉验证折叠的微平均值，标准差为±1。**c. 基于leave-one-KO-out交叉验证的方法与远程同源搜索方法的比较。**对用灰点表示的九个功能类别中的每一个都获得了评价指标。柱高为平均值，柱差为±1 SD。源数据作为源数据文件提供。

4.4 基因嵌入可以揭示未标注基因的功能

在对标注基因的模型进行训练后，我们开始为语料库中的假设基因分配功能类别。总的来说，我们的嵌入空间包括519,398个假设的基因家族，其中444,521个在NCBI的NR蛋白数据库中也没有注释。使用我们训练好的DNN，我们预测了56,617个假设基因家族的功能，跨越了总共超过2000万个基因**(图4a, b)。我们只考虑高度可靠的预测，考虑到预测分数和每个功能类别分类的可靠性(见方法)。(问题: 44万个基因家族缺少注释，本文只解决了其中5.6万的基因家族的注释问题，为什么其他的基因家族无法被注释？有什么办法能够对其他的基因家族进行注释?)**

为了进一步验证我们的发现，除了交叉验证分析，我们研究了56,617个基因，我们的分类器提供了高分预测。这些基因未被KEGG注释，但对于7691个基因，我们能够从针对NCBI NR的BLAST搜索中恢复注释**(图4c)。总的来说，基因描述与模型的预测非常吻合。例如，通过对AUPR最低的两个类别的检查，发现预测的“能量代谢”基因确实与KEGG层次结构中的这一类别及其子类别相关，其中包括大量与氮和甲烷代谢相关的基因。同样，在“氨基糖和核苷酸糖代谢”类别中，预测描述最多的基因是糖基转移酶和氨基葡萄糖-6-磷酸脱氨酶。在补充数据集7**中详细介绍了在每个功能类别中预测的基因恢复的完整NR注释。

许多原核防御系统最近被发现和实验验证。由于与这些系统相关的基因在KEGG中尚未被注释，并且在NR中大部分未被注释，因此我们可以利用这一庞大的信息来进一步验证我们的方法。因此，我们测试了我们检测与BREX、disarmament、Theoris、Durantia、Gabija、Hachiman、Kiwa、Lamssu、Septu、Wadjet和Zoria抗噬菌体系统相关基因的能力^8,39,40。我们的语料库包含1369个与这些系统相关的基因。令人欣慰的是，我们正确地将其中98.6%分类为原核防御基因**(图4c, d和补充表8)**。这表明我们的方法能够检测基因功能，即使它们与用于训练我们分类模型的注释基因没有同源性。

值得注意的是，除了上述系统的基因外，预计还有40,247个(即35292+4955=40247)未表征的基因家族与防御系统有关，这表明未发现的原核防御系统存在巨大的多样性。检查分类器的结果显示，在预测的功能类别中，基因家族的数量变化很大。为了估计每个类别中其他基因的“发现潜力”，我们对每个类别进行了调整后的稀疏分析**(图4e)**。我们反复对数据集中的基因进行亚采样，并测试在不同样本量的每个功能类别中预测了多少个基因家族。这使我们能够推断出在每个功能类别中仍有多少基因家族有望被发现。正如预期的那样，在与核心细胞功能相关的已得到充分研究的类别中，如核糖体基因和核苷酸代谢，增加基因数量并没有发现额外的基因家族(在稀疏图中由平台表示)。然而，原核防御系统、分泌系统和双组分系统的功能分类显示出很高的发现潜力，这一点从这些分类中预测的基因家族数量随着分析中基因总数的增加而不断增加可以看出。事实上，在这项研究的结论之后，两项大规模的研究发现了42种原核防御系统^41,42，证实了这一类别的巨大发现潜力。稀薄分析提供了一个系统的和定量的评估“可发现性”，但尚未表征的系统在不同的功能类别。

**图4 |假设基因家族的功能预测。a.完整预测空间。**基因家族根据预测的功能类别进行颜色编码。**b.每个功能类别的预测。**每个条形图代表一个类别中假设基因的总数。黑点表示接受分类预测的假设基因家族的数量。每个点旁边都明确地说明了预测的家庭数量。条形图由用于训练模型的每个功能类别的单词数量进行颜色编码。**c.在每个类别中获得可靠预测的基因数量，在NCBI NR中分为有和没有信息注释的基因。d.新近发现的防御系统中基因家族的功能预测。e稀薄分析。**每一行对应一个功能类别。x轴表示采样基因的数量，y轴表示子样本中具有预测功能类别的基因家族的数量。源数据作为源数据文件提供。

4.5 基于嵌入的预测识别未表征的细菌膜机制

接下来，我们不仅使用我们的预测来检测分离的基因功能，而且还使用具有类似预测功能的共同发生的假设基因簇。这些可以揭示以前未表征的系统或已知系统的扩展功能。我们最初专注于与分泌系统相关的基因，我们的**稀疏分析(rarefaction analysis)**表明这是一个具有很高发现潜力的类别。我们寻求可靠的预测，定义为那些在众多基因组中具有显著共现信号的预测。我们确定了两个这样的候选者:一个假设的分泌相关系统和一组与IV型菌毛系统高度相关的基因。

第一个与分泌有关的系统是一个由8到9个基因组成的操纵子，这些基因都被我们的分类器标记为分泌系统基因。这该系统普遍存在于梭菌纲的3个属:蔷薇菌Roseburia、反刍菌Ruminococcus、和真杆菌Eubacterium(检出总数超过2000次)。这些属的成员是人类肠道微生物群中丰富的革兰氏阳性菌，影响健康和疾病的多种代谢途径^{43 - 45}。我们观察到两种主要变异，一种存在于Roseburia，另一种存在于Ruminococcus和Eubacterium。这两个变异都有一个菌毛组合atp酶(CpaF, PilB/GspE家族成员)，一个IV型菌毛/ II型分泌系统膜蛋白(PilC/GspF)，一个具有未鉴定的DUF5411结构域的蛋白，以及一个具有细胞侵袭结构域的蛋白(与Internalin-B具有远程和部分同源性，图5a，补充图7和补充表9)。所有其他基因与任何已知的IV型菌毛/ II型分泌成分没有显着的相似性，但在某些情况下，位于编码PilA或粘附蛋白的基因附近。

第二个推测的分泌相关操纵子编码于细孔菌属(即Veillonella genus )，与IV型菌毛系统密切相关，与II型分泌系统具有密切的同源性。该系统在NCBI的WGS数据库(1280个基因组)中86%的细孔菌基因组上被发现。细孔菌是在哺乳动物口腔黏膜和肠道中发现的厚壁菌门的革兰氏阴性细菌，以其与其他生物体相互作用的能力而闻名，特别是在生物膜中，通过共聚集/共粘附。大多数细络菌与常见口腔疾病的发生有关^46,47。我们的分类器识别出五到六个与IV菌毛蛋白相邻的未注释基因**(图5b，补充图7和补充表10)**。这似乎是IV型菌毛的远缘变异:虽然有5个基因编码已被充分表征的Pil蛋白，但我们的分类器鉴定的基因未被注释，仅与IV型菌毛具有远缘同源性结构域，主要指向内膜辅助蛋白。值得注意的是，我们检测到所有成分的同源性，除了外膜PilQ/GspD蛋白，该蛋白与该基因组位点有数百个基因距离。细毛杆菌属属革兰氏阳性厚壁菌门中的革兰氏阴性菌，这种IV型菌毛系统的组织和表达可能反映了对其独特外膜的适应⁴⁸。

4.6 基于嵌入的分类器揭示了一个假定的原核防御系统

基于嵌入的分类器揭示了一个假定的原核防御系统。我们的分析表明，在分析的功能类别中，原核防御系统具有最大的发现潜力。因此，我们使用我们的分类来识别共同发生的非特征基因家族，这些基因家族被预测具有原核防御功能。这使我们发现了一个假定的防御系统，其中包含3到5个未表征的基因(图5c，补充图8和补充表11)。如下所述，该系统编码各种DNA结合和切割结构域，并且广泛存在于整个细菌界(图5d)。

我们确定了两种类型的这种假定的系统，在他们的基因含量不同，其中大部分是未表征的，尽管少数包含注释域。该系统的核心基因由所有变体共享，编码:(1)925 aa的大蛋白，具有单个未表征的z1结构域，据报道与限制性修饰系统相关⁴⁹;(2)具有PD-(D/E)XK基序和DUF4220结构域的蛋白，与核酸酶活性和II型限制性内切酶相关⁵⁰;(3)胞嘧啶甲基转移酶。

系统的I型(图5c)由6个基因组成: 3个核心基因和3个附加基因。这种类型与已知的DNA修复内切酶密切相关。II型分为II- a和II- b两个亚型。这两种亚型都包括该系统的core蛋白和一种额外的AIPR家族蛋白，这是一种流产感染蛋白⁵¹。这两种亚型与MORC家族cw型锌指蛋白和一种推定的转录调节因子高度相关(图5c)。II型让人想起最近发现的mzaABCDE系统⁵²，与mzaBCDE组件显示低序列相似性(~20%的同一性)，并且mzaA (sigma因子)仅存在于少数系统中(我们鉴定的系统的新颖性也通过PADLOC⁵³和DefenseFinder⁵⁴进行了验证)。需要进一步调查以确定该系统的功能和机制。然而，大量DNA结合和切割结构域的存在，以及与已验证的mzaABCDE系统共享的结构域，表明我们基于嵌入的分类器确实揭示了一个原核防御系统。这种操纵子结构在两种系统类型的广泛分类分布中的存在进一步证实了这可能是一个真正的防御系统。

**图5 |预测系统。**对于所有的基因操纵子，已知的注释在基因图的下方表示，未注释基因的结构域在基因的上方用箭头标记。a.在3个梭菌属中丰富的预测的与分泌有关的操纵子。通过我们的方法预测的基因被标记为黄色/橙色渐变着色。b.在细枝菌属的两个代表中，IV型菌毛系统的一个遥远变种。用我们的方法预测的基因用深浅不一的蓝色表示。pil基因为带注释的IV型菌毛基因。c.在多种细菌中发现的候选防御系统，具有来自四个基因组的代表。d. c中所示的系统的细菌分布。图中最上面的图包括细菌生命树70，用每种系统类型的存在进行颜色编码。下面的面板说明了每个系统在顺序级别上的分类分布。源数据作为源数据文件提供。

5. 讨论

我们提出了一种方法，旨在根据其基因组背景普遍表征微生物基因之间的功能关系。这些关系由“基因空间”捕获，即根据基因家族共现性计算的数学表示(嵌入)。根据这种方法的假设，我们发现大多数具有相似功能的基因倾向于在基因嵌入空间中聚集在一起。此外，在一些有趣的情况下，在嵌入空间的不同区域检测到具有不同功能的同一基因的变体。我们利用基因嵌入作为深度学习分类器的输入来预测基因功能，并发现它表现良好，具有一定的差异能力，可以更好地推断某些功能类别。这可能是由于更强的上下文信号和更少的基因共享在不同的系统/途径。我们的预测突出了具有高发现潜力的功能类别，并显示了未被探索的细菌膜结合机制在人类微生物组和微生物防御系统中运作，这在细菌王国中广泛存在。

NLP方法以前被用于表征基因和蛋白质，主要是在氨基酸或核苷酸k-mers^{16,17,20,22,23}的水平上。在这里，我们瞄准了一个更粗糙的生物学水平，将完整的基因表示为单词，并探索它们的普遍语义。我们分析了一个广泛的基因组语料库，该语料库基于NCBI’s和EBI’s的基因组和宏基因组数据库中所有组装的微生物组群。我们的方法依赖于对语言模型的直观适应，基因组是“句子”，基因是“单词”。将原始序列封装到基因家族中可以减少噪声，并将抽象级别提高到与基因功能建模相关的级别。

为了处理大量的基因组数据并确保有意义的上下文，我们过滤掉了不常见的基因(少于24个代表的家族)和短contigs (<10 kbp)。此外，我们使用了标准的基因预测，这意味着我们的分析不包括小ORFs55。这限制了我们的分析，因为稀有基因和短肽被忽略了，而短contigs的排除可能会阻碍小的移动遗传元件(例如，转座子和小质粒)上的基因检测。但是，可以针对旨在包含这些元素的应用程序调整相关参数。对于基因嵌入空间的计算，我们采用了word2vec算法，该算法相对简单、快速、直观。使用更先进的结构，如变压器(transformer)和长短期记忆递归神经网络(LSTM)可以帮助进一步改善嵌入，并可能更好地解决遗传数据的稀疏性，其中包括许多“未知单词”。这些模型尚未针对极长的基因组序列进行优化，并且需要大量的计算资源，但它们有望用于未来的基因语义模型。【当前研究的局限性：过滤掉不常见基因，短的contigs, 模型上可以进一步优化。】

我们的稀疏分析表明，不同功能类别的发现潜力是高度可变的。“饱和”主要出现在属于大多数微生物“核心基因组”的特征明确的类别中，如能量代谢和翻译。预计这些功能包含相对较少的未表征基因。然而，这种明显的“饱和”也可能是由于上面讨论的方法限制(即缺乏不同系统之间共享的上下文信号或基因)。值得注意的是，作为所谓“云”或“辅助”基因组一部分的功能类别似乎远未饱和。这些辅助基因，其中一部分与防御和分泌有关，是微生物之间以及与宿主、捕食者和环境相互作用的关键。因此，这些基因中的许多都处于进化军备竞赛的关系中，导致它们迅速多样化，这可以解释未知系统的高发现潜力。值得注意的是，辅助基因，特别是那些与原核生物防御和分泌有关的基因，经常参与DNA操纵和毒力^10,42,56,57，因此，可能对生物技术工具和临床应用的发展做出巨大贡献。

这里介绍的实现仅仅触及了使用NLP方法“读取”基因组的潜力的表面。我们使用嵌入技术将基因分类为一组预定义的、一般的、功能的类别。然而，专注于领域特定的功能可以用于深入研究特定的功能，而不是将其用于感兴趣的功能。这可以通过使用这里提供的嵌入来实现，并进行一些微调，作为在特定感兴趣的基因上训练的分类器的输入，或者通过创建基于相关语料库的从头嵌入来实现，例如病毒基因组，特定微生物组，或包含额外单词集(例如，短肽或非编码rna)的词汇表。功能分类器还可以将基因嵌入与基于序列和结构的特征相结合，从而创建考虑感兴趣基因的内容和上下文的模型。

我们提出的方法有能力显著丰富我们的微生物基因功能的知识。它特别适合于推断与特征蛋白没有序列相似性的基因的功能。它可以检测相似基因，在不共享序列相似性的情况下执行相似的功能，并突出同源基因的不同功能，揭示以前未知的基因专门化。描述“基因组语言”的NLP模型可以进一步丰富:未来的模型可以包括，除了共现，还包括关于基因之间的共向性和距离的信息，以及额外的元素，如启动子、终止子和调控元件作为基因组句子的“标点符号”。我们预计，这里提出的结果，加上改进和丰富的基因语义模型，将导致更好地理解基因功能和进化在广大的微生物宇宙。