提供代码 | 学习如何通过单细胞数据探索癌症纵向演化路径 | Nat.Genet

Basic Information

• 英文标题：Deciphering the longitudinal trajectories of glioblastoma ecosystems by integrative single-cell genomics
• 中文标题：通过整合单细胞基因组学解码胶质母细胞瘤生态系统的纵向轨迹
• 发表日期：09 May 2025
• 文章类型：Article
• 所属期刊：Nature Genetics
• 文章作者：Avishay Spitzer | Mario L. Suvà
• 文章链接：https://www.nature.com/articles/s41588-025-02168-4

Abstract

Para_01

1. 标准治疗后IDH野生型胶质母细胞瘤（GBM）的进化仍知之甚少。
2. 在这里，我们使用单核RNA测序和批量DNA测序分析了59名患者的匹配原发性和复发性GBMs，评估了GBM生态系统在细胞和分子异质性各层中的纵向进化。
3. 最一致的变化是复发时恶性细胞的比例降低，而肿瘤微环境（TME）中的胶质细胞和神经元细胞类型则出现相反的增加。
4. 大多数匹配对之间的主要恶性细胞状态不同，但没有任何状态在任一时间点上是独有的或高度富集的，也没有在整个队列中存在一致的纵向轨迹。
5. 然而，在部分患者中，特定的轨迹被富集。
6. 恶性状态丰度的变化反映了TME组成和基线特征的变化，这反映了GBM生态系统的共同进化。
7. 我们的研究提供了GBM多样化纵向轨迹的蓝图，并突出了影响这些轨迹的治疗和TME修饰因素。

Main

Para_01

1. 尽管基础和转化研究取得了显著进展，IDH野生型胶质母细胞瘤（GBM）仍无法治愈。
2. 新诊断的GBM的标准治疗包括最大手术切除，随后进行放疗，并同时使用烷基化剂替莫唑胺（TMZ）。
3. 这种治疗通常会引发部分反应，特别是当O6-甲基鸟嘌呤-DNA甲基转移酶（MGMT）基因的启动子被甲基化时，因为该基因编码一种可以移除TMZ添加到O6-甲基鸟嘌呤上的烷基基团的DNA修复蛋白。
4. 然而，疾病不可避免地复发，中位总生存期为14.6个月。
5. 为了探究复发机制，一些研究分析了GBM样本的纵向演变。
6. 胶质瘤纵向分析（GLASS）联盟报告称，原发性和复发性GBM在遗传驱动改变方面高度相似，表明标准治疗并未产生显著的克隆选择。
7. 相反，转录组和蛋白质组研究显示了更动态的进化模式，例如神经前体向间质（MES）的转变、增强的神经信号传导以及复发时细胞外基质（ECM）相关基因的上调。
8. 然而，已发表的纵向GBM研究通常依赖于整体组织分析（在恶性细胞和肿瘤微环境（TME）区域之间平均信号），缺乏独立探究恶性细胞和TME细胞类型及细胞状态演化的细节。
9. 单细胞纵向研究的队列规模有限，临床数据也有限，阻碍了对单细胞图谱与临床相关性的直接分析。
10. 这些限制，加上协调不同研究小组结果的挑战，这些研究采用了不同的分析和分类方案，突显出需要一个更大规模的单细胞纵向GBM数据集，具有详细的临床和基因组注释，并具备稳健的分析。

Para_02

1. 我们在这里提出了两项研究来解决这些挑战。
2. 我们建立了GBM细胞抗性和进化（CARE）联盟，以在单细胞分辨率下表征GBM的纵向进化轨迹。
3. 我们对59名患者的配对原发性和复发性GBM样本进行了单核RNA测序（snRNA-seq），这些患者有详细的临床注释和遗传特征。
4. 在参考文献15的配套研究中，我们对数据进行了独立于采样时间（原发性与复发性）的分析，表征了GBM转录生态系统在多个层面的情况。
5. 在这项工作中，我们研究了原发性与复发性之间这些生态系统的纵向变化，并识别出修饰其轨迹的分子和临床特征。
6. 我们发现，在复发时，恶性细胞比例较低，而胶质细胞和神经元（即胶质-神经）肿瘤微环境（TME）细胞类型的比例如同增加，这是复发中最一致的纵向变化，出现在66%的病例中。
7. 与之前的研究相反，我们在整个队列中未观察到复发时一致的恶性GBM状态富集。
8. 然而，特定的轨迹与某些病例子集相关，例如具有MGMT甲基化或放疗相关的小缺失表型的GBMs。
9. 值得注意的是，恶性状态丰度的变化与基线特征和TME组成的相应变化相呼应，反映了GBM生态系统的相互依赖性变化。
10. 我们的研究突出了由治疗反应和TME相互作用塑造的多样化纵向GBM轨迹。

Results

Longitudinal profiling of a clinically annotated GBM cohort

临床注释的GBM队列的纵向分析

Para_01

1. 我们从原发和复发手术中收集了匹配的冰冻GBM标本（图1a和补充表1）。
2. 肿瘤样本来自七家医院，总共对59名患者的121个肿瘤进行了snRNA-seq和全外显子组测序（WES）或全基因组测序（WGS；图1b；方法）。
3. 对于大多数患者，我们对两个时间点的样本进行了分析（标记为T1和T2）；对于三名患者，我们对三个时间点的样本进行了分析（T1、T2和T3）。
4. 大多数匹配的样本对来源于原发肿瘤和首次复发肿瘤，除了三对样本由两个复发样本组成，四对样本由原发肿瘤和第二次复发肿瘤组成。

Fig. 1: Longitudinal profiling of a clinically annotated GBM cohort.

图片

- 图片说明

◉ a，描述本研究工作流程的示意图。TMZ，替莫唑胺；RT，放疗；int.，中间。b，本数据集中所有59名患者的临床和技术信息，包括手术间隔、诊断时年龄、性别、复发部位（rec）、MGMT启动子状态、治疗和组学数据类型。c，基于我们队列中细胞的基因表达值进行降维的UMAP图，按细胞类型（顶部）和时间点（底部）着色。d，检测到遗传异常的患者在数据集中的比例（%），按时间点着色。顶部的环形图表示平均体细胞突变分布，无论是在两个时间点（T1 + T2）均检测到改变、仅在原发肿瘤（T1）或仅在复发肿瘤（T2）检测到改变。底部的柱状图用于已知的驱动基因组事件及其与基因组改变相关的代表性表型—臂级指整个染色体臂扩增（amp）（chr7q）/缺失（del）（chr10q），CNA指高程度的基因扩增或缺失，SNV指单核苷酸突变。新生高突变（HM）表示在两个时间点均具有高突变负荷，但没有治疗相关突变特征的证据。小缺失表示在复发时小缺失负荷增加的肿瘤。UMAP，统一流形近似和投影；chr，染色体；assoc.，相关。◉ ,

Para_02

1. 51名患者的肿瘤位置已知，其中45例在手术边缘局部复发，6例在切除腔以外的部位复发。
2. 大多数患者（51/59）在手术后接受了同步替莫唑胺的放疗以及辅助替莫唑胺治疗。
3. 与历史队列相比，我们队列中的临床结果有所改善（两次手术之间的中位间隔时间为12.4个月，总生存期中位数为24个月），这可能反映了适合接受多次手术的患者身体状况更好，且诊断时年龄更小（中位数为53岁，范围为19–73岁）。

Para_03

1. 在我们的10x Genomics snRNA-seq数据集中，429,305个细胞通过了质量过滤。
2. 这些细胞根据推断的染色体拷贝数异常（CNA）被分类为恶性（246,408个细胞）或非恶性（182,897个细胞）。
3. 非恶性细胞基于已建立的细胞类型特征和标记基因表达进行了进一步分类和注释。
4. 非恶性细胞在时间点之间聚集在一起。
5. 我们没有观察到不同时间点之间分析的细胞数量或每个细胞检测到的基因数量有显著差异。

