基于單細胞測序和轉錄組分析的胰腺癌壞死相關預后模型

胰腺癌作為一種死亡率高、療效差的腫瘤類型，其發(fā)病機制尚不清楚。目前有作者通過TCGA數(shù)據(jù)庫、ICGC數(shù)據(jù)庫和GEO數(shù)據(jù)庫（GSE85916），利用單細胞測序分析和轉錄組分析構建了胰腺癌壞死性凋亡相關的標記，構建胰腺癌預后模型，為胰腺癌的診斷、治療及預后評估提供參考。該研究發(fā)表在《Frontiers in Immunology》，IF：8.786。

實驗方法：

將TCGA隊列設置為訓練隊列，將ICGC和GSE85916隊列設置為驗證隊列。胰腺癌單細胞測序數(shù)據(jù)來自GEO數(shù)據(jù)庫中的GSE154778。通過加權共表達網絡分析和單細胞測序分析，確定了與壞死最相關的基因。對這些基因進行COX回歸和Lasso回歸，建立預后模型。通過計算風險評分，將胰腺癌患者分為NCPTS_high組和NCPTS_low組，進行組間生存分析、免疫浸潤分析和突變分析。通過基因敲除、CCK-8實驗、克隆形成實驗、transwell實驗和創(chuàng)面愈合實驗等細胞實驗，探討關鍵基因EPS8在胰腺癌中的作用。對臨床樣本進行PCR檢測，驗證EPS8的表達。

技術路線：

主要研究結果：

1. 加權基因共表達網絡分析

該研究的流程圖如圖1所示。如圖2A, B所示，采用ssGSEA對每個樣本進行壞死性凋亡（necroptosis）評分。根據(jù)壞死評分中位數(shù)將患者分為高壞死組和低壞死組，較高的壞死評分胰腺癌患者死亡率很高。此外，高壞死組患者預后較差（P <0.05），提示壞死性凋亡是胰腺癌的危險因素。使用WGCNA進一步尋找與壞死性凋亡共變的基因集。如圖2C所示，當軟域值為7，R^2>0.8時，數(shù)據(jù)更符合冪律分布，平均連通性趨于穩(wěn)定，適合后續(xù)分析。如圖2D所示，共得到27個非灰色模塊。然后，將模塊的相似度域值設置為0.4，對低于該值的模塊進行合并，最終得到16個非灰度模塊。作者發(fā)現(xiàn)，MEgreenyellow模塊與壞死性凋亡密切相關，包含3352個基因（COR = 0.8, P <0.001）（圖2E, F）。

圖1實驗流程圖

圖2單樣本基因集富集分析（ssGSEA）和加權基因共表達網絡分析（WGCNA）

2. 單細胞測序分析

如圖3A所示，研究共納入14個樣本，各樣本間細胞分布較為均勻，說明樣本間不存在明顯的批量效應，可用于后續(xù)分析。隨后，所有細胞聚類為14個簇（圖3B），分為8種細胞類型（圖3C），如上皮細胞、巨噬細胞和T細胞。根據(jù)每個細胞中壞死基因的特點，確定每個細胞中壞死基因的百分比。根據(jù)中位數(shù)將細胞分為高、低兩組，即高NCPTS組和低NCPTS組（圖3D）。經P <0.05差異分析，兩組間有2518個差異表達基因。如圖3E所示，變化最顯著的前10個基因分別是RPL8、AIF1、HLA-DQB1、LST1、HLA-DPA1、RPL7、MS4A7、SPINK1、RPL13A、RPLP0。同時，基因本體富集分析表明，這些基因主要與mRNA代謝和蛋白質定位有關（圖3F）。最后，將差異表達基因與綠綠黃模塊基因交叉得到805個與壞死性凋亡密切相關的基因。

圖3單細胞測序分析

3. TCGA隊列中預后模型的構建與評價

