scGPT方法解读

scGPT

方法概述

单细胞测序能够在单个细胞水平上分析分子特征。例如，scRNA-seq 可以测量RNA转录本的丰度，有助于深入了解细胞身份、发育阶段和功能。scGPT是单细胞领域采用生成式预训练方法的基础模型。其核心模型包含堆叠的带有多头注意力机制的Transformer层，可同时生成细胞和基因嵌入。

scGPT包含两个训练阶段：首先在大型细胞图谱上进行通用预训练，然后在较小的数据集上进行特定应用的后续微调（图1a - c）。在预训练阶段，引入专门设计的注意力掩码和生成式训练流程，以自监督的方式训练scGPT，共同优化细胞和基因的表征。这项技术解决了基因表达的非顺序性问题，使其能够适应顺序预测的自然语言生成（NLG）框架。在训练过程中，模型逐渐学会根据细胞状态或基因表达线索生成细胞的完整基因表达。在微调阶段，预训练模型可以适配新的数据集和特定任务。scGPT提供了灵活的微调流程，适用于单细胞研究中的各种关键任务，包括带批次校正的scRNA-seq数据整合、细胞类型注释、多组学整合、扰动预测和基因调控网络（GRN）推断。

为了收集丰富多样的测序数据用于scGPT的自监督预训练，scGPT从CELLxGENE数据库（https://cellxgene.cziscience.com/，图1d）中收集了3300万个处于正常（非疾病）状态的人类细胞的单细胞RNA测序（scRNA-seq）数据。这个综合数据集涵盖了来自51种器官或组织、441项研究中的多种细胞类型，充分展现了人体细胞的异质性。预训练完成后，用UMAP可视化方法 ，对3300万个细胞中10%的人类细胞的scGPT细胞嵌入进行可视化（图1e）。得到的UMAP图清晰得令人瞩目，不同颜色准确地代表了局部区域和簇中的各种细胞类型。鉴于数据集中纳入了400多项研究，这显示出预训练在提取生物学差异方面的卓越能力。

输入embedding

单细胞测序数据被处理为cell-gene矩阵，$X\in R^{N\times G}$，每个元素$X_{i,j}\in R^{+}$表示分子的丰度，$i$表示细胞$i$。在接下来的章节中，把这个矩阵称为原始计数矩阵。scGPT的输入主要由三个部分组成：（1）基因（或峰值）token；（2）表达值；（3）条件token。对于每一个建模任务，基因token和表达值会根据原始计数矩阵$X$进行相应的预处理。

基因token
在scGPT框架中，每个基因都被视为最小的信息单元，类似于自然语言处理（NLG）中的一个单词。因此，我们将基因名称用作token，并为每个基因$g_j$分配一个唯一的整数标识符$id(g_j)$。这些标识符构成了scGPT中使用的token词汇表。这种方法具有很大的灵活性，能够整合来自不同基因集（即由不同测序技术或预处理流程生成）的多项研究数据。具体而言，通过取所有研究中基因的并集，可以将不同的基因token集整合到一个通用词汇表中。此外，我们在词汇表中加入了特殊标记，比如用$<cls>$将所有基因聚合为一个细胞表征，用$<pad>$将输入填充到固定长度。从概念上讲，我们将基因标记与自然语言处理中的单词标记进行类比。因此，每个细胞$i$的输入基因标记由一个向量$t_g^{(i)}$表示（$M$是预定义的最大输入长度）。

表达值
基因表达矩阵$X$在作为建模输入之前需要额外处理。基因表达建模中的一个基本挑战是不同测序协议下绝对数值的可变性。测序深度的差异以及稀疏表达基因的存在，导致不同批次测序样本的数据尺度存在显著差异。这些差异很难通过常见的预处理技术（如每百万转录本归一化和log1p转换）来缓解。即使经过这些转换，相同的绝对值在不同测序批次中可能传达不同的 “语义” 含义。为了解决这种尺度差异，我们采用了值分箱技术，将所有表达计数转换为相对值。对于每个细胞中的每个非零表达计数，我们计算原始绝对值，并将它们划分为$B$个连续区间$[b_k,b_{k+1}]$，其中$k\in [1,2,\dots ,B]$。需要注意的是，每个细胞都要计算一组新的分箱边界。

比如，值$x_j^{(i)}=B$始终表示基因中最高的表达水平。值得注意的是，对于微调任务，在值分箱步骤之前还进行了log1p转换和高变基因（HVG）选择。为简化表示，混合使用$X_{i,j}$来表示分箱前的原始数据矩阵和预处理后的数据矩阵。因此，细胞$i$的分箱表达值的最终输入向量表示为：

条件token
条件标记涵盖了与单个基因相关的各种元信息，例如扰动实验的变化（由扰动标记指示）。为了表示按位置排列的条件标记，将使用一个与输入基因维度相同的输入向量。这个向量表示为：

其中，$t_{c,j}^{(i)}$表示一个条件的整数index。