Para_04

1. 突变和CNA谱型是从批量DNA测序数据中获得的（图1b）。
2. GBM驱动改变的频率与先前的研究一致，且在时间点上基本保持稳定（图1d和扩展数据图1c）。
3. 然而，在复发时，26%的病例（54名患者中的14例）出现了与治疗相关的DNA高突变或小缺失负担增加（图1d，右侧）。

Longitudinal changes in cell-type composition

细胞类型组成的纵向变化

Para_01

1. 全局（非匹配）的细胞类型丰度比较显示，时间点之间总体相似，除了复发肿瘤中恶性细胞的丰度显著较低，而胶质神经TME（少突胶质细胞、神经元和星形胶质细胞（ACs））的丰度较高（图2a和扩展数据图2a）。
2. 我们进一步使用非负矩阵分解（NMF；方法）将广泛的非恶性细胞类型（例如，髓系细胞）细分为更具体的亚型（例如，小胶质细胞和骨髓来源的巨噬细胞）。
3. 对这些高分辨率细胞状态的分析表明，复发肿瘤中富集了反应性ACs以及表达干扰素程序的少突胶质细胞（扩展数据图2b）。
4. 为了进一步评估细胞类型组成的差异，我们根据恶性细胞和TME细胞类型的相对丰度将肿瘤分为六种组成组（方法）。
5. 简而言之，根据恶性细胞的比例，肿瘤被分为以下三组：高恶性比例（HMF）、中间恶性比例（IMF）和低恶性比例（LMF）。
6. 后一组根据主导的TME细胞类型进一步细分（例如，LMF-肿瘤相关巨噬细胞（TAM）、少突胶质细胞（OC）和胶质神经（GN）；扩展数据图2c，左侧）。
7. HMF组在原发肿瘤中富集，而具有胶质神经（LMF-GN）或少突胶质细胞TME（LMF-OC）的LMF组在复发时富集（两组P=0.005，双侧Fisher检验；图2b，左侧）。

Fig. 2: Prognostically significant longitudinal cell-type changes in the TME, global similarity of cell states between primary and recurrent tumors.

图片

- 图片说明

◉ 在两个时间点（原发和复发，n = 118）中，每个细胞类型的占比（%），其中观察到统计学显著差异。P 值通过双侧 Wilcoxon 秩和检验计算。◉ b，区分组别的特征及其在原发/复发状态下的分布情况，包括组成（左）、恶性状态频率（中）和基线特征（右）。P 值通过 Fisher 精确检验计算。◉ c，Kaplan–Meier 曲线显示了根据第二次手术时的转录组特征，第二次手术后的生存时间。左图，12 名具有胶质神经组成组复发的患者与 47 名其他类型复发的患者。右图，6 名具有 MES/缺氧细胞状态组复发的患者与 53 名其他类型复发的患者。X 轴反映从第二次手术到结果的时间。每张图中两组之间生存差异的统计学意义通过对数秩检验计算。◉ d，热图显示了转录组特征之间的皮尔逊相关性。相关性分别在原发样本（左上）和复发样本（右下）中计算，并按照复发样本的层次聚类模式排序。NS，无显著性。

Para_02

1. 为了验证我们的发现，我们重新分析了一个独立的纵向GBM单细胞RNA测序队列。一致地，我们观察到复发时恶性细胞数量较低，少突胶质细胞数量较高，以及LMF组成组的富集程度相似。

Para_03

1. 按组成组分层的生存分析显示，第二次手术后的生存时间存在显著差异（扩展数据图2g，顶部），但第一次和第二次手术之间的时间间隔没有显著差异。
2. 复发时具有胶质神经TME的患者，第二次手术后的生存时间明显长于其他复发的患者（中位数为19个月与8个月，P = 0.012，对数秩检验；图2c，左侧）。
3. 这些临床过程的差异与其他临床特征无关，包括MGMT甲基化状态（扩展数据图2h，顶部）。

Longitudinal stability of transcriptional features

转录特征的纵向稳定性

Para_01

1. 接下来，我们研究了恶性细胞状态在不同时间点的丰度变化。
2. 如附带论文中更详细描述的那样，我们对恶性细胞应用了NMF，并识别出十个GBM元程序（MPs），代表肿瘤内细胞异质性（方法）。
3. 这包括细胞周期以及之前识别出的AC样、MES样、应激、缺氧、神经前体细胞样（NPC样）和少突胶质前体细胞样（OPC样）状态。
4. 此外，我们还识别出以下三种之前未描述的状态：分化神经元状态（NEU样）、胶质前体样状态（GPC样）和纤毛相关状态（Cilia样），这些状态通过标记基因、通路富集和与正常神经细胞类型的相似性进行表征。
5. 除了细胞状态外，我们定义了三种基线基因表达谱（BP），反映了不受恶性细胞状态丰度差异驱动的肿瘤间异质性。
6. 这三个BP被注释为胶质-BP、ECM-BP和神经元-BP。

Para_02

1. 在我们的队列和参考文献14的队列中，原发样本和复发样本之间的恶性细胞状态丰度比较未显示显著差异（扩展数据图3a,b）。
2. 此外，当我们根据主要恶性状态或BPs将肿瘤分为六组时，在时间点上没有检测到显著差异（图2b（中间和右侧），扩展数据图2c（中间和右侧）和扩展数据图3c,d）。
3. 然而，按主要恶性细胞状态分层的生存分析显示，复发时MES样/缺氧组的预后明显较差（第二次手术后的中位生存时间为5个月与10.3个月，P = 0.00086，log-rank检验；图2c（右侧）和扩展数据图2g,h），这与最近的一项研究一致。

Para_03

1. 细胞类型组成的多个转录组层次、恶性细胞状态和恶性生物过程之间存在显著关联，并定义了三种典型的胶质母细胞瘤生态系统（即胶质型、神经元型和细胞外基质型）。
2. 我们接下来研究了胶质母细胞瘤生态系统中各组分之间的相互作用是否在原发性和复发性肿瘤之间保持一致。
3. 我们分别计算了每个时间点的以下相关性：（1）11种肿瘤微环境细胞类型的比例，（2）十种恶性细胞状态的比例以及（3）三个生物过程评分。
4. 在两个时间点，我们发现了三组可比较的相关特征（图2d和扩展数据图3e）。
5. 这些保守的关联（通过置换检验P < 0.01）基本保持稳定。

Para_04

1. 然而，我们确实观察到一些差异，例如神经元高生态位与OPC样恶性细胞（在原发样本中）和NEU样恶性细胞（在复发样本中）的偏好性关联（图2d和扩展数据图3e）。
2. 有趣的是，使用CellChat进行的配体-受体分析显示，在不同时间点预测的相互作用存在一些差异。例如，由少突胶质细胞和神经元表达的神经元粘附分子（NLGN）和神经调节蛋白（NRG）基因家族，被预测与NEU样和NPC样恶性细胞状态中的神经元素（NRXN）基因家族或ERBB4发生相互作用，优先在复发时出现（扩展数据图3f）。
3. 总体而言，我们队列的全局分析（未按样本对匹配）显示了细胞类型组成的差异，这些差异与恶性细胞比例的降低以及复发性GBM细胞与脑实质的进一步整合一致，但并未揭示恶性细胞状态或BPs方面的显著全局差异。

Divergence of individual sample pairs

个体样本对的差异

Para_01

