bulk和scRNA揭示肝细胞癌中的细胞异质性和免疫浸润
全文地址：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9168757/
摘要
肝细胞癌（HCC）免疫治疗的有效性因其高度的肿瘤间和肿瘤内异质性以及免疫抑制性肿瘤微环境而受到阻碍。然而，HCC中肿瘤异质性和免疫抑制性微环境之间的相关性尚未得到很好的解决。 在这里，我们剖析HCC中肿瘤间和肿瘤内的异质性，并揭示它们如何对免疫抑制性微环境做出贡献。我们通过非负矩阵分解（NMF）聚类进行共识分子亚型，以分层HCC肿瘤的异质性概况。
我们将来自癌症基因组图谱（TCGA）患者的HCC肿瘤分为三个亚型（S1、S2和S3），
其中S1被描述为具有高T细胞基因表达水平和免疫评分的“热肿瘤”。 S2被描述为肿瘤纯度得分最高的“冷肿瘤”， S3被描述为预后最差的“免疫抑制肿瘤”，免疫抑制基因如细胞毒性T淋巴细胞相关蛋白4、TIGIT和PDCD1的高表达水平。
此外，我们在S3样亚型（CS3）的单细胞数据集中结合加权基因共表达网络分析和单细胞调节网络推理和聚类（SCENIC），并确定了一种转录因子BATF，它可以上调免疫抑制基因。最后，我们确定了一个细胞相互作用网络，其中髓系来源的抑制细胞样巨噬细胞亚型可以促进免疫抑制T细胞的形成。
关键词：肝癌、多组学、单细胞、转录因子、肿瘤异质性
在这项研究中，通过整合批量和单细胞 RNA 测序，我们定义了一个由 108 个基因组成的分类器，并研究了肝细胞癌（HCC）中的细胞异质性以及免疫浸润。在癌症基因组图谱和国际癌症基因组联盟队列中，细胞毒性 T 淋巴细胞相关蛋白 4、TIGIT 和 LAG3 水平高的患者预后最差。此外，单细胞 RNA 测序分析显示，在 CS3 簇中这些免疫抑制基因的表达水平很高，患者生存率低。我们的研究表明，BATF 可能调节肝细胞癌中免疫抑制基因的表达。
介绍
肝细胞癌（HCC）是原发性肝癌最常见的形式，5年生存率为12-18%。肝细胞癌的特点是肿瘤间和肿瘤内的异质性。揭示肝细胞癌异质性背后的分子机制对于靶向治疗的发展至关重要。近年来，对大量肝细胞癌样本的mRNA表达谱进行了全基因组分析。然而，大多数研究都集中在bulk序谱上，阻碍了肿瘤内和肿瘤间异质性的研究。
HCC中免疫治疗的效果通过整体免疫浸润和免疫亚群的丰富共生来调和。肿瘤微环境（TME）中的复杂细胞组成和特征突出了癌症免疫抑制的多种非冗余机制的存在。为了提高免疫治疗疗效，我们应该进一步调查不同患者亚型中免疫细胞的异质性，并确定适合特定免疫治疗的患者。我们还应该研究肿瘤抑制的机制，丰富免疫治疗疗效。近年来，单细胞RNA测序（scRNAseq）已成为揭示肿瘤免疫微环境中细胞异质性的有力工具。例如，据报道，肝内胆管癌（CCA）可以通过TIGIT—PVR的配体受体（L—R）对与调节性T细胞（Treg）相互作用，导致肝内CCA（ICC）中的免疫抑制。据报道，耗尽的CD8+T和Treg细胞在[肝癌中优先富集并可能克隆扩增。然而，尚未评估肿瘤间和肿瘤内异质性之间的相关性和一致性。
在这项研究中，我们对单细胞 RNA 测序（scRNAseq）和多组学数据进行整合分析，揭示了肝细胞癌（HCC）中的肿瘤异质性和免疫抑制机制。这些发现可以促进临床诊断并丰富肝细胞癌的免疫治疗。
材料与方法
一个下载HCC数据的好网址
http://lifeome.net/database/hccdb/download.html
肝癌亚类的识别-NMF
选择平均绝对偏差>1TOP基因的基因进行NMF聚类。随后，使用NMF的r包对归一化表达数据进行无监督NMF聚类方法。选择共遗传相关系数幅度开始下降时的k值作为最佳聚类数。
多组学数据采集与处理
TCGA中“Masked Somatic Mutation”类别的所有肝癌患者的体细胞突变数据使用varscan软件（https://portal.gdc.cancer.gov/）进行处理。使用maftools（版本2.10.0）[13分析和可视化突变。使用单样本基因集富集分析（GSEA，ssGSEA）使用基因集变异分析（“GSVA”）r包（版本1.42.0）[14]评估标志基因的富集分数。标志基因集是从GSEA软件（http：//www. gSeamsigdb.org/gsea/downloads.jsp）的MSigDB数据库中获得的。
差异表达基因分析
“Limma”包被用于执行差异表达基因（DEG）分析。应用经验贝叶斯方法估计NMF聚类方法使用适度t检验鉴定的两个簇之间的差异基因[15]。考虑到S1和S3中的高免疫细胞浸润以及S2中的高肿瘤纯度，我们随后将S1和S3中的差异表达基因与ImmPort（https://immport.niaid.nih.gov/home）的免疫基因重叠，并排除S2中的免疫基因以获得最终的候选基因。 使用Benjamini-Hochberg校正计算多重测试的调整P值。log2倍变化绝对值大于1且FDR<0.05的基因被鉴定为两个簇之间的特征基因集。我们对每个集群进行差异分析，与其他两个集群进行比较，以选择显著上调（log2FC>1；FDR<0.05）或显著下调（log2FC<−1；FDR<0.05）基因。
富集分数
计算富集分数以评估每个簇中的亚型分布。首先，我们计算每个簇中每个亚型的频率。然后，我们将gP除以簇频率（簇中的细胞数除以总细胞数），并获得每个簇中每个亚型的富集分数。
共表达网络分析
使用归一化表达矩阵使用r包构建加权共表达网络（WGCNA）。为了减弱噪声和异常值的影响，对假细胞-pseudocells进行了分析，计算为在每个细胞类型内随机选择的10个细胞的平均值。使用默认参数的block wiseModules函数构建了共表达网络。模块特征烯基因和细胞类型信息之间的相关性使用Pearson测试确定了模块的意义。之后，根据每个基因在中心模块中的模块化连通性和表型性状关系选择中心基因。模块连通性被定义为基因之间Pearson相关性的绝对值（模块成员资格）。临床性状关系被定义为每个基因与细胞类型之间Pearson相关性的绝对值（基因显著性）。我们将模块成员资格设置为>0.7，候选枢纽基因的基因显著性设置为>0.6。
SCENIC分析
我们使用R包SCENIC（https://github.com/aertslab/SCENIC，版本1.1.1–10，RcisTarget1.6.0和AUCell1.8.0）来分析细胞亚型中转录组因子的富集。SCENIC中每个样本的输入矩阵是来自Seurat的原始UMI计数。我们保留在至少0.5%的细胞中检测到的表达总和>3×0.005×细胞数的基因。按照标准的SCENIC程序，我们使用GENIE3方法（对于单个样本）和GRNBoost（对于组合样本）来识别潜在的转录因子（TF）目标。此外，使用AUCell评估每个细胞中每个调节子的活性，AUCell计算曲线下的面积，并整合调节子中所有基因的表达等级。