Embedding层
分别使用常规的嵌入层（即PyTorch嵌入层，https://pytorch.org/docs/stable/generated/torch.nn.Embedding.html ）$em{b}_g$和$em{b}_c$对基因标记和条件标记进行处理，以将每个标记映射到维度为$D$的固定长度嵌入向量。对于分箱后的表达值，使用全连接层（记为$em{b}_x$）来增强其表达能力。细胞$i$的最终嵌入$h^{(i)}\in R^{M\times D}$定义为：

细胞和基因表达建模

scGPT transformer
scGPT使用自注意力Transformer架构对输入嵌入$h^{(i)}$进行编码。自注意力机制作用于由$M$个嵌入向量组成的序列，这使得它特别适合捕捉基因之间的相互作用。堆叠的Transformer模块的输出可以定义如下：

我们将得到的表征$h_n^{(i)}$用于基因层面和细胞层面的任务。对于细胞任务，首先聚合$h_n^{(i)}$为细胞embedding，然后再进行后续操作。

输入维度$M$所涉及的基因数量可达数万，这大大超过了自然语言处理（NLG）中常用的传统Transformer的输入长度。为应对这一挑战并确保自注意力机制高效运行，采用了FlashAttention所实现的加速自注意力机制。这一实现方式有效提升了模型容量，使模型能够对大输入维度进行有效处理。尽管此处采用了FlashAttention，但原则上任何高效的Transformer架构都可用于scGPT，比如具有线性复杂度的Transformer（Linformer）和核化自注意力机制（KSA）。

细胞表征
每个细胞类似于一个由基因组成的 “句子”，其表征$h_c^{(i)}\in R^D$是通过聚合学习到的基因层面的表征$h_n^{(i)}$获得的。**在这种情况下，各种池化操作，如逐元素平均池化或加权池化，都可以轻易应用。**在scGPT中，选择使用特殊标记$<cls>$来表示细胞，使模型能够在Transformer模块中学习池化操作。$<cls>$标记被添加到输入标记的开头，该位置的最终嵌入被提取出来作为细胞的表征。因此，细胞embedding $h_c^{(i)}$可以从堆叠的最终层嵌入$h_n^{(i)}[<cls>]$提取，其中$[<cls>]$操作是在$<cls>$ token位置的索引处检索。

批次和模态的表征
使用额外的标记集来表示不同的测序批次和测序模态（例如，来自RNA测序的基因、来自ATAC测序的峰值等），这在单细胞RNA测序（scRNA-seq）和单细胞多组学整合任务中尤为重要。这与“输入嵌入”部分介绍的条件标记类似，并且同样通过标准嵌入层来实现。模态标记$t_m^{(i)}$ 与单个输入特征$g_j$相关联（例如，用于指示它是基因、区域还是蛋白质）。批次标记原本作用于细胞层面，但也可以传播到单个细胞的所有特征上。换句话说，相同的批次标记$t_b^{(i)}$ 可以在单细胞$i$的长度为$M$的输入特征中重复出现：

在“输入embedding”中描述的标记与批次和模态标记之间的区别在于：这些批次和模态标记的嵌入并不作为Transformer模块的输入。相反，在进入特定的微调目标之前，它们会在特征或细胞层面与Transformer的输出进行连接。

• 这是为了防止Transformer在放大相同模态特征之间的注意力的同时，低估不同模态特征之间的注意力。此外，了解模态和/或批次身份有助于在下游微调目标中进行基因表达建模。随着模型学习基于模态和/或批次身份来预测表达值，这种偏差会从基因和细胞表征本身中被隐含地消除。这是一种有助于批次校正的技术。

例如，在单细胞多组学整合任务中，scGPT将Transformer的输出与批次嵌入和模态嵌入之和进行连接。这将作为下游用于表达建模的微调目标的输入：

其中，$em{b}_b$和$em{b}_m$表示批次和模态embedding层，$h_n^{(i)}$表示scGPT transformer 层$n$的输出。

或者，在单细胞RNA测序（scRNA-seq）整合任务中，将批次嵌入与细胞表征连接起来，得到以下表示作为输入，这可以用于直接的整合优化：

其中，$t_b^{(i)}$表示细胞$i$的batch ID，$h_c^{(i)}$是细胞embedding。注意，批次整合与多组学整合都是下游微调才使用的目标

生成式预训练

基础模型预训练
基础模型旨在成为一种可通用的特征提取器，能够为各种下游任务带来帮助。预训练中使用的 token 词汇表包含人类基因组中的全套基因。在模型预训练前，对表达值进行了分箱处理（见“输入embedding”部分）。为加快训练速度，将输入限制为每个输入细胞中具有非零表达的基因。为有效训练模型以捕捉基因 - 基因关系和基因 - 细胞关系，scGPT引入了一种带有专门注意力掩码的生成式训练策略。

注意力掩码策略
自注意力机制已被广泛用于捕捉token之间的共现模式。在自然语言处理中，主要通过两种方式实现：