1. 尽管细胞类型和恶性细胞状态在全球范围内是保守的，但在检查个体匹配样本对时，我们发现存在广泛的差异（图3a）。
2. 大多数（68%）样本对在诊断时属于一个组成组，而在复发时则属于不同的组成组。
3. 同样，64%和56%的样本分别从主导恶性状态组和BP组切换到了不同的组。
4. 我们在参考文献14的队列中观察到了同样的现象，其中68%的患者在复发时切换了组成组，72%切换了恶性组（扩展数据图4a）。
5. 虽然没有一组显著富集，但复发时神经元-BP的比例有增加的趋势，而胶质-BP的比例减少（图3a（右））。
6. 值得注意的是，我们观察到所有可能的轨迹至少出现了一次（63/88条轨迹），包括相反的轨迹。
7. 例如，患者P1在复发时AC样恶性细胞增加而GPC样恶性细胞减少，而患者P2则表现出相反的轨迹（扩展数据图4b）。
8. 此外，那些在三个时间点进行分析的患者在两个治疗间隔（T1-T2和T2-T3）之间表现出相反的轨迹，尽管复发发生在相同部位（图3b）。
9. 最后，局部复发和远处复发之间的保守频率没有显著差异（P = 0.66，双侧Fisher检验）。

Fig. 3: Divergence of individual sample pairs.

图片

- 图片说明

◉ a，桑基图展示了时间点之间在组成（左）、恶性状态（中）和基线特征组（右）之间的转换。特征的宽度表示具有指定分类的肿瘤的比例。浅灰色链接反映了从一个组（T1）到T2不同组的转换。◉ b，纵向转录组变化在两个在三个时间点进行分析的特征中突出显示了不可预测的基因表达轨迹。详细的临床信息包括治疗方案和肿瘤MRI（顶部），三种转录组层的比例和评分热图（中部）以及时间点之间的突变分数（底部）。INF，干扰素；BEV，贝伐珠单抗；T1Gd，T1加权增强图像。◉ c，观察到的（彩色条）与预期的（黑色点）在时间点上的保守性。P值通过二项式检验计算，并使用Holm方法进行多重检验校正。误差条代表来自二项式检验的置信区间。星号表示Holm校正后的P < 0.01。MRI，磁共振成像。

Para_02

1. 这种个体样本对的差异促使我们比较观察到的保守模式与一个零模型中的预期模式，该模型中匹配的原始-复发对彼此之间的相似性与随机选择的样本对相同。
2. 这项分析显示，在所有三个层次上，匹配的对比样本比随机预期的更加保守（每个层次的P值均小于0.004，二项式检验；图3c）。
3. 因此，匹配的样本通常在转录组分类上出现差异，并且可能在患者之间和患者内部呈现出广泛的轨迹差异，但它们仍显著更有可能保留任何特定的分类，相比随机选择的样本对。
4. 这一观察结果促使我们进一步对队列进行分层，并尝试识别进化轨迹的预测因子/修饰因子。

MGMT expression in GBM single cells influences recurrence trajectories

GBM单细胞中MGMT的表达影响复发轨迹

Para_01

1. 先前的研究表明，甲基化的MGMT启动子患者对TMZ的反应更好，预后更佳，但这种差异是否在细胞结构或治疗后的细胞轨迹中有所体现尚不清楚。
2. 我们根据原发肿瘤的MGMT启动子甲基化状态将患者分为甲基化（MGMT-MET，n = 12）、未甲基化（MGMT-UM，n = 19）或未知（n = 28）三组。
3. 正如预期的那样，MGMT-MET患者的手术间隔时间比MGMT-UM患者更长，这在我们的队列和GLASS队列中均得到证实。

Para_02

1. 我们接下来评估了MGMT启动子甲基化状态是否与单细胞中的MGMT表达相关，因为这种关联在批量分析中一直难以确定。
2. 在原发性MGMT-MET样本中，恶性细胞的MGMT表达显著低于原发性MGMT-UM样本（P = 5.3 × 10−6，双侧Wilcoxon秩和检验；图4a），但在TME细胞类型中未观察到这种效应（图4b）。
3. 在部分MGMT-MET病例中，复发时MGMT表达增加，这可能反映了低MGMT表达细胞在治疗过程中处于竞争劣势。

Fig. 4: MGMT expression in GBM single cells influences recurrence trajectories.

图片

- 图片说明

◉ 每份样本中恶性细胞的平均MGMT表达量随时间点的变化。只有在甲基化检测中验证了MGMT启动子状态，并根据原发样本的状态标注为甲基化/非甲基化的样本才被纳入分析（n = 12个甲基化和n = 19个非甲基化）。虚线连接匹配的样本。统计显著性通过双侧Wilcoxon秩和检验估计；显示的是未调整的P值。在a中，组内比较使用了Wilcoxon符号秩检验（配对）。◉ 甲基化和非甲基化样本中每种细胞类型的MGMT表达水平。统计显著性通过双侧Wilcoxon秩和检验估计；显示的是未调整的P值。◉ 根据MGMT表达水平分层的本分析中患者的生存曲线。X轴表示第一次和第二次切除之间的时间。各图中两组生存差异的统计显著性通过log-rank检验计算。◉ 根据MGMT表达水平分层的MGMT-MET患者的生存曲线。X轴表示第一次和第二次切除之间的时间。各图中两组生存差异的统计显著性通过log-rank检验计算。◉ 根据MGMT启动子甲基化状态分层的高MGMT表达组的生存曲线。X轴表示第一次和第二次切除之间的时间。各图中两组生存差异的统计显著性通过log-rank检验计算。◉ NOS，未另行说明。◉ 本研究及参考文献14中，根据MGMT表达水平分层的两个时间点（原发和复发）中MES样细胞状态比例（%）的差异。所有样本点均以抖动点形式展示。统计显著性通过双侧Wilcoxon秩和检验估计；显示的是未调整的P值。◉ 纳入分析的患者的MGMT特异性特征图。数值表示每个患者的T2–T1差异（基因表达的log2(FC)）。患者按特征评分差异排序。基因按其与特征评分差异的相关性（在患者中）排序。FC，倍数变化。

Para_03

1. MGMT启动子甲基化与低MGMT表达之间的相关性促使我们构建了一个逻辑回归模型，该模型使用恶性细胞的平均MGMT表达水平来推断MGMT甲基化状态（扩展数据图5c）。
2. 在应用此模型时，我们发现根据MGMT表达水平进行分类与手术切除间隔有关（图4c），事实上，其关联性比仅根据MGMT启动子甲基化进行分类更好——在MGMT-MET组中，与MGMT高表达/中等表达组相比，MGMT低表达组的手术间隔显著更长（图4d）；相反，在MGMT高表达组中，进一步根据MGMT启动子甲基化状态分层时，手术间隔没有显著差异（图4e）。
3. MGMT表达水平表现的改善可能是由于排除了可能影响MGMT启动子甲基化评估的非恶性细胞；或者，这可能表明治疗靶向的恶性细胞中的MGMT表达水平更能直接反映MGMT蛋白活性，因此对治疗反应更具预测性。

Para_04

1. MGMT低表达和高表达样本之间的细胞结构比较显示，在细胞类型组成和恶性状态分布上没有显著差异。
2. 然而，MGMT表达水平与复发时MES样状态的丰度变化显著相关，这在我们的CARE数据集和参考文献14中均有体现。
3. 具体来说，与原发样本相比，MGMT低肿瘤复发时MES样状态的丰度较低，而MGMT高肿瘤则表现出相反的模式。
4. 基线特征以类似的方式发生变化——在MGMT低肿瘤中趋向于减少细胞外基质，这突显了细胞状态和BPs的协调进化。

Para_05