為了進一步鑒定與預后相關的壞死基因，對TCGA和PA-AU組805個壞死基因進行單因素Cox分析，鑒定出48個壞死基因，（p<0.05）。TCGA隊列采用forest-map顯示48個基因的單變量Cox分析結果（圖4A），然后進行LASSO回歸。如圖4B, C所示，最優(yōu)Lambda為0.0678，最終得到由8個基因構建的模型。根據(jù)中位數(shù)將患者分為NCPTS_high組和NCPTS_low組。如圖4D所示，NCPTS_high組患者預后不良（P<0.001）。ROC曲線顯示，模型值NCPTS預測患者預后的AUC維持在0.75左右，明顯優(yōu)于性別、年齡、分期等臨床特征（圖4E）。同樣，決策曲線分析發(fā)現(xiàn)，與臨床特征相比，患者從基于NCPTS的臨床干預中獲益最多（圖4F）。多變量Cox分析顯示，NCPTS是胰腺癌患者的獨立預后因素（圖4G）。此外，該模型可以更好地將患者區(qū)分為NCPTS_high和NCPTS_low組（圖4H）。

圖4 TCGA隊列預后模型的建立

4. 預后模型的外部驗證和評估

在外部驗證隊列PA-AU（圖5A-E）、PA-CA（圖5F-J）和GSE85916（圖5K-O）中也觀察到了類似的結果。在PA-AU隊列（圖5A-E）和PA-CA隊列（圖5F-J）中，NCPTS_high組患者的預后較NCPTS_low組患者差。此外，預后ROC分析顯示，曲線下AUC波動在0.85 ~ 0.7之間，對患者的預后評估有較好的效果，優(yōu)于性別、年齡等臨床特征（圖5B, G）。決策曲線分析也顯示，以NCPTS為干預依據(jù)的患者獲益最大（圖5C, H）。多因素Cox分析顯示，NCPTS是影響患者預后的獨立因素（圖5D, I），并且該模型能夠很好地將患者區(qū)分為NCPTS_high和NCPTS_low組（圖5E, J）。同樣，在GSE85916隊列中，發(fā)現(xiàn)NCPTS_high患者的預后較差，通過預后ROC分析發(fā)現(xiàn)曲線下AUC較為穩(wěn)定，波動為0.65（圖5K, L）。通過決策曲線分析也表明基于NCPTS的臨床干預具有合理的效益效果（圖5M）。單因素Cox分析顯示，NCPTS是一個預后因素（圖5N）。此外，該模型仍然可以將患者區(qū)分為NCPTS_high和NCPTS_low組（圖5O）。

圖5預測模型的外部驗證和評估

5. 免疫浸潤分析和識別更適合免疫治療的患者

免疫微環(huán)境在患者預后中起著重要作用。分析NCPTS_low組和NCPTS_high組的免疫浸潤情況。如圖6A-D所示，與NCPTS_high組相比，NCPTS_low組的估計評分、免疫評分和間質評分較高（P<0.001），但腫瘤純度較低。與上述結果類似， NCPTS_low組的T細胞、B細胞等免疫細胞浸潤水平更高（圖6E）。免疫原性細胞死亡基因（ICDs）、人類白細胞抗原（HLAs）和免疫檢查點（ICPs）也與腫瘤的發(fā)生發(fā)展密切相關（圖6F）。在NCPTS_low組中，ICDs表達整體呈下降趨勢，HGF、TLR4、P2RX7、FPR1表達上調 （圖6G）。NCPTS_low組HLAs整體表達上調，而HLA-G表達下調（圖6H）。NCPTS_low組ICPs表達普遍上調，如PDCD1、CD27、CTLA4等。然后，根據(jù)模型中的基因對患者進行無監(jiān)督聚類分析，發(fā)現(xiàn)患者可以聚為2類（圖6I-J）。與Cluster1相比，Cluster2的預后較差（P<0.05）。隨后，探究NCPTS分組與聚類以及免疫亞型之間的關系（圖6K）：Cluster1主要與NCPTS_low組相關，免疫亞型C3和C6主要分布在NCPTS_low組。Cluster2主要與NCPTS_high組相關，免疫亞型C1和C2主要分布于NCPTS_high組。隨后，進行免疫亞型之間的生存分析（圖6L）：發(fā)現(xiàn)免疫亞型C3組預后最好，而C3主要分布在NCPTS_low組，這與前期結果一致。綜上，NCPTS_low組有較高水平的免疫細胞浸潤，ICPs和HLAs的表達也較高，且免疫亞型C2主要分布在該組。作者假設NCPTS_low組可能免疫“熱”腫瘤，而NCPTS_HIGH組可能免疫“冷”腫瘤。然后觀察NCPTS組之間的免疫治療反應情況（圖6M, N），發(fā)現(xiàn)NCPTS_low組的TIDE較低，說明免疫逃逸腫瘤的可能性較低，患者獲得免疫治療獲益的可能性較大（圖6M）。此外，TIDE與NCPTS呈強正相關關系（圖6N），即隨著NCPTS的增加，TIDE值也隨之增加，腫瘤免疫逃逸的概率也隨之增加，免疫治療的獲益更少。

