bulk and sc-RNA seq剖析HCC中肿瘤异质性，揭

Integrating bulk and single-cell RNA sequencing reveals cellular heterogeneity and immune infiltration in hepatocellular carcinoma

整合批量和单细胞RNA测序揭示肝细胞癌的细胞异质性和免疫浸润情况

发表期刊：Mol Oncol

发表日期：2022 Jun

影响因子：7.449

DOI:  10.1002/1878-0261.13190

一、研究背景

        肝细胞癌（HCC）是原发性肝癌中最常见的形式， HCC的特点是肿瘤间和肿瘤内的异质性。揭示HCC异质性的分子机制对于靶向治疗的发展至关重要。到目前为止，已经有几种药物被批准用于治疗HCC，但结果并不令人满意。

        免疫治疗在HCC中的效果是通过整体的免疫浸润和丰富的免疫亚群的共同出现来调节的。肿瘤微环境（TME）中复杂的细胞组成和特征突出了癌症免疫抑制的多种非冗余机制的存在。为了提高免疫疗法的疗效，应该进一步研究不同患者亚型中免疫细胞的异质性，并确定适合特定免疫疗法的患者。

二、材料与方法

1、数据来源

1） TCGA获得了353个样本

2） ICGC获得了232个样本用于后续分析

3） 单细胞数据集：GSE149614、GSE115469、GSE98638、GSE146409

4） BATF敲除RNA数据：GSE149197

2、分析流程

1） HCC亚类的鉴定：选择平均绝对偏差>1的基因进行NMF聚类

2） 多组学数据的采集和处理：使用'GSVA'r包，单样本基因组富集分析（GSEA，ssGSEA）评估标志基因的富集分数,标志性基因集来自GSEA软件中的MSigDB数据库

3） 差异表达基因分析：使用 "Limma "软件包来进行差异表达基因（DEG）分析；考虑到S1和S3的高免疫细胞浸润，以及S2的高肿瘤纯度，将S1和S3的差异表达基因与ImmPort的免疫基因重叠，排除S2的免疫基因；对每个聚类进行差异分析，与另外两个聚类的两个聚类进行比较

4） 免疫浸润和肿瘤纯度的估计：CIBERSORT、Estimate

5） 单细胞RNAseq数据处理

6） 富集得分：计算富集分数以评估每个簇中的亚型分布；计算每个簇中每个亚型的频率，将gP除以集群频率（集群中的细胞数除以总细胞数），得到每个集群中每个亚型的富集得分

7） 共表达网络的构建：WGCNA

8） SCENIC分析：使用R软件包SCENIC来分析细胞亚型中转录组因子的富集情况；使用AUCell评估每个细胞中每个调节子的活性

9） 免疫抑制评分、肝脏评分和活化T细胞评分：采用压抑T标志基因以及肝脏标志基因的表达来进一步计算免疫抑制评分和肝脏评分；肝脏评分是按照Kim等人的24个肝脏标记基因的平均表达量计算的；免疫抑制评分和激活的T细胞评分是根据Guo等人的35个已知被抑制的标记基因和28个激活的GZMK-CD8基因定义的

三、实验结果

01 - 非负矩阵分解确定了HCC的三种亚型

        首先，从TCGA获得353名HCC患者的转录组和体细胞突变信息。利用NMF算法的共识聚类分析，确定了三个不同的修饰模式集群，包括集群S1的120例，集群S2的144例和集群S3的89例。共识矩阵的热图显示出明显边界，表明聚类结果的准确性和稳健性（图1B）。还在ICGC队列中验证了聚类结果，并使用SubMap进行了复制分析。发现在三个亚型中，S3的患者预后最差（图1D）。同时，分析了这三个亚型的临床指标的异质性，发现S3患者的肿瘤分期明显较高，这可能部分解释了该亚型的预后不良。

图1 用NMF共识聚类法识别HCC亚类