1. 为了进一步识别与MGMT高和MGMT低样本相关的纵向基因表达模式，我们比较了它们的恶性伪批量表达（每个样本中所有恶性细胞的基因表达汇总）在不同时间点的变化。
2. 这揭示了四组基因——那些在复发过程中持续上调或下调的基因（无论MGMT表达如何）（CONS-UP和CONS-DN），以及仅在MGMT低样本中上调的基因（low-UP）或在MGMT高样本中上调的基因（high-UP；扩展数据图5h）。
3. CONS-UP富含与突触后活动相关的基因（例如，GRIN1/GRIN2A/GRIN2B，GABRA1/GABRA4），这与复发中神经元-BP趋势增加一致（图3a），而CONS-DN则富含与氧化磷酸化相关的基因（例如，NDUFB4/NDUFB11，UQCRB/UQCRQ）和干细胞特性相关的基因（例如，CD24，DCX，SOX4/SOX11和SOX11；扩展数据图5i）。
4. 与这些MGMT无关的基因集相反，low-UP基因富集于与胶质发生相关的基因（例如，ACSL1，OLIG1/OLIG2和SOX8），而high-UP基因则富集于间质发育和细胞外基质重塑相关基因（例如，COL1A1/COL1A2，COL1A2，DCN，FN1和LOXL1；图4g）。
5. 为了验证我们的发现，我们在参考文献14的数据集中重复了此分析，并发现显著的相似性，特别是在high-UP特征中（Jaccard指数8.5%对比预期1.2%，P < 0.0001，自举法）。

Para_06

1. 总之，我们的单细胞分析表明，恶性细胞中MGMT的表达可以预测复发时间并影响纵向细胞轨迹——MGMT高（可能无反应）患者的复发中，MES样状态和与ECM相关的基因往往富集，而在MGMT低（可能有反应）患者中，复发与MES样状态的丧失以及与胶质发生相关的基因富集有关。
2. 我们的结果表明，之前未能将MGMT表达与治疗反应相关联的尝试可能是由于整体分析以及非恶性细胞中MGMT的表达而受到干扰。

Mutational signatures associated with recurrence trajectories

与复发轨迹相关的突变特征

Para_01

1. 治疗可能影响GBM的遗传进化，这得到了关于TMZ化疗后体细胞超突变和/或放疗后小缺失表型的报告支持。
2. 为了分析CARE队列中纵向突变模式，我们检查了治疗后的突变负担和突变特征。
3. 我们观察到突变负担（每兆碱基突变数）在复发时通常增加（P = 0.01，双侧配对Wilcoxon检验；扩展数据图6a），并识别出四个与治疗相关的高突变样本（复发时突变负担 > 10突变每兆碱基（Mut per Mb）；方法）。
4. 两个肿瘤对在两个时间点均表现出高突变负担，但没有治疗相关突变特征的证据（COSMIC单碱基替换（SBS）11；扩展数据图6b）。

Para_02

1. 为了研究复发时新的突变特征关联，我们检查了COSMIC v3 SBS特征的频率，包括我们的CARE队列（n = 46名患者）和GLASS队列（n = 138名患者）。
2. 在两个队列中，我们观察到SBS1（衰老/时钟样病因）显著减少，而SBS21（错配修复缺陷）增加（图5a和扩展数据图6c）。
3. 值得注意的是，突变归因于SBS21的患者并未倾向于获得错配修复基因的突变。
4. 虽然突变特征在时间上有所不同，但我们未观察到复发时富集的个体突变（调整后的Fisher精确检验P > 0.05）。
5. 在小插入/缺失变异中，我们观察到删除负担增加（图5b），但插入没有增加（扩展数据图5d），包括十次复发病例在无治疗相关高突变的情况下获得了小缺失表型（方法；图5b）。
6. 综上所述，我们的分析表明，部分肿瘤在治疗后会获得一致的遗传改变模式。

Fig. 5: Mutational signatures associated with recurrence trajectories.

图片

- 图片说明

◉ a，针对选定的SBS特征（SBS1，时钟样；SBS11，烷化剂；以及SBS21，DNA错配修复缺陷），在每个时间点上分配到COSMIC v3突变特征的突变相对比例的纵向差异。异常点显示出来。虚线连接同一患者在T1和T2的时间点值。P值表示双侧配对Wilcoxon秩和检验。◉ b，T1和T2之间小缺失负担的纵向差异，虚线连接患者值，红色表示被归类为小缺失负担增加组的患者。分析仅限于非高突变病例。P值反映双侧配对Wilcoxon秩和检验。◉ c，按纵向变化递减排序的缺氧相关恶性状态的瀑布图。患者由一个单独的条形图（x轴）表示。增加的小缺失负担或SBS21特征以黑色注释条形图下方展示。报告的P值表示对于有小缺失或SBS21增加的患者的缺氧恶性状态纵向变化的双侧配对Wilcoxon秩和检验。◉ d，来自GLASS联盟的批量RNA数据对缺氧snRNA恶性元程序特征进行了评分，y轴上展示了相对ssGSEA得分。这些面板显示了那些肿瘤有小缺失增加（左）和没有小缺失负担增加（右）的患者的缺氧特征的纵向变化。双侧配对Wilcoxon秩和检验的P值。◉ e，散点图表示标准化的遗传距离（批量DNA—私有突变的比例和拷贝数改变区域的差异）与标准化的转录距离（snRNA），蓝色线条表示线性回归线。显示了皮尔逊相关系数及其相关的P值（n = 38）。ssGSEA，单样本基因集富集分析；diff.，差异。

Para_03

1. 为了了解治疗后的基因改变是否会影响细胞状态的演变，我们研究了具有高突变、SBS21和小缺失表型的肿瘤的特定轨迹。
2. 对于高突变肿瘤，我们在复发时没有观察到一致的细胞表型，尽管该分析受到样本量的限制（n = 4）。
3. 相反，获得小缺失的患者富集于缺氧和MES样细胞状态（图5c和扩展数据图6e；配对Wilcoxon秩和检验，P < 0.05）。
4. 我们观察到具有增加的SBS21信号的肿瘤也有类似的效果，尽管在排除同时具有小缺失表型的样本后，这种效果较弱（图5c和扩展数据图6f）。
5. 复发特异性小缺失与缺氧状态之间的关联在GLASS9的批量样本分析中得到重现（图5d）。
6. 相反，MES样观察结果在GLASS队列中未得到验证（扩展数据图6g），这可能是因为MES样和非恶性特征之间存在相似性，使得在批量数据集中难以区分这些信号。
7. 缺氧细胞已被证明与表观遗传学异常和应激反应有关36，并因此与对放射治疗的抗性增加有关37。
8. 因此，缺氧区域的恶性细胞可能更有能力修复放疗相关的DNA损伤，导致复发时缺氧状态的扩增。

Para_04

1. 最后，我们检查了总体遗传距离（由纵向单核苷酸变异（SNV）/拷贝数变异（CNA）差异定义）与总体恶性转录距离之间的相关性（方法）。
2. 我们发现存在显著的相关性（皮尔逊相关系数 = 0.42，P = 0.009；图5e），在仅分析SNV或CNA时该相关性仍然存在（扩展数据图6h,i），这表明疾病进展期间肿瘤内遗传进化对转录进化有贡献。

Discussion

Para_01

1. GBM中的肿瘤间和肿瘤内异质性已被认为与治疗耐受和疾病复发有关，但特定的GBM细胞亚群如何对治疗产生反应仍知之甚少。
2. 在此，我们利用了来自59名患者的121个肿瘤的纵向匹配的多组学GBM单核RNA测序和批量DNA数据集，包含约43万个细胞核，以全面表征GBM在标准治疗后经历的组成、转录和遗传变化。
3. 我们发现不同患者的纵向轨迹存在广泛的变异（图6）。
4. 尽管许多轨迹不可预测，但我们识别出与特定轨迹相关的患者特征——MGMT甲基化和小片段缺失负担增加——并扩展了我们对复发的理解。