• （1）在诸如BERT和RoBERTa等Transformer编码器模型中使用的掩码标记预测，即在模型输出中预测输入序列中随机掩码的标记；
• （2）在诸如OpenAI GPT系列等因果Transformer解码器模型中进行具有顺序预测的自回归生成。OpenAI GPT-3和GPT-4中使用的生成式预训练采用统一框架，模型从由已知输入token组成的“提示”中预测最可能的下一个token。这个框架在各种自然语言处理应用中具有很大的灵活性，并且在零样本和微调设置中展现出上下文感知等能力。

scGPT认为，生成式训练可以以类似的方式使单细胞模型受益。具体而言，关注两项任务：

• （1）基于已知基因表达生成未知基因表达值，即通过“基因提示”进行生成；
• （2）在给定输入细胞类型条件下生成全基因组表达，即通过“细胞提示”进行生成。

尽管标记和提示的概念相似，但由于数据的非顺序性，对基因读数进行建模与自然语言处理有着本质的不同。与句子中的单词不同，细胞内基因的顺序是可以互换的，并且不存在可预测“下一个基因”这样的对等概念。这使得将GPT模型中的因果掩码公式直接应用于单细胞数据变得具有挑战性。为应对这一挑战，scGPT开发了一种专门的注意力掩码机制，该机制基于注意力分数来定义预测顺序。

注意力掩码通常可以应用于Transformer模块中的自注意力图：对于由$M$个基因标记组成的输入，第$(l+1)$个Transformer模块对其输入的$M$个标记$h_l^{(i)}\in R^{M\times D}$ 应用多头自注意力机制。具体而言，每个自注意力操作的计算方式如下：

其中，$Q,K,V\in R^{M\times d}$为query，key，value。$W_q,W_k,W_v\in R^{D\times d}$是可学习的权重。注意力掩码矩阵 A m a s k ∈ { 0 , − i n f } M × M A_{mask}\in\left\{0,-inf\right\}^{M\times M} Amask​∈{0,−inf}M×M 通过像在 $\frac{Q{K}^T}{\sqrt{d}}$ 中那样修改查询（queries）和键（keys）之间的原始注意力权重，来划定自注意力的作用范围。具体来说，在矩阵中的位置$(i,j)$处加上$-inf$，会使经过softmax操作后的注意力权重变为零，从而阻止第$i$个查询和第$j$个键之间的注意力计算。另一方面，加上 0 则意味着注意力权重保持不变。这种注意力掩码技术能让模型专注于特定的上下文元素。

scGPT的注意力mask经过特殊设计，可以用统一的方式支持gene提示和细胞提示生成。注意力掩码$A_{mask}$如补充图1a所示，其中查询按行排列，键按列排列。

如图底部所标注的那样，输入嵌入$h_l^{(i)}$中的每个标记可以属于以下三组之一：（1）用于细胞嵌入的保留$<cls>$标记；（2）带有标记嵌入和表达值嵌入的已知基因；（3）其表达值有待预测的未知基因。scGPT注意力掩码的经验法则是只允许在“已知基因”的嵌入与查询基因自身的嵌入之间进行注意力计算。这是通过使用掩码矩阵$A_{mask}$中的元素$a_{i,j}$ 来实现的，具体如下：

在每次生成迭代中，scGPT会预测一组新基因的基因表达值，而这些基因反过来会在下一次迭代中成为用于注意力计算的“已知基因”。这种方法通过在非顺序的单细胞数据中进行顺序预测，体现了传统Transformer解码器中带有下一个标记预测的因果掩码设计理念。

如补充图1a所示，在训练过程中，scGPT会随机选取一定比例的基因作为未知基因，这样在输入中会省略它们的表达值。注意力仅应用于已知基因与作为查询的未知基因本身之间，而不会应用于其他未知基因的位置。例如，位于位置$j$的待预测基因仅与细胞嵌入、已知基因及其自身有注意力分数，而与其他未知基因没有，如注意力掩码的最后一行所示。scGPT模型通过上述带有掩码注意力图的堆叠Transformer模块来预测这些未知基因的表达。

推理步骤在补充图1b中展示。在进行细胞提示生成的推理时，scGPT会基于特定的细胞类型生成全基因组范围内的基因表达。在第一个位置输入经过训练的细胞嵌入，以代表细胞类型条件。数千个基因表达值的整个生成过程分K个迭代步骤进行（即，在补充图1b中K = 3步）。例如，在某一次迭代 $i\in (1,2,\dots ,K)$中，注意力掩码机制允许与之前从第0到第 $i-1$ 次迭代中所有已预测的基因产生注意力关联。在每次迭代中，scGPT会从未知基因集中选择预测置信度最高的前$1/K$个基因，将它们作为已知基因纳入下一次迭代 $i+1$中。直观地说，这个工作流程以自回归的方式简化了基因表达的生成过程，即首先生成预测置信度最高的基因表达值，并利用这些值来辅助后续轮次的生成。基因提示生成的工作方式类似，也是以迭代的方式进行。不同之处在于，它是从一组具有观测到表达值的已知基因开始，而不是从细胞嵌入开始。