02 - HCC的肿瘤间TME异质性

        接下来，用ESTIMATE算法分析了三个亚型的免疫浸润和肿瘤纯度的异质性。结果显示，S2的免疫和基质得分明显低于S1和S3，而S2的肿瘤纯度得分最高（图2A，图S2A）。由于S3在三种亚型中预后最差，用CIBERSORT进一步描述了他们在22种免疫相关细胞类型中的免疫学状况。结果显示，与S3相比，S1具有较高的活化NK细胞（aNK）、CD8+T细胞、M1细胞和CD4记忆静止T细胞的丰度，但Treg和M0细胞的丰度较低。S2的M2细胞的丰度较低（图2B）。特别是，T细胞标记基因，如CD3E，在S1和S3中高度表达（图2C）。

        作者进一步调查了三种亚型中免疫抑制标志基因的概况，并计算了免疫抑制分数。结果表明，免疫抑制评分（图S2B）和标记基因如VEGFA、CTLA4、HAVCR2和TIGIT在S3中高表达（图2C，图S2D）。高比率的免疫抑制性T细胞可能与预后不良有关。此外，S3患者的巨噬细胞标记基因CD68和EMT标记基因如MMP2和MMP9的水平很高，表明肿瘤微环境中的这些细胞可能在肿瘤进展中发挥重要作用（图S2D，E）。相比之下，S1患者的预后结果最好，其激活的T细胞标志物如CD3E、PRF1和GZMK的水平高于S2，T细胞的标志物基因如PDCD1和HAVCR2也耗尽了，而且免疫抑制的分数也低于S3（图2C，图S2C，D）。

图S2 三种HCC亚型的免疫学和突变相关性的肿瘤间异质性

        作者还通过GSVA对三个亚型的不同基因进行了GO_BP和HALLMARKER通路的富集。结果显示，免疫相关途径如IL6-STAT3和胸腺T细胞选择，在S1和S3中富集，而CD4激活、炎症反应和T细胞分化途径在S1中显示出较高的富集分数。此外，在S3中发现B细胞凋亡和WNT途径的富集分数较高，而在S2中发现脂肪酸、脂质代谢和氧化磷酸化的富集分数较高（图2D）。综上所述，S2因免疫浸润率较低而表现出 "冷肿瘤 "的特征，S1则表现出"热肿瘤 "的特征，而S3则因高表达的免疫抑制基因而表现出 "免疫抑制的肿瘤"。

图2 三种HCC亚型中的免疫学肿瘤间异质性的调查

        作者分析了HCC亚型中体细胞和拷贝数变异（CNV）突变的频率差异和不同亚型的特异性突变特征（图2A）。具体来说，S2的CTNNB1和ARID1A的突变频率明显高于S1或S3。与S1和S2相比，S3表现出更高的TP53和BAP1的频率。此外发现CTNNB1-非突变组的免疫评分和CD3D和CTLA4的表达较高。

        此外，三个亚型中的CNV谱存在明显的异质性。S1有最多的扩增变异样本，而S2有最多的删除样本（图3B）。在CNV变异区域，S1患者主要在1q、5p、8q、6p和11q等区域被扩增，S2患者在1q、11q、1q和2q等区域被扩增，S3患者在S2、8q和13q等区域被扩增（图3C，图S4B）。值得注意的是，在S3中高表达的YEATS4和VIMP在S1中被删除，但在S3中被放大（图3D，图S4A）。相比之下，在S1中高表达的CYFIP2和ABLIM3在S1中同时被扩增，但在S3中被删除（图3D，图S4A）。综合来看，BAP1的最高突变、YEATS4和VIMP的扩增以及CYFIP2和ABLIM3的缺失，可能会诱发S3的免疫抑制环境，而CTNNB1的高突变可能会抑制S2的免疫浸润。

图3 三种HCC亚型的突变情况的肿瘤间异质性特征 图S4 三种HCC亚型中CNV突变情况的肿瘤间异质性