Fig. 6: Longitudinal GBM cell state trajectories.

图片

- 图片说明

◉ 对GBM CARE数据集的分析揭示了治疗后GBM进化的关键模式。◉ （1）最常见的观察结果是在细胞类型组成层，许多样本在复发时表现出恶性细胞比例降低。◉ （2）在恶性状态层面，许多肿瘤以不可预测的轨迹复发，这可能反映了患者特异性的进化（顶部）。◉ 一小部分肿瘤以更可预测的轨迹复发，例如MGMT-MET与MGMT-UM原发性肿瘤相比，MES样比例减少（中间），或在治疗后小缺失负荷增加的样本中缺氧比例增加（底部）。◉ 分为恶性和非恶性比例的列显示了每种肿瘤的相对比例。◉ 为清晰起见，仅展示了选定的非恶性和恶性细胞状态。◉ Hyp，缺氧。

Para_02

1. 先前的研究强调了GBM进化过程中的特定转录轨迹。
2. 最近的GLASS联盟研究强调了向神经元、MES和增殖特征的三种轨迹9，其他研究也强调了前神经元向MES表型的转变11,14,38。
3. 虽然我们同样识别出这些复发相关模式，但我们的分析表明它们的整体频率是有限的，且这些轨迹仅在特定患者亚组中富集，而原发性和复发性样本之间细胞状态的总体比较显示整体相似性。
4. 特别是，在治疗后向MES表型的转变并未在我们的数据集中作为全局富集趋势观察到，仅在特定患者亚组中观察到。
5. 这种与先前研究的不一致可能是由于以下三个因素：(1) 技术差异（例如，批量测序与单细胞测序），(2) 队列/采样（例如，单中心采样，组织消化与未消化），以及(3) 分析方法。
6. 此外，还可能存在由于未解决的额外治疗导致的偏倚（例如贝伐珠单抗、类固醇等）。

Para_03

1. 我们的数据集包含来自六个国家七个医疗中心的肿瘤样本，这降低了因机构特定程序导致结果偏向某种复发轨迹的可能性。
2. 此外，我们数据集中的样本由三个不同的实验室处理和分析，这降低了因实验室特定批次效应导致偏差的可能性。
3. 我们实验室用于处理和分析肿瘤样本的单细胞核测序协议是消化独立的，并针对脑肿瘤进行了优化，捕获的转录本数量比以往的单细胞GBM研究更多（扩展数据图2e）。
4. 这很可能使肿瘤组织的表示更加准确，这从我们数据集中TME细胞类型的更高多样性以及发现新的恶性细胞状态中可以得到体现。

Para_04

1. 神经元通路在复发性GBM中被一致发现是上调的。
2. 然而，需要注意的是，这仅反映了每项研究中的部分患者，并且我们观察到从原发性到匹配的复发样本存在广泛的轨迹差异。
3. 利用我们数据集的分辨率和规模，我们表明这种现象与恶性细胞比例的降低以及肿瘤微环境中神经元比例的增加相关；只有在少数病例中，这伴随着NEU样恶性状态和肿瘤的神经元基线表达特征的富集。
4. 这些结果表明，在治疗后，残留的GBM细胞与脑实质的整合可能增强。

Para_05

1. 除了定义的子集外，我们的研究提出了一个具有挑战性的观察结果，即大多数GBMs的纵向细胞轨迹是多样且不可预测的。
2. 我们推测这种多样性至少反映了以下三种现象：首先，鉴于原发肿瘤中的遗传和表观遗传异质性，可以预期具有不同潜在生物学特性的肿瘤也会有不同的复发轨迹。
3. 其次，即使在单一GBM中，也存在广泛的区域异质性；由于原发和复发样本是从不同位置获取的，它们的差异可能反映了这种区域异质性，并因此依赖于取样。
4. 第三，患者在治疗史和对治疗的反应方面存在显著差异，这可能导致不同的复发轨迹。
5. 许多GBM患者接受的额外治疗方式多种多样，超出了标准治疗方案，这是本研究的一个局限。
6. 此外，即使是标准治疗也会引发不同的患者反应，我们识别出以下两种情况，其中治疗反应与复发轨迹相关：(1) MGMT-MET（MGMT低表达）样本（更可能对烷基化剂产生反应）与复发时MES样状态的耗竭相关；(2) 小片段缺失表型，之前与放疗相关，与复发时缺氧状态的增加相关。

Para_06

1. 总之，我们对GBM复发进行了单细胞分辨率的全面分析。
2. 我们的分析表明，总体而言，复发性GBM在细胞状态、组成和其他GBM生态系统特征方面通常与原发性GBM相似。
3. 因此，观察到广泛的复发轨迹，突显了GBM进化中的广泛多样性。
4. 然而，有特定治疗反应证据的肿瘤亚群倾向于偏向特定的转录轨迹。
5. 这些发现提供了比以往观察更详细的见解，并为更精确地分层复发患者提供了机会，可能有助于制定具体的治疗策略。

Methods

Experimental method details

实验方法细节

Human participants and ethical approval

人类受试者和伦理批准

Para_01

1. 根据世界卫生组织2021年分类，被诊断为‘胶质母细胞瘤，IDH野生型’的冷冻GBM标本在七家机构收集，如下所示。
2. 收集和基因组分析已获得各机构的机构审查委员会批准，所有患者均相应提供了知情同意。
3. 这些肿瘤来自MD安德森癌症中心，在MD安德森癌症中心机构审查委员会协议2012-0441下获得批准。
4. 杜克大学队列（n = 12）——这些肿瘤来自杜克大学医院的普雷斯顿·罗伯特·蒂施脑肿瘤中心生物库。该研究在杜克大学机构审查委员会协议Pro0007434下获得批准。
5. 东京大学队列（n = 30）——这些肿瘤来自东京大学医院神经外科。该研究在东京大学医院机构审查委员会协议G10028下获得批准。
6. 皮蒂耶-萨尔佩特里埃医院队列（n = 28）——这些肿瘤来自皮蒂耶-萨尔佩特里埃医院。该研究在医院皮蒂耶-萨尔佩特里埃的Onconeurotek脑肿瘤库认证96-900下获得批准。
7. 圣迈克尔医院队列（n = 14）——这些肿瘤来自多伦多统一健康圣迈克尔医院神经外科部。该研究在圣迈克尔医院、多伦多统一健康的研究伦理委员会协议REB 13-141下获得批准，所有患者均相应提供了书面知情同意。
8. 首尔国立大学（SNU）队列（n = 20）——这些肿瘤来自SNU医院神经外科部。该研究在SNU医院机构审查委员会批准（批准号H-2004-049-1116）下进行。
9. NORLUX神经肿瘤实验室队列（n = 13）——这些肿瘤来自卢森堡中心医院（CHL，神经外科部）。该研究在卢森堡国家研究伦理委员会协议201201/06下获得批准。
10. 队列被加入到达纳-法伯/哈佛癌症中心10-417机构审查委员会协议中。
11. 患者包括男性和女性。
12. 队列的临床信息总结在补充表1中。

Statistics and reproducibility

统计学与可重复性

Para_01

1. 没有使用统计方法来确定样本量。
2. 当样本中恶性细胞的数量在至少一个时间点较低时，数据会被排除在某些分析之外。
3. 在本文的图表中（图2–5和扩展数据图1–3、5和6），箱形图反映了以下汇总统计信息：(1) 分割箱体的线代表中位数；(2) 箱体的下边缘和上边缘分别对应第一四分位数和第三四分位数（第25和75百分位数）；(3) 上须线从边缘延伸到不超过边缘1.5倍四分位距（IQR）的最大值；(4) 下须线从边缘延伸到不超过边缘1.5倍四分位距的最小值；(5) 如果图表中没有显示代表实际数据的抖动点，则超出须线末端的数据会单独作为异常点绘制。
4. 本研究中描述的统计分析是使用R版本4.0.1及以上进行的。