• 补充图1：scGPT注意力掩码。被掩码的位置用蓝色表示，未被掩码的位置用白色表示。这些被掩码和未被掩码的位置对应于$M$个输入标记的$M\times M$注意力图。行对应于查询Q，列对应于键K。在Transformer的自注意力计算中，被掩码位置上的注意力分数将被去除。
• 与每一列相关联的标记标识标注在下方，即“cell emb$<cls>$”表示细胞嵌入，“genes & expression”表示已知基因，“genes to predict”表示未知基因。
• （A）训练中的scGPT注意力掩码，其中只有查询基因和已知基因参与注意力计算。（B）训练之后，在scGPT细胞提示生成的迭代过程中每一步的注意力掩码。

预训练的学习目标
scGPT采用了一个基因表达预测目标来优化模型，使其能够预测未知基因的表达值。具体而言，使用一个多层感知器网络（MLP）来估计未知的表达值，并计算均方损失$L$：

其中，$U_{unk}$ 表示未知基因的输出位置的集合，$x_j^{(i)}$ 是待预测的实际基因表达值。$|\cdot |$ 运算用于获取集合中元素的数量。

如“注意力掩码策略”中所述，基因提示生成和细胞提示生成这两种方式均被支持。在训练过程中，这两种模式是连续进行的。对于给定的一个细胞的输入基因标记，会选取一定比例的基因作为“未知”基因，并且它们的表达值会被省略。首先，在基因提示步骤中，模型的输入包含$<cls>$标记嵌入、已知基因嵌入和未知基因嵌入。使用上面模型的输出计算损失$L$。其次，在细胞提示步骤中，前一步输出的细胞嵌入（即“细胞表征”$h_c^{(i)}$ ）会被用来替换$<cls>$位置处的嵌入。其他计算保持不变。最后，将这两个步骤的损失值相加，然后用于计算梯度，以优化模型参数。

微调目标

scGPT利用多种微调目标，既有助于学习细胞和基因在生物学上有效的表征，也用于诸如批次校正等正则化目的。

基因表达预测
为了促进对基因间相互作用的学习，scGPT纳入了基因表达预测（GEP）。这一微调目标的工作方式与预训练中的目标类似，但它应用于被掩码的位置。具体来说，对于每个输入细胞，一部分基因标记及其相应的表达值$x^{(i)}$会被随机掩码。对scGPT进行优化，使其能够准确预测被掩码位置处的表达值。这一微调目标有助于模型有效地对数据集中基因之间的共表达进行编码。该目标使被掩码位置处的均方误差最小化，用 $M_{mask}$ 表示。基因表达预测（GEP）的工作原理如下：

其中，${\tilde{x}}^{(i)}\in N^M$表示估计的细胞$i$的基因表达量。注意，如果批次和模态也提供了，使用$h_n^{{}^{\prime }(i)}$代替$h_n^{(i)}$。

基因表达预测（GEP）提出了一个通用的自监督微调目标，旨在预测基因表达值。在某些下游任务中，比如扰动预测，模型需要预测的是受到扰动后的基因表达值，而非原始值。scGPT将这种变体称为扰动基因表达预测（perturb-GEP）。保留上面公式中的多层感知器（MLP）估计器，但使用扰动后的基因表达值作为目标值$x_j^{(i)}$。在扰动基因表达预测（perturb-GEP）中，模型应预测所有输入基因在受到扰动后的表达值。

基因表达预测用于细胞建模
这个微调目标的运作方式与基因表达预测（GEP）类似，但它是基于细胞表征$h_c^{(i)}$来预测基因表达值，从而明确地助力于细胞表征的学习。对于输入细胞$i$中的每个基因$j$，创建一个查询向量$q_j$，并使用$q_j$与细胞表征$h_c^{(i)}$的参数化内积作为预测的表达值：

用于细胞建模的基因表达预测（GEPC）继承了基因标记嵌入 $em{b}_g(t_g^{(i)})$。在整合任务中，使用$h_c^{{}^{\prime }(i)}$代替$h_c^{(i)}$。在实验中观察到，与单独使用基因表达预测（GEP）或用于细胞建模的基因表达预测（GEPC）这两种方法中的任何一种相比，将它们结合起来会显著提升性能。

弹性细胞相似性
这个微调目标通过利用一种相似性学习损失来增强细胞表征：

其中，$sim$表示余弦相似度函数，而$i$和$i^{\prime }$ 指的是mini-batch中的两个细胞。此外，$\beta$表示一个预定义的阈值，$ECS$是弹性细胞相似度（Elastic cell similarity）。这种方法背后的基本思想是增强那些余弦相似度值高于$\beta$的细胞对之间的相似度，使它们变得更加相似。相反，不相似的细胞对则会被促使变得更加不同。