03 - 整合bulk和单细胞转录组数据得到的新型基因分类器

        为了整合scRNAseq和批量RNA样本，作者首先构建了一个分类器参考集。最初根据基因的平均表达值对GSE149614数据集中的10个scRNAseq样本进行了分组。接下来，计算免疫抑制、活化T细胞（aT）和肝脏的分数，并将其上四分位数定义为阳性样本的截止值。同时，用ESTIMATE软件评估了10个样本的免疫和肿瘤纯度分数。在这10个scRNAseq样本中，发现个别样本HCC02T、HCC03T、HCC04T和HCC05T拥有最高的肝脏和最低的免疫抑制分数（图4A），所以它们被定义为'冷肿瘤'样本。相反，HCC08T、HCC09T和HCC10T具有最高的免疫抑制分数，表明它们是 "免疫抑制肿瘤 "样本（图4B）。同样，HCC01T、HCC06T和HCC07T的活化T细胞得分高，是"热肿瘤 "样本（图4C）。

        接下来，通过选择最特殊的基因来建立一个分类器，以识别上述参考集中的所有阳性样本。使用从三个亚型（S1、S2和S3）的签名中获得的差异表达基因的序列号在参考集中进行分类分析。结果显示，当选择前108个基因时，所有假阳性率（FPR）为0，所有真阳性率（TPR）为1（图4D）。当包括更多的基因时，FPRs和TPRs值对参考集的分类没有改变。

图4 RNA数据与单细胞数据整合，获得基因分类器

        接下来，在TCGA和ICGC队列中使用涵盖这108个基因的面板验证了分类效果。在这两个大队列中，样本也被分成了三个亚型（图S5A,B）。S3在TCGA中仍有最差的预后（图S5A），在ICGC中具有免疫抑制基因的高表达水平、高免疫评分、低肿瘤纯度评分以及最差的预后（CTLA4，TIGIT；图S5B,C,E,F）。此外，从108个基因分类器得到的亚型中的SubMap分析验证了TCGA和ICGC队列之间的一致性（图S5D）。富集的途径在两个队列中也是一致的（图S6）。这些结果进一步证实，108基因分类器可以将单细胞样本映射到三个亚型中，并将HCCs归入不同的免疫状态亚型。

图S5 对分类器的验证

04 - 三个HCC亚型中TME的异质性

        作者研究了免疫抑制亚型的调节机制，分析了三种单细胞亚型（CS1、CS2和CS3）细胞水平的肿瘤免疫微环境的异质性。CS3样本的免疫和免疫抑制得分最高，CS2的肿瘤纯度和肝脏得分最高（图5A）。因此，CS1类似于S1，CS2类似于S2，而CS3类似于S3。经过基因表达归一化、降维、聚类和基于细胞系特异性标记基因的表征，这些细胞被分为16种类型，包括6007个T细胞、1845个B细胞、14 552个上皮细胞、350个NK细胞、1850个内皮细胞和1548个成纤维细胞（图5C）。然后，得到四个T细胞、两个B细胞和三个巨噬细胞亚细胞类型。

        此后，通过GSVA分析，利用获得的免疫学基因组，对免疫细胞进行分类，分析其基因表达谱和功能。在T细胞和NK细胞中，mT是一种记忆性T细胞类型，过量表达IL7R、CCR7和CD69等基因，以及GSVA中幼稚或mT信号的富集。 tT是一种免疫抑制性T细胞类型，HAVCR2、TIGIT、BATF和CTLA4表达水平高（图5D，E，图S7A，B）。具有高水平的MKI67和TIGIT的T细胞被指定为增殖型T细胞（pT）细胞。细胞毒性基因如PRF1和GZMA在aT和NK细胞中高度表达，CD8A在aT中表达，表明它们分别是有效的CD8 T和NK细胞（图S7A）。在GSVA中富含幼稚型或记忆型B细胞信号并过量表达CCR7和CD69的B细胞被指定为记忆型B细胞（mB）。同样地，巨噬细胞也被分为两个亚型。TAM样巨噬细胞（TAM-Mφ）在GSVA中具有高表达水平的M2特征，以及TAM样基因，如APOE、C1QA、SLC40A1和GPNMB（图S7A,C）。然而，MDSC样巨噬细胞（mMφ）同时显示出M1和M2的特征，表现出高表达水平的促炎症基因，如FCN1和VCAN，以及免疫抑制基因IL10（图S7A,C）。此外，CAF表达高水平的αSMA（ACTA2），被指定为肌成纤维细胞（myCAF；图S7C）。组别分布富集结果显示，NK和aT细胞富集在CS1、mT、mB、tT和mMφ，myCAF细胞富集在CS3，上皮细胞（H1、H2、H3、H4）富集在CS2，这进一步验证了基因分类器的分类效果（图5F）。