Nuclei isolation from frozen tissue

从冷冻组织中分离细胞核

Para_01

1. 方案1（杜克大学、MD安德森癌症中心、东京大学医院、圣迈克尔医院和皮蒂耶-萨勒皮特里埃医院的队列）—从冷冻肿瘤组织中分离出细胞核，方法如下：组织解冻后，在盐-Tris缓冲液（10 mM Tris-HCl（pH 7.5），146 mM NaCl，1 mM CaCl2和21 mM MgCl2）中进行机械解离，加入0.49% CHAPS（3-((3-胆酰胺丙基)二甲基铵)-1-丙磺酸盐）（Millipore，28300）。单细胞核悬液通过40 μm筛网过滤，以500g离心5分钟，然后在补充了0.01% BSA（NEB，B9000S）的盐-Tris缓冲液中重悬。细胞核悬液通过显微镜检查，使用血球计数板计数，并用于10x Genomics工作流程。

Para_02

1. 方案2（NORLUX和SNU队列的实验室3）—组织样本在核EZ裂解缓冲液（Millipore Sigma，NUC101）中通过Dounce匀浆进行机械解离。溶液在冰上孵育5分钟，并在孵育过程中搅拌一到两次。
2. 单核悬浮液通过70微米滤网过滤，并在4°C下以500g离心5分钟，重新悬浮在核EZ裂解缓冲液中并在冰上孵育5分钟。
3. 溶液在4°C下以500g离心5分钟，重新悬浮在1% BSA/0.2 U μl RNase抑制剂/PBS缓冲液中（三次）。
4. 最终重悬时，向缓冲液中加入DAPI，溶液通过40微米滤网过滤，使用Countess II自动细胞计数器（Thermo Fisher Scientific）进行细胞计数，并将核转入10x Genomics工作流程。

10x Genomics for single-nucleus sequencing

10x基因组学用于单核测序

Para_01

1. 根据制造商的说明书使用了10x Chromium Single Cell 3′ Reagent Kit v3（10x Genomics，PN1000128）。
2. 简而言之，核被加载到Chromium Chip（10x Genomics，PN12000120）上，目标核回收数为6,000–8,000个，并在Chromium Controller中进行处理。
3. 单个核被分配到乳剂中的凝胶珠（GEMs）中，随后通过条形码进行RNA逆转录。
4. 通过破裂GEMs并合并条形码部分、cDNA扩增、片段化以及样品索引和接头的连接来构建文库，并在NextSeq500或NovaSeq（Illumina）上进行测序。
5. 对于NORLUX和SNU队列，核按照制造商的说明书，以目标细胞回收数为6,000个核的方式加载到Chromium芯片上，用于单细胞多组学ATAC和基因表达。
6. 10x多组学的基因表达化学反应与10x v3 3′ snRNA-seq的化学反应相同。

WES

WES

Para_01

1. 杜克大学、MDACC、东京大学医院和圣迈克尔医院的样本进行了全外显子组测序，方法如下：使用DNeasy Blood & Tissue试剂盒（Qiagen，69504）从每个冷冻肿瘤样本和对应的血液样本中提取DNA。
2. 基因组DNA（100–250 ng）通过超声波破碎仪（Covaris）进行声波剪切，目标片段大小为150 bp。
3. 使用KAPA HyperPrep试剂盒（KAPA Biosystems，KK8504）进行文库制备，随后使用AMPure XP磁珠（Beckman Coulter，A63882）进行酶促纯化。
4. 使用定制的外显子捕获探针（由Twist Biosciences根据Broad Institute的要求制造）进行外显子捕获。
5. 捕获后的文库在NovaSeq 6000（Illumina）上进行150碱基对配对末端测序。
6. 对于Pitié-Salpêtrière医院的样本，在DNA经LE220超声波破碎仪（Covaris）剪切并进行片段大小选择后，使用Twist Human Core Exome试剂盒（Twist Bioscience）按照制造商的说明书进行文库制备和捕获。
7. 测序在NovaSeq 6000（Illumina）上进行，采用200碱基对配对末端测序。

WGS

WGS

Para_01

1. 为NORLUX神经肿瘤实验室的一组冷冻样本收集了新的全基因组DNA测序数据。
2. 简要来说，对于有足够肿瘤组织的样本以及在可用时匹配的正常血液样本，使用AllPrep DNA/RNA Minikit（Qiagen, 80204）提取DNA。
3. 需要注意的是，用于批量DNA分离的组织是用于单核数据生成的组织旁边的组织。
4. 简要来说，使用LE220超声仪（Covaris）将DNA剪切成400 bp，并使用AMPure XP磁珠（Beckman Coulter, A63882）进行大小选择。
5. 全基因组文库使用NovaSeq 6000（Illumina）进行制备和150碱基对配对末端测序。
6. SNU队列的WGS数据以相同的方式制备，并且之前在参考文献9中报道过。
7. SNU队列的数据可在Synapse上获得（https://www.synapse.org/glass）。

Somatic variant detection and analysis

体细胞变异检测与分析

Para_01

1. DNA测序比对、指纹分析、体细胞变异检测（Mutect2）和拷贝数分割是根据基因组分析工具包（GATK）的最佳实践使用GATK 4.0.10.1进行的，如之前所述。
2. 简要来说，全外显子组和全基因组读段是使用BWA MEM 0.7.17与b37基因组（human_g1k_v37_decoy）进行比对的。
3. 使用‘CrosscheckFingerprints’（Picard）进行DNA指纹分析，确认所有属于同一患者的样本均来自同一个体，表明本研究中没有样本错配。
4. 使用Mutect2（v4.1.0.0）检测肿瘤样本中的体细胞突变，并且需要匹配的正常血液样本。
5. 为每个测序批次构建了正常样本的面板，以考虑DNA文库制备的差异（例如全基因组与全外显子组），并用于过滤常见的人工假阳性变异和胚系变异。
6. 对于没有匹配正常血液的患者肿瘤样本，进行了拷贝数变异分析，并根据GATK最佳实践使用正常样本参考面板进行肿瘤单样本SNV检测。
7. 肿瘤单样本体细胞变异用于识别高突变肿瘤。
8. 为了识别小缺失突变负担增加的样本，应用了以下阈值：复发特异性小缺失突变负担需要大于每兆碱基0.2个突变，并且在比较复发和原发肿瘤的所有小缺失变异时，每兆碱基的增加需要大于0.1个突变。
9. 突变特征估计是根据先前发表的方法，使用COSMIC v3特征对GLASS和CARE数据集进行的。
10. 简要来说，所有全基因组测序样本被用来生成突变矩阵，并使用SigProfilerExtractor（v.1.1.4）提取新的突变特征，最大特征数设为10，NMF重复次数设为100。
11. 然后将这些新特征去卷积为COSMIC v3.3特征。
12. 最后，使用R包Palimpsest（v.2.0.0）将每个样本的突变分配到去卷积后的COSMIC特征中最可能的特征。
13. 获得SBS21的样本被定义为在T2时SBS21变化处于前四分位数的样本。
14. 对于有匹配血液的样本，根据每兆碱基10个突变的截止值确定高突变，并进一步根据COSMIC突变特征分为新生或治疗相关。
15. 对于有匹配正常血液的样本，治疗相关的高突变通过突变负担的纵向增加（即从<10个突变/兆碱基到>10个突变/兆碱基）和SBS11信号（烷化剂相关特征）进行分类。
16. 对于没有匹配血液的样本，一个复发样本的突变负担相比初始肿瘤显著增加（18.3倍）并伴有SBS11特征，提示为治疗相关的高突变。
17. 同一患者时间分离样本之间的遗传距离定义为私有突变（SNV）的总DNA比例以及两个样本之间CNA区段的差异。