圖6 TCGA隊列中更適合免疫治療的患者免疫浸潤分析及鑒別

6. 突變景觀分析

基因突變在腫瘤發(fā)生發(fā)展和患者預后中也起著重要作用。在NCPTS_high組和NCPTS_low組中，突變頻率最高的5個基因分別是KRAS、TP53、SMAD4、CDKN2A和TTN（圖7A，B）。NCPTS_high組TMB較高，且TMB與NCPTS呈正相關（P<0.001）（圖7C，D）。TMB常與患者預后不良相關，這可能是NCPTS_high組預后不良的原因。隨后，對兩組通路進行富集分析，發(fā)現(xiàn)NCPTS_high組激活了與細胞周期相關的G2M_CHECKPOINT和E2F_TARGETS通路（圖7E，F(xiàn)）。

圖7突變景觀分析

7. 模型的臨床意義

對不同隊列NCPTS的預后影響進行了薈萃分析。在圖8A中，發(fā)現(xiàn)合并的HR值表明NCPTS仍然是患者預后不良的一個促成因素。為了進一步促進臨床胰腺癌患者風險組的分類和管理，對NCPTS和臨床特征進行了決策樹分析。根據(jù)NCPTS水平、性別、N期可將患者分為4個危險組：RIS1、RIS2、RIS3、RIS4（圖8B）。RIS4患者均分布在NCPTS_high組，而RIS1、RIS2、RIS3患者分布在NCPTS_low組（圖8C）。生存分析發(fā)現(xiàn)，RIS4組預后不良（P<0.001，圖8D）。最后構建nomogram圖， TCGA-2J-AABK患者1年、3年和5年的死亡率分別為18.5%、62.5%和75.2%（圖8E）。預后ROC曲線分析和決策曲線分析進一步評價Nomogram應用價值（圖8F）。Nomogram預測患者預后的AUC值保持在0.8左右，明顯高于其他臨床特征?；贜omogram臨床及時治療患者的獲益率高于其他臨床特征（圖8G）。

圖8模型的臨床意義

8. 體外驗證關鍵基因EPS8在胰腺癌細胞株中的作用

在CAPAN-1和PANC-1細胞系中，EPS8顯著減少（圖9A）。在CAPAN-1和PANC-1細胞系中，EPS8敲除后胰腺癌細胞的活性顯著降低（圖9B）。隨后，克隆形成試驗顯示，EPS8敲除后，CAPAN-1細胞系（圖9C）和PANC-1細胞系（圖9D）產生菌落的能力大大降低。在transwell實驗中，EPS8的敲除顯著降低了CAPAN-1細胞系和PANC-1細胞系中胰腺癌細胞的遷移和侵襲能力（圖9E, F）。創(chuàng)面愈合實驗中，si-EPS8組EPS8在CAPAN-1細胞和PANC-1細胞中的遷移能力較弱（圖9G, H）。

圖9體外驗證了關鍵基因EPS8在胰腺癌細胞系中的作用

9. 通過對臨床標本進行PCR檢測，驗證EPS8在胰腺癌組織中的表達

為了進一步驗證EPS8在臨床標本中的表達，進行PCR檢測。結果發(fā)現(xiàn)，與癌旁組織相比，胰腺癌組織中EPS8表達明顯上調（圖10）。

圖10 EPS8表達的臨床驗證

結論：

通過結合單細胞測序和轉錄組分析的結果，作者構建了胰腺癌的壞死性凋亡相關預后特征，該的預測特征可以有效評估胰腺癌患者的預后，在一定程度上為胰腺癌的治療提供參考。