图5 HCC的肿瘤内异质性 图S7 GSE149614中标记基因的表达情况

05 - 转录因子BATF和MDSC样巨噬细胞可促进免疫抑制性细胞的形成

        作者分析了促进该亚型免疫抑制性T细胞形成的调控机制和细胞相互作用。首先，获得了与CS3特定亚型（mT、tT、myCAF和mMφ）相关的基因模块。将WGCNA与用SCENIC计算的TF结果相结合，挖掘出促进免疫抑制的关键TF调控。在WGCNA分析中，将细胞类型信息作为一种表型，从WGCNA中得到12个基因模块（图6A），其中绿色模块与mT相关，品红色模块与tT相关，黄色模块与myCAF相关，紫色模块与mMφ相关。之后，得到了每个基因模块的枢纽基因。mT枢纽基因在T细胞选择和T细胞差异途径中富集，而tT枢纽基因对T细胞激活和白细胞介素10分泌有负向调节作用（图6B）。随后，在tT枢纽基因和TF调控结果之间进行了超几何检验，并选择了关键的TF；TF，BATF，是被选中的对象之一。值得注意的是，BATF可以调节TIGIT和CTLA4以及共刺激基因ICOS（图6C）。SCENIC中TF激活分数的热图也证实了BATF在TOT中的细胞类型特异性。

        作者进一步证实了BATF、CTLA4和TIGIT在TCGA和ICGC队列中的共同表达。结果显示，这些基因对之间存在明显的共现性（co-occurrence ）（图6D）。在这两个数据集中，BATF的高表达水平与预后不良相关（图6E）。关键的是，来自GSE149197的BATF基因敲除Treg细胞的基因表达数据显示BATF、CTLA4、TIGIT和FOXP3的表达明显降低（图6F），而BATF在健康肝脏单细胞数据集（GSE115469）中几乎没有表达。然后在另外两个单细胞数据集（GSE98638和GSE146409）中验证了BATF的功能。这两个数据集使用了相同的数据处理管道。GSE98638中有九个T细胞亚型，包括两个免疫抑制性T-reg亚型（CD4-CTLA4和CD4-FOXP3），GSE146409中有一个T细胞簇（图S9A,B）。在这两个数据集中，BATF、TIGIT、LAG3和CTLA4也有共表达关系（图S9C,D）。此外，BATF在两个数据集中都能调节免疫抑制基因（图S9E,F）。这表明，TF BATF可以通过上调免疫抑制基因的表达，在形成免疫抑制细胞中发挥关键作用。

图S9 BATF和regulon在其他两个单细胞数据集中的表达情况

        接下来分析了肿瘤微环境细胞在形成肿瘤免疫抑制微环境中的作用。细胞相互作用的结果显示，tT的细胞可以通过趋化因子CXCL12_CXCR4、CCL4_CCR5和CCL3_CCR1与CS3特异性mMφ相互作用。mMφ的特点是免疫抑制基因IL10的过度表达，并且可以通过NECTIN2_TIGIT与tT相互作用。tT可以进一步抑制T细胞的免疫反应，最终促进HCC中产生免疫抑制环境。此外，mMφ经常通过趋化因子如CXCL12_CXCR4和生长因子VEGFA_FLT1与myCAF和内皮型（End）的内皮细胞相互作用（图6G）。同时，End型的内皮细胞也可以通过TIGIT_PVR与tT相互作用，这也可能促进免疫抑制细胞的形成。因此，mMφ可以直接或间接地促进S3类HCC亚型的免疫抑制状态。

图6 单细胞的免疫抑制机制

四、结论

        这项研究利用批量和单细胞转录组数据调查了HCC的肿瘤间和肿瘤内异质性。研究发现了三种不同的亚型，即S1、S2和S3。在这些亚型中，预后最差的S3具有高度的巨噬细胞和免疫抑制性T细胞浸润。结果表明，Treg细胞以及MDSC样巨噬细胞中的TF BATF可以促进免疫抑制细胞的形成并影响HCC患者的预后。这些发现可以促进HCC的临床诊断和治疗。