通过反向传播进行领域自适应
细胞表征学习会受到批次效应的阻碍，这种效应源于测序技术引入的非生物学批次差异。为缓解这一问题，scGPT采用一个独立的多层感知器（MLP）分类器，依据每个输入细胞的细胞表征 $h_c^{(i)}$ 来预测其对应的测序批次，并通过在模型内反转梯度来修改反向传播过程。此方法借鉴了加宁（Ganin）和伦皮茨基（Lempitsky）提出的稳健领域自适应方法的思路。

细胞类型分类
这个微调目标旨在利用学习得到的细胞表征来对单细胞进行注释。scGPT使用一个独立的多层感知器（MLP）分类器，根据细胞的表征 $h_c^{(i)}$ 来预测细胞类型。这个微调目标通过预测的细胞类型概率与真实标签之间的交叉熵损失 $ce$ 来进行优化。

在下游任务上的微调

细胞类型注释
对于细胞类型注释任务，scGPT在带有真实标签的参考数据集上对模型进行微调，并在预留的查询数据集上验证注释性能。预训练基础模型与参考数据集之间共有的基因标记集得以保留。在模型微调之前，基因表达值先进行了归一化、对数转换和分箱处理。除了输出细胞类型分类器是随机初始化的之外，所有预训练模型的权重都用于初始化微调后的模型。训练过程中使用了所有表达值为零和非零的基因标记。

扰动响应预测
为了针对扰动预测任务对模型进行微调，scGPT在模型训练前选择了高变基因（HVGs）并对表达值进行了预处理。预训练模型中嵌入层和Transformer层的参数被用于初始化微调后的模型。在微调过程中，纳入了所有表达值为零和非零的基因标记。针对扰动预测任务的输入，scGPT做了两个显著的改变：首先，使用经过 “log1p” 变换的表达值作为输入和目标值，而不是分箱后的值，以便在这项任务中更好地预测扰动后的绝对表达值。其次，scGPT在每个输入基因的位置附加了一个二元条件标记，用以表明该基因是否受到了扰动。scGPT采用了扰动基因表达预测（perturb-GEP）的微调目标，并对训练设置做了进一步修改。scGPT没有使用同一细胞的掩码和未掩码表达值作为输入和学习目标，而是使用一个对照细胞作为输入，将受扰动的细胞作为目标。这是通过将每个受扰动的细胞与一个未受扰动的对照细胞随机配对来构建输入 - 目标对实现的。输入值由对照细胞中的所有基因表达值组成。因此，该模型学会了基于对照基因的表达情况和扰动标记来预测扰动后的响应。

整合多个 scRNA-seq 数据集时的批次校正
当输入的原始计数值矩阵包含来自不同测序批次或技术的多个数据集时，批次效应可能会成为细胞类型聚类中的一个主要干扰因素。因此，scGPT的目标是在整合多个 scRNA-seq 数据集的同时，校正批次效应，且保留生物学差异。为了针对这项整合任务进行微调，预训练基础模型与当前数据集之间共有的基因标记集被保留下来。进一步从这个共有集中选择了一部分高变基因（HVGs）作为输入。在模型训练前，对表达值进行了预处理，这与细胞类型注释任务中的处理方式类似。所有预训练模型的权重都被用于初始化微调后的模型。默认情况下，训练过程中使用了所有表达值为零和非零的基因标记。除了基因表达预测（GEP）和用于细胞建模的基因表达预测（GEPC）之外，弹性细胞相似度（ECS）、通过反向传播进行的领域自适应（DAR）以及特定领域批量归一化（DSBN）这些微调目标也同时进行了优化，以便通过反向传播和特定领域的归一化来增强细胞的对比学习并实现明确的批次校正。

单细胞多组学数据的整合表征学习
单细胞多组学数据可能在不同实验批次中包含不同的测序模态。对于单细胞多组学数据，scGPT研究了两种数据整合设置，即配对设置（paired）和镶嵌设置（mosaic）。在配对设置中，所有样本（细胞）的模态是联合测量的。在镶嵌设置中，一些批次共享几种常见的数据模态，但并非全部。由于存在额外的染色质可及性（ATAC）或蛋白质token，scGPT仅继承了针对RNA数据训练好的基因嵌入，而对额外的标记嵌入以及模型的其余部分进行从头训练。如果数据集包含额外的蛋白质数据，那么训练过程中仅使用表达值非零的标记。否则，默认情况下会使用表达值为零和非零的标记。scGPT使用了一组额外的模态标记来表明每个标记的数据类型（即基因、区域或蛋白质），并在基因表达预测（GEP）和用于细胞建模的基因表达预测（GEPC）微调目标中辅助掩码基因和值的预测。默认情况下，使用基因表达预测（GEP）和用于细胞建模的基因表达预测（GEPC）微调目标对模型进行优化。如果存在多个批次，则会纳入通过反向传播进行的领域自适应（DAR），以助力多模态的批次校正。