Single-cell/single-nucleus data processing of human glioma samples

人脑胶质瘤样本的单细胞/单核数据分析

错误！！！ - 待补充

Clustering and cell-type annotation

聚类和细胞类型注释

Para_01

1. 为了便于聚类、细胞类型注释和基于微滴的数据评分，我们使用了 R 包 Seurat（v4.0.4）。
2. 唯一分子标识符计数数据被归一化和缩放（使用 NormalizeData 和 ScaleData 函数），然后进行聚类（通过主成分分析将数据投影到低维空间，然后运行 FindNeighbors 和 FindClusters 函数）。
3. 肿瘤相关巨噬细胞、T 细胞、少突胶质细胞和内皮细胞的初步分类是基于已知标记基因的表达。
4. 而其余细胞则被注释为推测性的恶性细胞。
5. 在将细胞分类为恶性或非恶性之后（参见章节‘从基因表达数据推断 CNA’和‘将细胞分类为恶性或非恶性细胞类型’），将非恶性细胞进行评分（使用 AddModuleScore 函数）以检测已知的细胞类型特征，并使用‘将细胞分配到状态’章节中解释的方法分配细胞类型。
6. 最后，使用 R 包 scDblFinder 过滤掉异型双细胞。

Inferring CNAs from gene expression data

从基因表达数据推断拷贝数变异

Para_01

1. CNAs 之前文所述方法46,47,48,49进行估计，使用 R 包 infercna 中的 infercna 函数（该方法在参考文献50中首次提出，可在 https://github.com/jlaffy/infercna 获取）并采用默认参数。
2. 简而言之，该算法根据每个样本中基因的染色体位置对分析的基因进行排序，并在每个染色体中使用包含 100 个基因的滑动窗口对相对表达值进行移动平均。
3. 被分类为非恶性细胞（少突胶质细胞、巨噬细胞和内皮细胞）的评分定义了正常核型的基线，它们的平均 CNA 值用于将所有细胞的值居中。

Classification of cells into malignant or nonmalignant cell types

将细胞分类为恶性或非恶性细胞类型

Para_01

1. 在使用单核液滴测序分类的细胞中，我们通常遇到较低的CNA信号和与以下因素相关的连续信号分布：(1) 相对于单细胞数据而言数据质量较低，以及(2) 在假定的恶性群体中包含了非恶性细胞类型，如ACs和神经元，因为这些细胞类型由于与恶性细胞的转录相似性，无法事先被分类为非恶性。

Para_02

1. 因此，我们为基于液滴的单核测序数据使用了多步骤分类方案。

Step 1: cell classification by detection of copy-number events

第一步：通过检测拷贝数事件进行细胞分类

错误！！！ - 待补充

Step 2: exclude questionable cells based on CNA correlation

步骤2：根据CNA相关性排除可疑细胞

Para_01

1. 我们将CNA相关性定义为每个细胞与其来源肿瘤样本的平均CNA图谱之间的皮尔逊相关性。
2. 我们仅将计算限制在检测到‘真实’事件的染色体以及来自其他染色体的约1,000个基因的对照面板上。
3. 基于分类为恶性细胞的CNA相关性分布，我们将恶性性的下限阈值定义为0.5。
4. 为了确定非恶性的上限阈值，我们从snRNA-seq数据中进行了SNV调用，以估计算法误分类率。
5. 简而言之，对于每个具有可用WES数据的样本，我们使用Vartrix工具（https://github.com/10XGenomics/vartrix）从FASTQ文件中对每个细胞进行SNV调用。
6. 误分类的细胞被定义为携带恶性突变但在算法第一步中被分类为非恶性的细胞。
7. 经过成本效益分析（误分类率与排除细胞百分比），我们将非恶性的上限阈值设为0.35。
8. 最后，将CNAcor < 0.5的‘恶性’细胞或CNAcor > 0.35的‘非恶性’细胞重新分类为‘未解决’，并从下游分析中排除。

Step 3: refine classification by setting confidence levels

步骤3：通过设置置信度级别来优化分类

Para_01

1. 为了细化分类，我们定义了置信度等级。置信度等级1包括被分类为恶性且携带恶性SNV的细胞，检测到至少50%预期CNA事件的细胞，包括染色体7扩增或染色体10缺失，以及被分类为非恶性的CNAcor<0.25的细胞。
2. 置信度等级2包括被分类为恶性且检测到至少50%预期CNA事件但不包括染色体7扩增或染色体10缺失的细胞，被分类为恶性但检测到少于50%预期CNA事件的细胞，包括染色体7扩增或染色体10缺失，以及被分类为非恶性的CNAcor<0.5的细胞。
3. 其余被分类为恶性或非恶性的细胞被分配为置信度等级3。

Para_02

1. 在本研究中提出的分析中，我们仅使用了置信度等级1和2的细胞。

Deriving metaprograms from gene expression data

从基因表达数据中推导元程序

Para_01

1. 为了捕捉同一细胞类型之间的异质性，我们利用了NMF51。
2. NMF是在将负值设为零后，对每个样本的相对表达值独立进行的。
3. NMF算法需要事先定义‘k’参数，反映数据中预期的潜在特征数量。
4. 由于‘k’在样本之间有所不同且大多未知，我们使用不同的值（3–10）对每个样本运行NMF算法，从而为每个样本生成52个程序。
5. 每个NMF程序通过基于NMF系数的前50个基因进行总结。
6. 从NMF程序推导出元程序的详细方法在参考文献24中有所描述，并在此简要说明。
7. 该方法首先过滤掉不稳健的NMF程序（在它们产生的样本内或多个样本中不重复）或在样本内冗余的程序（即与其他NMF程序高度重叠）。
8. 随后，根据Jaccard相似度对稳健的NMF程序进行聚类，使用自定义的聚类算法，该算法迭代考虑稳健NMF程序之间的重叠程度，并将高度重叠的程序组合成一个簇。
9. 每个簇中前50个最常出现的基因定义了一个元程序。

Para_02

1. 总体而言，该算法生成了16个元程序，这些元程序基于功能富集分析进行注释（使用基因本体（GO）术语、mSigDB标志性基因集合以及从正常大脑发育数据集中衍生的基因集合）。
2. MP 16包含来自各种细胞类型的基因，导致一个几乎适用于所有细胞的元程序，因此没有反映任何异质性，并因此被排除在分析之外。
3. 此相同的过程对每种非恶性细胞类型分别重复进行。

Assignment of cells to states

细胞到状态的分配

Para_01

1. 恶性细胞通过使用 Seurat 函数 AddModuleScore 对 NMF 代谢程序进行评分。
2. 为了便于细胞分类，我们通过每次采样 5,000 个细胞并打乱每个基因的表达值，生成了 20 个随机表达矩阵。
3. 然后我们对每个随机矩阵对 NMF 代谢程序进行评分，从而为每个表达程序得到 100,000 个正态分布的分数。
4. 这些作为细胞状态分类的零分布。
5. 对于每个原始细胞，我们使用之前生成的零分布统计信息，通过 z 检验（R 的 pnorm 函数）计算每个表达程序的 P 值。
6. 我们使用 BH 方法对所有 P 值进行多重检验校正。
7. 如果某个状态的校正后 P 值小于 0.05，则将该细胞分类到特定状态。
8. 那些在多个状态中达到校正后 P < 0.05 的细胞被分配到得分最高的状态。
9. 那些在任何状态中都没有达到校正后 P < 0.05 的细胞被分配到"未解决"状态。

Para_02

1. 如果细胞在细胞周期元程序中达到调整后的P值小于0.05，则被分配为‘循环’状态，否则为‘非循环’状态。

