1. 引言
肝细胞肝癌(HCC)被认为是世界第六大常见肿瘤和第三大癌症死亡原因，占到了原发性肝癌的80%以上[1]，同时由于肝细胞癌在早期的隐匿性，许多患者在确诊时已经发展到了中、晚期甚至终末期，从而失去根治性手术时间。靶向治疗与免疫治疗近年来在肝细胞癌的治疗领域受到越来越多的关注。然而，由于免疫治疗反应的低应答，很大一部分患者并不能从这种治疗中获益，一些患者在初始应答后很快出现耐药[2]，已有研究证明这与肿瘤的抑制性免疫微环境(TIME)有关[3]。
MDSC细胞(Myeloid-Derived Suppressor Cells)，即髓系免疫抑制细胞，MDSC是具有异质性和免疫抑制特性的未成熟骨髓源性细胞群[4]，是抑制性肿瘤免疫微环境的重要组成部分。MDSC不仅可以抑制NK细胞等免疫杀伤细胞的作用，同时促进调节T细胞(Treg)、肿瘤相关成纤维细胞等其他免疫抑制细胞的产生[5]。已有研究证明MDSC细胞与肝细胞癌的不良预后及免疫治疗的耐药性有关[6]。但如何评价MDSC在肝细胞癌患者中的浸润情况尚未解决。
目前，人类MDSC细胞被标记为CD33+/CD11b+细胞[7]，但MDSC具有高度的异质性，随着技术的发展，MDSC再特定肿瘤中进一步鉴定出了特异的标记与亚群，例如：CD84也被鉴定为乳腺癌患者MDSC的细胞表面标志物[8]，在黑色素瘤或非小细胞肺癌患者中观察到发挥免疫抑制作用的CCR5 + MDSC亚群蓄积[9]。在本文中，我们希望通过单细胞转录组分析联合二代测序数据，鉴定MDSC在肝细胞癌中的特征标记，并开发一种评分工具能够代表肝细胞癌患者MDSC细胞的浸润情况，进一步探究MDSC细胞在肝细胞癌中的免疫特征。
2. 数据与方法
2.1. 数据来源
单细胞测序数据来自GSE189903和GSE156625数据集。该数据集可在GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)中下载，单细胞数据集CNP0000650数据集包括了19名肝细胞癌患者的单细胞券窗转录组数据，该数据集可在中国基因库数据库(https://db.cngb.org/)中公开获取。肝细胞癌的二代测序数据集联合了TCGA数据库中肝细胞癌患者数据和GTEx数据库中的正常肝组织数据。TCGA数据由UCSC Xena (https://portal.gdc.cancer.gov/)下载，GTEx中正常肝组织的数据可在(https://www.genome.gov/)中下载。
2.2. 单细胞数据的清洗质控与降维聚类
按照细胞内检测到的分子总数不大于40,000、每个细胞中检测到的基因数量不大于5000、检测到的线粒体基因不大于20%、检测到的核糖体基因不大于50对数据进行质控。数据质控后利用NormalizeData函数对数据进行标准化，并使用ScaleData函数对单细胞数据进行归一化处理。最后，利用RunHarmony函数对每个样本进行去批次效应。使用RunPCA函数对数据进行线性降维，选择前50个PC对细胞进行聚类，设置分辨率为0.41对细胞进行分群降维。利用常见细胞标记进行人工注释，通过执行Seurat包的FindClusters函数鉴定差异基因。
2.3. 二代测序数据处理
从TCGA数据库下载TCGA-LIHC的TPM表达数据，筛选原发肿瘤数据371例，正常对照数据50例。从GTEx数据库提取正常肝组织表达数据110例。合并TCGA，GTEx数据，利用sva包的Combat函数，设置TCGA及GTEx为批次信息，正常及肿瘤分组为协变量信息，对合并后的数据进行去批次，并使用limma包的normalizeBetweenArrays函数对数据进行归一化。利用秩和检验进行差异分析，并筛选FDR < 0.05，|log2foldChanges| > 1的基因作为差异基因。
2.4. 基因筛选与拟合评分
将随机种子设置为“set.seed(123)”，利用creat包的train函数，使用gbm方法对数据进行Gradient Boosting Machine (GBM)分析，筛选重要性评分大于1的基因作为特征基因，利用randomForest包的randomForest函数对数据进行随机森林分析，设置循环树为500，筛选重要性大于1的基因作为特征基因，利用creat包的rfe函数，使用svmRadial方法进行SVM分析，利用glmnet包的glmnet函数，使用binomial方法进行Lasso分析，使用cv.glmnet进行交叉验证。使用“GSVA”R软件包对TCGA-LIHC数据集进行了单样本GSEA评分。
2.5. 免疫相关分析
利用estimate包对模型进行TME分析，并对高低风险组绘制差异小提琴图。使用CIBERSORT包对数据进行免疫浸润分析，绘制高低风险组之间的免疫细胞浸润差异箱线图。利用cor函数，使用pearson检验分析风险分数与免疫检查点基因之间的相关性，并绘制相关性热图。TIDE数据库是一个用于预测癌症患者对免疫检查点抑制剂(如PD-1/PD-L1抑制剂)反应的算法，将TCGA数据中的371例肝细胞癌患者导入上传网站分析，获得对免疫治疗的反应评分和T细胞排斥评分。利用ggviolin函数绘制高低风险组之间的免疫治疗差异情况。利用oncoPredict包进行药物敏感性分析，使用boxplot函数绘制高低风险组之间药物敏感性差异情况。
3. 结果
3.1. 单细胞分析鉴定MDSC的特征基因
为了鉴定肝细胞癌患者中MDSC细胞的特征基因，我们选择了两个肝细胞癌肿瘤组织单细胞数据集，经过降维聚类后，我们共计得到了6.9w个细胞，这些细胞聚集成了21个簇(见图1(a))，我们选择常用的免疫细胞标记进行人工注释：Mast cell (CD63)，Neutrophil (FCGR3A, ITGAM)，cDendritic cell (FCER1A, CST3)，pDendritic cell (GZMB)，Monocyte (LYZ)，B cell (MS4A1, CD79A)，Plasma cell (MZB1, IGKC, JCHAIN)，NK cell (GNLY, NKG7, KLRD1)，T cell (CD3D, CD3E)，MDSCs (ITGAM, CD33)，得到注释后的TSNE图(见图1(b))。我们发现MDSC细胞主要聚集在第4簇中。通过MDSC与其他免疫细胞的差异分析，我们获得了MDSC的单细胞差异表达基因共计548个。为了进一步验证单细胞差异分析的结果，我们从CNSA数据库中获取了另外一组肝细胞癌患者的单细胞数据集，按照相同的步骤进行数据质控后以相同标记(ITGAM, CD33)注释MDSC细胞，通过差异分析获得差异表达基因共计1274个。我们将两个数据集中MDSC的差异表达基因取交集，单细胞数据分析获得MDSC的差异表达基因共计388个。
Figure 1. (a) TSNE Plot after clustering at resolution choice 0.41; (b) TSNE Plot after cell annotation
图1. (a) 分辨率选择0.41降维聚类后TSNE图；(b) 细胞注释后TSNE图
3.2. 二代测序数据筛选肿瘤特异MDSC基因
我们希望标记肿瘤组织中的MDSC细胞，为了能筛选肿瘤差异表达的基因，我们选择TCGA联合GTEx数据中的肝细胞癌患者与正常肝组织的数据进行差异分析，利用秩和检验的方法筛选出HCC患者差异表达的基因共计2348个。我们将这些基因与MDSC的差异表达基因取交集得到了19个基因。为了进一步筛选MDSC特征性基因，我们利用多种机器学习的方法缩小基因的范围，经过GBM分析、随机森林分析、SVM分析、Lasso分析四种方法(见图2(a))，我们最终获得了11个MDSC的特征基因(见图2(b))。这些基因在肿瘤组织中差异表达同时可以作为MDSC的特征标记物。
Figure 2. (a) Screening for relevant genes through four machine learning methods; (b) Genes after intersection of machine learning results
图2. (a) 四种机器学习方法筛选高相关性基因；(b) 机器学习结果取交集后获得基因
3.3. ssGSEA拟合MDSC的特征评分与预后有关
为了能表征肝细胞癌患者中MDSC细胞的浸润程度，我们希望借用MDSC的特征基因开发一种评价工具。我们使用ssGSEA的方法拟合MDSC的公共标记物和获得的11个特征基因，绘制TCGA数据库中肝细胞癌患者的评分热图(见图3(a))。为了验证MDSC细胞浸润与HCC预后的关系，我们以ssGSEA的中位数为区分值将患者分为MDSC高/低分组，我结合TCGA的临床信息进行生存分析，通过K-M曲线我们发现MDSC高得分组的患者预后生存时间较低得分组明显降低(见图3(b))，具有统计学意义(P < 0.05)，说明MDSC的浸润与HCC患者的不良预后有关。
Figure 3. (a) Heatmap of ssGSEA scoring of TCG-LIHC data; (b) Survival K-M curves for patients in high/low MDSC subgroups
图3. (a) TCG-LIHC数据ssGSEA打分热图；(b) 高/低MDSC分组患者生存K-M曲线
3.4. MDSC的免疫相关分析
为了探究MDSC细胞与其他免疫细胞的关系，我们将患者区分为MDSC高/低得分组，利用cibersort进行免疫浸润分析并绘制箱线图，我们发现在MDSC高得分组中调节T (Treg)细胞、M0髓系细胞和树突状(DC)细胞出现差异高浸润；在MDSC的低得分组中NK细胞、肥大细胞和单核细胞差异高浸润(见图4(a))。这一结果验证了MDSC发挥免疫抑制作用时通过联合Treg等免疫细胞构成免疫抑制微环境，降低NK细胞等免疫杀伤细胞的作用来促进肿瘤的进展。为了评估MDSC不同分组之间免疫微环境的差异。为了评估MDSC得分与常见免疫检查点基因表达的关系，我们绘制了与免疫检查点基因的相关性热图(见图4(b))，我们发现MDSC的评分与免疫检查点基因的表达大多表现为正相关，尤其与HAVCR2、CD86等免疫抑制基因的表达具有较高的相关性。且在MDSC高/低组间，以PD-L1为代表的免疫检查点基因表达存在明显的差异性(见图4(c))。
“***”P < 0.01; “**”P < 0.05; “*”P < 0.1.
Figure 4. (a) Box plot of MDSC score high/low intergroup in immune infiltration analysis; (b) Heatmap of MDSC score correlation with immune checkpoint gene expression; (c) Box plot of differential expression of immune checkpoint genes between high/low MDSC score groups
图4. (a) MDSC得分高/低组间免疫浸润分析箱式图；(b) MDSC得分与免疫检查点基因表达相关性热图；(c) MDSC得分高/低组间免疫检查点基因差异表达箱式图
为了探究MDSC细胞浸润对免疫检查点抑制剂(ICB)治疗的关系，我们借助TIDE数据库进行了免疫治疗预测分析，通过比较MDSC高/低表达组间对抗PD-1和抗CTLA-4的应答效率，我们发现在抗PD-1应答组(见图5(a))、抗CTLA-4应答组(见图5(b))和两种治疗均应答组(见图5(c))中，MDSC低浸润组均表现出了差异性的高得分(P < 0.05)。这一结果表明MDSC的浸润与肝细胞癌患者免疫治疗的抵抗有关。为了探究MDSC浸润与常见靶向和化疗药物的敏感性是否存在相关，我们借助oncoPredict包对MDSC的高/低得分组进行了药敏分析。我们发现在肝细胞癌常用TKI靶向药物索拉菲尼的分析中，MDSC高得分组的IC50值反而更小且与低得分组差异具有统计学意义(P < 0.05) (见图5(d))。这提示我们MDSC高浸润的患者在靶向治疗中可能会获得更好的结果。总之，以上结果表明MDSC的浸润在肝细胞癌的免疫治疗中表现出了不良影响。
Figure 5. (a) Violin plot of treatment responsiveness between MDSC high/low groups in anti-PD-1 positive conditions; (b) Violin plot of treatment responsiveness between MDSC high/low groups in anti-CTLA-4 positive conditions; (c) Violin plot of treatment responsiveness between MDSC high/low groups in anti-CTLA-4 and PD-1 positive conditions; (d) Box plot of differences in sorafenib drug sensitivity between high/low MDSC scores
图5. (a) 抗PD-1应答组MDSC组间治疗反应性；(b) 抗CTLA-4应答组MDSC组间治疗反应性；(c) 抗PD-1/CTLA-4应答组MDSC组间治疗反应性；(d) MDSC得分高/低间索拉菲尼药敏差异性箱式图
4. 讨论
骨髓生成是健康人体中的常见过程。从造血干细胞到骨髓前体细胞再到成熟免疫细胞的转变受到严格调控[10]。然而，在患有肿瘤或慢性炎症的患者中，连续刺激通常导致未成熟骨髓细胞的分化缺陷[11]，这些细胞倾向于分泌高水平的精氨酸酶1 (ARG 1)、活性氧(ROS)、髓过氧化物酶(MPO)、一氧化氮(NO)和其他产物以发挥免疫抑制作用[10]-[12]，因此被命名为髓源性抑制细胞。MDSC是高度异质的群体，基于形态学特征MDSC细胞常被分类为单核细MDSC (M-MDSC)和多形核MDSC (PMN-MDSC) [13]。但在具体肿瘤组织中的表型和特征有待进一步探究，我们利用单细胞转录组分析联合bulk测序数据，鉴定了肝细胞癌中MDSC的特征基因，并利用ssGSEA方法开发了一种能够评价MDSC细胞浸润的评分工具。已有研究证明MDSC的浸润与多种癌症的不良预后有关，我们在肝细胞癌中也验证了这一结果。
近年，ICB等免疫治疗药物在肝细胞癌中的应用愈加广泛，但并非所有患者对免疫治疗都表现出积极的反应，且许多接受免疫治疗的患者很快出现耐药性[14]。已有研究证明，MDSC在肿瘤免疫抑制微环境中起主导作用，可能是CPI治疗耐药的主要原因。在本文的研究中，我们通过评分工具有效地判断患者MDSC的浸润程度，并通过预测免疫治疗的反应性验证了MDSC在肝细胞癌患者抗PD-1/CTLA-4治疗中具有抵抗作用。这为解决免疫治疗的耐药提供了新的思路，ICB联合其他免疫制剂的用药方法已成为近年研究的热点[15]。同时，靶向MDSC的治疗药物近年来亦有报道，例如全反式维甲酸(ATRA)通过诱导MDSC的分化减少MDSC的浸润[16]，NLRP通路抑制剂可以通过阻止MDSC的招募减少MDSC的浸润[17]。未来，抗MDSC药物联合ICB治疗的方案可能是免疫治疗的新希望。
尽管我们探究了肝细胞癌患者中MDSC的特征标记，但由于条件限制，我们的研究缺少实验进一步的验证，同时对MDSC在肝细胞癌中的作用机制仍缺少探讨。我们希望在未来的研究中能够进一步探究MDSC在肝细胞癌患者中的作用。
5. 结论
总而言之，我们利用生信分析工具，找到了肝细胞癌中MDSC细胞的特征标记基因，并基于此开发了一种评价MDSC细胞浸润的打分工具。同时，验证了MDSC与免疫检查点基因的表达正相关而对免疫检查点抑制剂的治疗具有抵抗作用。
NOTES
*通讯作者。