基因调控网络推断
对于基于基因嵌入的基因调控网络（GRN）推断，在零样本设置下，根据KNN算法，利用scGPT预训练的基因嵌入构建了基因相似度网络。在微调设置下，以类似的方式，从在人类免疫数据集上微调后的scGPT模型构建了基因相似度网络。遵循塞利亚（Ceglia）等人的方法，进一步在相似度图上进行了莱顿（Leiden）聚类分析，并从由五个及以上基因组成的基因簇中提取了基因程序。

对于基于注意力机制的目标基因选择，在亚当森（Adamson）扰动数据集上对scGPT血液模型进行了微调，该数据集包含对一种白血病细胞系进行的87项CRISPR扰动实验。在图6a中展示了目标基因选择流程。

• 图6a：基于注意力机制的针对扰动数据的基因调控网络（GRN）发现工作流程。获取处于对照细胞状态和受扰动细胞状态下的注意力分数，并依次按行和按列进行排序。相应地，在“对照”“受扰动”和“差值”设置下选择受影响最大的基因。

对于每个感兴趣的受扰动基因，scGPT首先分别向模型输入受扰动细胞集和对照细胞集，从而获取两组注意力图，即受扰动的注意力图和对照的注意力图。请注意，原始注意力分数是从模型最后一个注意力层的所有八个注意力头中获得的。然后，原始注意力分数要经过两轮rank归一化处理，首先按行进行，然后按列进行。接着，对经过rank归一化的注意力分数在八个注意力头之间取平均值，以输出一个聚合注意力图。这就得到了用于选择受影响最大基因的最终注意力图。

对于每个感兴趣的受扰动基因，通过对最终注意力图中该受扰动基因所在列的分数进行排序，来选择受其影响最大的基因。这体现了这样一种直觉，即注意力图中的列表明了感兴趣的基因对其他基因的影响程度。scGPT提供了三种受影响最大基因的选择设置，即从对照注意力图中选择“对照”情况的基因，从受扰动注意力图中选择“受扰动”情况的基因，以及从两者的差值中选择“差值”情况的基因。从对照注意力图中选择的目标基因应反映出感兴趣的基因所参与的基础通路，而受扰动注意力图则反映了扰动后的影响。这两个注意力图之间的差值应突出显示基因网络中从扰动前到扰动后变化最大的边。

scGPT benchmark设置

单细胞RNA测序（scRNA-seq）的细胞类型注释

在髓系细胞、多发性硬化症（MS）以及人类胰腺数据集上，将scGPT与两种近期基于Transformer的细胞类型注释方法（scBERT和TOSICA）进行了基准测试对比。对于每个数据集，使用参考数据分区来进行模型训练和验证。获取查询集上预测的细胞类型标签以进行评估。基于四个标准分类指标来评估细胞类型分配的性能：准确率、精确率、召回率和宏平均F1值。准确率、精确率和召回率是从全局角度计算总体性能，而宏平均F1值是按类别进行平均，以增加稀有细胞类型的权重。还报告了一个归一化的混淆矩阵，其中包含按细胞类型划分的“精确率”以提供更多细节。

单细胞RNA测序（scRNA-seq）的扰动预测

将scGPT与近期的扰动预测方法GEARS以及一个线性回归模型进行了比较。线性回归模型将来自对照细胞的基因表达值以及每个基因上扰动状态的二元编码作为输入特征。该模型通过最小化回归误差，使用输入特征的线性组合来估计每个基因在扰动后的表达值。为确保一致性，遵循了鲁哈尼（Roohani）等人在其基准测试中概述的预处理步骤。25%的扰动被划分到测试集中，在训练期间不被模型看到。虽然在GEARS的研究中，CPA方法被认为效果欠佳，但我们在实验中证实了这些发现。然而，我们认为这种性能差距在很大程度上是由于实验设置的差异造成的。CPA主要是为了学习适用于未见过的细胞类型的共享扰动嵌入而设计的，这与我们关注的完全未见过的扰动有所不同。为确保公平的基准测试，我们在最终的比较中排除了CPA方法。我们在相同的设置下训练了所有模型，并报告了评估结果。最初，基因表达值是根据所有基因的总计数对每个细胞进行归一化处理，然后应用对数变换。随后，我们选择了5000个高变基因（HVGs），并将任何最初未考虑的受扰动基因纳入基因集。在实验中，对于所有三个数据集中的单基因扰动，对扰动进行划分以确保测试扰动在训练中未出现，也就是说，训练集中没有细胞经历过任何测试扰动。对于诺曼（Norman）等人数据集中的双基因扰动，训练-测试划分包括三种难度逐渐增加的场景：（1）训练集中两个未见过的基因都没有；（2）训练集中两个未见过的基因中有一个；（3）训练集中两个未见过的基因都有。为评估扰动预测的准确性，检查了扰动后与对照细胞状态之间的表达变化（“delta”）。我们计算了数据中预测的和观察到的表达变化之间的皮尔逊相关系数，记为Pearson delta。我们还报告了在差异表达最显著的前20个基因上的这些皮尔逊指标。因此，我们给出了针对差异表达条件的额外评估指标，即在差异表达基因上的Pearson delta。此外，我们测试了预测的基因表达与真实表达值之间的皮尔逊相关系数（Pearson）的评估。我们发现基于delta的评估比仅仅与真实表达水平的相关性更能真实地反映模型的性能。