Hybrids detection and classification

混合检测与分类

Para_01

1. 我们定义杂合细胞为至少在两个细胞状态中得分显著（即调整后的P值小于0.05），且最高两个显著得分之间的差异小于随机预期的细胞。
2. 为了量化随机预期的得分差异，我们利用了状态分类算法生成的空分布，并计算了每个人工细胞中最高两个状态之间的差异。
3. 该得分差异分布的95百分位数（等于0.24的差异）被设定为判断状态是单一还是杂合的阈值。

Para_02

1. 为了估计技术杂合体的预期频率，我们首先计算了每对状态（A，B）的杂合体预期频率，定义为Exp(A,B)=Freq(A)×Freq(B)，以及O/E比率，定义为log2(Obs(A,B)/Exp(A,B))。
2. 由于这是由于技术效应导致的预期频率的高估，我们通过平均那些由随机检测到的杂合体对（即O/E < 0）的O/E比率来估计这个技术因素。
3. 随后，我们通过将预期频率乘以技术因子，计算出每对杂合体的技术杂合体的预期频率，并将O/tE比率定义为log2(Obs(A,B)/(TF×Exp(A,B)))。
4. 最后，我们使用单侧Fisher检验过滤掉O/tE比率较大的不显著杂合体对（归因于数量较小）。
5. 杂合谱系模型是使用R包igraph生成的，其中细胞状态作为节点，连接形成显著杂合对的状态作为边。
6. 我们没有将"循环"状态纳入此分析，因为我们不认为"循环"是一个独立的状态，而是认为它是在神经发育身份（如AC样、OPC样等）之上的额外特征。

Baseline profile derivation

基线配置文件推导

Para_01

1. 每个肿瘤样本通过平均标准化的唯一分子标识符计数并进行log2转换，分解成七个伪批量表达谱（每个常见状态一个——AC-like、MES-like、缺氧、GPC-like、OPC-like、NPC-like和NEU-like——每个状态至少有25个细胞被分类到该状态）。
2. 如果基因在所有伪批量中的平均log2表达值至少为1，并且其中位方差（在每个状态内单独计算）至少为2.5，则将其包含在内。
3. 总体而言，有1,005个基因通过了这些过滤条件。
4. 接下来，对每个状态内的伪批量进行分析，以主成分分析去除状态特异性信号。
5. 然后从每个状态中推导出六个基因程序（前三个主成分的顶部和底部50个基因），计算每对程序之间的Jaccard相似性指数，并使用层次聚类对相似性矩阵进行聚类。
6. 这揭示了五个不同的聚类，我们通过包含每个聚类中至少25%的程序中重复出现的基因（最低三次）推导出五个共识特征。
7. 通过人工检查和GO富集分析，我们对共识特征进行了注释。
8. 由于长度较短或高度怀疑反映技术性偏差，有两个特征被排除。
9. 为了排除环境RNA可能驱动基础表达谱的可能性，我们使用R包SoupX估计了每个样本的环境RNA污染。
10. 不同基础表达谱之间的估计污染水平没有差异，且每个样本中污染最严重的基因与元程序或基础表达谱特征之间也没有显著重叠（由于缺乏显著性未显示数据）。

Assigning tumors to composition groups

将肿瘤分配到组合组

Para_01

1. 为了衡量肿瘤组成的差异，我们为每个肿瘤样本生成了一个组成向量，反映了肿瘤样本中每种细胞类型的比例。
2. 这使得可以将肿瘤分为以下三个主要组成组：HMF（恶性细胞百分比>75%）、IMF（恶性细胞百分比为50–75%）和LMF（恶性细胞<50%），其中LMF组可以根据主导的TME细胞类型进一步细分（TAM/少突胶质细胞/胶质神经元百分比>40%，或在没有主导TME细胞类型的情况下为混合型）。
3. 同样，我们为恶性细胞状态定义了比例向量，并根据主导细胞状态（至少25%的细胞被分配到特定状态）将肿瘤分为六个组：AC、MES/缺氧、GPC、OPC/NPC、神经元和混合型（在没有主导状态的情况下）。
4. 肿瘤根据基线特征的最高得分进行分类，如果最高得分小于0.25，则归类为混合型。

Multilayer group definition

多层组定义

Para_01

1. 每对特征（A，B）之间的关联是通过二项式检验计算的，将该对出现的次数作为成功次数，肿瘤样本的数量作为实验次数，特征A的频率乘以特征B的频率作为成功的期望概率（假设特征是独立的）。我们然后通过将每个特征定义为节点，并在观察到两个特征之间存在统计学显著关联（未经多重检验校正）时用边连接节点，生成了一个无向图。

Measuring transcriptional distance between matched pairs

测量匹配对之间的转录距离

Para_01

1. 为了衡量匹配对之间的转录距离，我们从每对中生成了状态控制的伪批量表达谱。
2. 状态控制是通过将两个样本都下采样到每个状态的25个细胞来实现的，前提是两个样本都至少包含25个该状态的细胞，从而得到具有相同细胞数且在不同细胞状态的基因表达谱上平衡的样本。
3. 随后，通过在选定的细胞中对每个基因进行平均计算出伪批量表达谱，从而得到两个伪批量表达谱。
4. 最后，计算这两个向量之间的欧几里得距离，并作为转录距离的度量。

Survival analysis

生存分析

Para_01

1. 使用 R 包 survival v.3.3-1 和 survminer v.0.4.9 进行了生存分析。
2. Kaplan–Meier 估计使用 survfit 函数计算，并使用 ggsurvplot 函数（survminer 包）进行绘制。
3. 生存曲线之间差异的统计显著性通过 log-rank 检验计算（在函数 survdiff 中实现）。
4. Cox 回归分析和森林图使用相同的包（coxph 和 ggforest 函数）进行并绘制。

Differentially expressed genes between MGMT-low and MGMT-high pairs

MGMT低表达与MGMT高表达配对中的差异表达基因

Para_01

1. 为了计算MGMT-low和MGMT-high组之间的差异表达基因，我们首先计算了每组在时间点上的相对基因表达（即纵向表达差异），该差异对于每个基因i定义为该组所有T2样本中基因的平均表达减去该组所有T1样本中基因的平均表达。
2. 接下来，我们计算了这两个向量之间的欧几里得距离。
3. MGMT-low基因被定义为在组间距离最大的前100个基因，并且在MGMT-low组中具有正向的纵向表达差异，在MGMT-high组中具有负向的纵向表达差异。
4. MGMT-high基因被定义为在组间距离最大的前100个基因，并且在MGMT-high组中具有正向的纵向表达差异，在MGMT-low组中具有负向的纵向表达差异。

Reporting summary

报告摘要

Para_01

1. 有关研究设计的更多信息可在与本文相关的Nature Portfolio报告摘要中找到。

Data availability

Para_01

1. 本研究生成的处理后的snRNA-seq数据可在Gene Expression Omnibus中通过访问编号GSE274546（10x Genomics）获得。
2. 所有去标识化的体细胞突变调用、拷贝数变异调用和基因组分析表均可通过Synapse（https://www.synapse.org/care_glioblastoma）获取。
3. 原始测序数据可根据Data Use Oversight System（DUOS）的同意书条款，在https://duos.broadinstitute.org/下以有限的方式获得，相关ID为：DUOS-000475、DUOS-000476、DUOS-000477、DUOS-000478、DUOS-000479和DUOS-000480。

Code availability

Para_01

1. 本研究使用的分析脚本可在 GitHub 上获取（https://github.com/dravishays/GBM-CARE-WT）和 Zenodo52（https://doi.org/10.5281/zenodo.14966015）。处理 DNA 测序数据的脚本可访问 https://github.com/Kcjohnson/care-glass。