基于聚类的生物学验证

我们首先从scGPT模型中获取每个扰动条件下的代表性基因表达谱。scGPT通过对数据集中所有对照细胞的基因表达求平均得到的单个样本对照基因表达向量（即大小为1×M个基因），预测每个扰动条件下的代表性扰动响应。诺曼等人的数据集包含105个独特的受扰动基因，这总共产生了5565个独特的扰动组合需要预测。将高维的预测扰动响应投影到二维的UMAP图上。首先将UMAP图与诺曼等人在原始论文中发现的功能组进行比较，在该论文中，236个扰动实验是基于真实的扰动响应进行聚类的，并通过标记基因表达对其功能角色进行注释。我们检查了scGPT预测的UMAP投影与原始论文中发现的功能分组之间的一致性。然后，分析了scGPT预测的UMAP图中存在的子簇。使用莱顿聚类分辨率为0.5时，在预测扰动响应的UMAP图中识别出了54个子簇。用出现次数最多的受扰动基因标注每个簇，并将其作为该簇的主导基因。

---

使用诺曼（Norman）数据集进行了聚类分析，以验证与生物学相关的功能信号。最初的Perturb-seq研究涵盖了针对105个基因的236种扰动。然而，考虑到这些目标基因的所有可能组合，总共有5565种潜在的扰动，这表明实验性的Perturb-seq数据仅占整个扰动空间的5%。因此，应用经过微调的scGPT在计算机模拟中扩展扰动，并使用均匀流形近似与投影（UMAP）在图3d中可视化了每种扰动的预测平均响应。利用原始研究中的注释，发现相同功能组的扰动条件聚集在相邻区域（补充图4）。

接下来，我们使用莱顿（Leiden）聚类算法对预测的表达进行聚类，并且观察到这些聚类与扰动组合中的“主导基因”具有高度相关性。例如，与KLF1基因相关的圈定聚类表明，该聚类中的数据点经历了涉及KLF1和另一个基因（即KLF1 + X）的组合扰动。以KLF1和CNN1聚类为例，进一步验证了相应的预测表达仅在这些区域中较高（图3e），这与诺曼数据集中CRISPRa（CRISPR介导的转录激活）Perturb-seq实验的预期结果一致。主导基因聚类展示了scGPT揭示扰动组合之间关联的能力。

• d，扰动条件下预测的基因表达谱的UMAP图。该UMAP图按莱顿聚类进行着色，并标注有每个簇的主导基因。e，在扰动条件的UMAP图上，两个选定的受扰动基因（KLF1和CNN1）的表达模式。

---

反向扰动预测任务

从诺曼数据集中选择了20个基因，构建了一个扰动用例，以微调并测试一个新的扰动模型。选择这20个基因的子集是为了在基于scGPT的训练-测试划分的基础上，使训练和测试用例中真实扰动数据的比例最大化。这个选定的子空间包含210个独特扰动组合中的39个训练用例、3个验证用例和7个测试用例。其余的是没有实验结果的未见过的用例。反向扰动预测遵循前K个检索任务设置：

我们使用所有210个扰动条件下的预测响应作为参考数据，将7个测试用例的真实响应作为查询集。目标是检索出产生与查询结果最相似响应的前几个扰动条件。对于参考数据，我们不是为每个扰动条件获取单个代表性基因表达谱，而是从30个随机采样的对照细胞中获取预测响应，以增加多样性。这产生了一个包含6300个预测的扰动后基因表达谱的参考数据库。对于X + Y扰动的每个测试用例，我们使用所有经历了X + Y扰动的细胞的真实基因表达谱作为查询集。对于前K个检索，我们设计了一个集成投票策略，涉及两轮选择：

• 在第一轮中，每个单独的查询细胞通过与参考数据集的欧几里得距离选择其最相似的前K个表达谱。
• 在第二轮中，我们根据所有查询细胞的投票数对候选扰动条件进行排名。将集成投票后第二轮排名中得票数最多的前K个扰动条件报告为预测的源扰动条件。

单细胞RNA测序（scRNA-seq）的批次整合

在这项工作中，将scGPT的性能与其他三种方法（即Seurat、Harmony和scVI）的性能进行了比较。评估内容涵盖了对三个整合数据集（COVID-19、外周血单个核细胞（PBMC）10k和鼻周皮质）的批次校正和细胞类型聚类。Harmony和scVI在近期的整合基准测试中被突出显示为表现最佳的方法。为确保公平比较，所有方法都使用相同数量的1200个高变基因（HVGs）作为输入。基因表达值通过考虑所有基因的总计数对每个细胞进行归一化处理，然后进行对数变换。在训练完成后获取整合的细胞嵌入，并用于评估。

使用生物学保守性指标对整合的细胞嵌入进行评估。这些指标包括细胞归一化互信息（NMI）、细胞调整兰德指数（ARI）和细胞轮廓系数（ASW）。这些分数衡量了导出的细胞类型簇与真实标签之间的一致性。批次校正性能使用批次聚类的平均轮廓宽度的倒数（记为ASWbatch）和图连通性度量（记为GraphConn）进行量化。

单细胞多组学整合

我们在配对和镶嵌两种整合设置下，将scGPT与近期的单细胞多组学整合方法Seurat（v.4）、scGLUE和scMoMat进行了基准测试对比。在配对数据整合实验中，我们以10x Multiome PBMC43数据集（包含RNA和染色质可及性测序（ATAC-seq）数据）为例，将scGPT与scGLUE和Seurat（v.4）进行了基准测试对比。所有方法都使用相同设置的1200个高变基因和4000个高变峰作为输入。进一步在BMMC数据集（包含配对的RNA和蛋白质）上，将scGPT与Seurat（v.4）进行了基准测试对比。为公平比较，在这种情况下我们没有对scGLUE进行基准测试，因为该方法并非专门为蛋白质数据建模而设计。同样，使用相同的1200个高变基因和所有134种蛋白质作为输入。在镶嵌数据整合实验中，在ASAP PBMC数据集上，将scGPT与scMoMat进行了基准测试对比。总共使用1200个高变基因、4000个高变峰和所有216种蛋白质特征作为两种方法的输入。在保持输入特征集一致的同时，使用每种方法的自定义预处理管道对表达值进行归一化处理。在训练后获取整合的细胞嵌入用于评估。

基因调控网络推断-embedding

针对已知的人类白细胞抗原（HLA）和分化簇（CD）基因网络，对scGPT基因嵌入相似度网络进行了验证。对于每个网络，首先通过过滤带有集合前缀（即HLA-和CD-）的基因名称来定义相关基因集。然后，从Reactome 2022数据库（https://reactome.org/）中过滤出参与免疫系统R-HSA-168256通路的基因。对于CD基因，为便于比较预训练模型和微调模型，使用了与来自人类免疫数据集的高变基因的共同基因集。然后，我们从scGPT模型中提取这些选定基因的基因嵌入，并构建一个k近邻（kNN）相似度网络。我们通过选择余弦相似度大于特定阈值（即HLA网络为0.5，CD基因网络为0.4）的边来突出显示强连接的子网。然后，将这些子网与免疫系统中已知的功能组进行比较。

此外，通过通路富集分析验证了scGPT模型提取的基因程序的质量。scGPT将每个基因程序作为输入基因列表，并选择具有统计学显著性的通路作为“通路命中”。根据执行的测试总数（即基因程序的数量乘以通路测试的数量），使用邦费罗尼校正（https://mathworld.wolfram.com/BonferroniCorrection.html）将P值阈值调整为0.05。报告了在Reactome 2022数据库（https://reactome.org/）中通路命中的数量。

作为基准，将结果与从基线共表达图中提取的基因程序进行了比较。共表达图是通过人类免疫数据集中基因的归一化基因表达之间的皮尔逊相关性来定义的。为确保与scGPT网络具有相似的模块性，我们将该图稀疏化为一个kNN相似度网络（k = 15）。按照与scGPT相同的流程，通过莱顿聚类从基因簇中识别基因程序。作为敏感性分析，我们报告了在莱顿分辨率为1、5、10、20、30、40、50和60时，scGPT和共表达方法的通路命中情况。我们进一步检查了在莱顿分辨率为40时，每种方法识别出的共同通路和独特通路，以更深入地了解性能差异。

基因调控网络推断-attention

在ChIP-Atlas数据库中验证了scGPT基于注意力机制的受影响最大基因选择方法，该数据库包含已知转录因子的经过实验验证的基因靶点。首先通过将亚当森（Adamson）扰动数据集中的受扰动基因列表与ChIP-Atlas进行交叉核对，选择了两个示例转录因子（由DDIT3和BHLHE40编码）。对于每个转录因子，通过将其与经过验证的基因靶点进行比较，验证了在“差值”设置下通过注意力机制选择的前20个受影响最大的基因。请注意，在“差值”设置下，通过检查受扰动注意力图和对照注意力图之间的差异，基于扰动后的变化来选择前20个基因。真实基因靶点列表是从ChIP-Atlas中通过过滤转录起始位点位于转录因子峰调用区间10千碱基对距离内的人类基因（hg38）获得的。我们报告了前20个通过注意力机制选择的基因靶点与真实靶点基因之间的重叠数量。

随后，通过检查所选择的前100个基因之间的重叠情况，比较了三种受影响最大基因选择方法（即对照、受扰动和差值）。这三个前100个基因集的重叠和差异情况在维恩图中进行了可视化展示。我们进一步在Reactome数据库中验证了这些前100个基因与转录因子一起参与的通路。通路命中情况和基因重叠百分比在热图中进行了可视化展示。