今天給同學們分享一篇關于腫瘤+免疫治療構建預后風險模型的單細胞生信文章 “ Integration of scRNA-Seq and Bulk RNA-Seq to Analyse the Heterogeneity of Ovarian Cancer Immune Cells and Establish a Molecular Risk Model ”，這篇文章于2021年9月發(fā)表在 Front Oncol 雜志上，影響因子=6.244。在這項研究中，作者構建了一系列組織特異性簇來預測來自GEO數據庫中兩個OC scRNA-seq數據集的免疫細胞組成。進行免疫細胞簇分析和單變量Cox回歸分析和逐步回歸以建立特征，并進行qPCR和免疫組織化學分析以評估特征基因在OC中的表達。最后，基于OC免疫細胞的特異性建立了雙基因特征預后分層系統(tǒng)（CXCL13和IL26），以識別潛在的免疫治療靶點并準確評估預后風險。

1.scRNA-Seq數據的整合和聚類

兩個scRNA-seq數據集用于表明GEO數據庫中的OC特異性。為了整合兩個單細胞轉錄組數據集，使用Seurat包SCTransform函數進行預處理并降低批處理效應。均勻流形近似和UMAP用于非線性降維（圖 2A）。FindCluster函數用于對細胞進行聚類，獲得20個簇（圖 2B）。

T細胞或NK細胞、B細胞和骨髓細胞根據免疫細胞標志物進行聚類（PTPRC是immune cell marker；EPCAM是epithelial cell marker；COL1A2是fibroblast marker；IL7R是naive T cell marker；CD8A和NKG7是CD8+ the T cell and NK cell markers） （圖 2C）。 

圖2 OC scRNA-seq的降維

2.免疫細胞生信分析

T細胞或NK細胞、B細胞和骨髓細胞的聚類分析基于免疫細胞標記物（圖 3A）。首先，作者對T細胞進行分類和識別，然后根據每個細胞樣本的GSVA富集評分進行聚類分析。根據T細胞功能狀態(tài)，如調節(jié)性、共刺激性、初始性、細胞毒性和窮竭性，鑒定了T細胞或NK細胞、B細胞和骨髓細胞的基因表達特征（圖 3B）。 

其次，分析了B細胞的功能狀態(tài)，如抗凋亡、幼稚記憶、細胞因子、增殖和生發(fā)中心基因表達特征（圖3C）。對于腫瘤浸潤的骨髓細胞，分析了M2樣和M1樣骨髓細胞的活性，結果表明M1和M2相關基因在GSE154600的P3和P4患者中顯著上調（圖 3D）。 

圖3 免疫細胞的GSVA富集分析

3.CIBERSORT免疫浸潤

基于單細胞測序數據和免疫細胞類型分析的結果，作者使用批量數據進行臨床分析和預后模型構建。由于批量RNA-seq數據具有樣本多、臨床信息多的優(yōu)勢，為了進一步分析OC浸潤的免疫細胞的臨床意義。CIBERSORT可以根據RNA-seq計數數據預測22個免疫細胞的比例，并用于計算378名TCGA-OC患者的大量 RNA-seq數據中M1樣TAM（腫瘤相關巨噬細胞）的豐度 (圖 4A)。生存分析結果顯示，M1樣TAMS豐度高的患者生存率更高（圖 4B、C）。具有M2樣TAMS比例的患者之間沒有顯著的生存差異，因此，作者對M1樣TAMS進行了深入的生信分析。

圖4 免疫細胞比例

4.WGCNA分析和免疫治療預測

為了進一步探索M1樣TAM在OC中的潛在作用，作者基于TCGA數據進行了WGCNA分析。中值絕對偏差(MAD)≤0.01的基因被過濾掉，留下35,165個基因。軟閾值設置為5（圖 5A、B），構建了一個無標度共表達網絡，以識別與M1樣 TAM相關的基因特征。共生成了7個模塊（圖 5C，D），其中棕色模塊（3213個基因）與M1樣TAM評分的相關性最高（r=0.42，P=2e-17，圖 5E）。如圖5F所示，基因表示為點，橫坐標模塊表示基因與模塊特征基因之間的相關性，縱坐標表示基因表達與OS之間的相關性。結果表明棕色模塊代表OS相關基因，最終從模塊中獲得45個中心基因（MM>0.5和GS>0.1）。

IMvigor210CoreBiology軟件包包含348名PD-L1免疫治療腫瘤患者的RNA-seq數據，分為inflamed、immune excluded和immune desert三種表型。作者使用e1071包構建了支持向量機來預測三種表型（圖5G），表明用于區(qū)分inflamed和immune desert時預測效果更好（圖5H）。這些結果表明，45個M1樣骨髓細胞相關基因是免疫浸潤的潛在預測因子。

圖5 Hub 基因篩選和免疫治療預測

5.基于M1相關基因的分子分型

首先，從TCGA數據庫中提取45個M1相關Hub基因的表達，并使用NMF包基于非負矩陣分解將TCGA樣本劃分為不同的亞組。Cluster = 2作為最佳參數（圖 6A），訓練集分為兩個組。然后，建立一致性矩陣（圖6B），一致性矩陣的值為[0,1]，1表示多次聚類且兩個數據點都屬于同一類，0表示不多次聚類在同一個類。熱圖顯示45個M1相關基因的表達（圖 6C），簇2的預后優(yōu)于簇1（圖 6D）。小提琴圖顯示簇2中M1的比例高于簇1（P=2.865e-07，Wilcox檢驗）（圖 6E）。一般來說，簇1患者的預后較差。通過limma包鑒定了簇2和簇1中差異表達的基因（|logFC|>2和adj.P.val<0.05），獲得658個差異基因，其中39個基因下調和619個基因在簇2中上調。Clusterprofiler包用于在簇2中執(zhí)行KEGG富集分析（圖 6F，G）。

圖6 M1相關分子分型

6.構建遺傳預后模型

為了便于后續(xù)驗證，從658個差異基因中選擇了101個蛋白質編碼基因進行后續(xù)分析。TCGA按照1:1的比例隨機分為訓練集和測試集，并構建預后模型，每個數據集中有189個樣本。Cox分析用于識別訓練集中的四個生存相關基因（CXCL13、PLA2G2D、IL26、CARD17）。然后，使用lasso回歸分析解決過程中的多重共線性問題，減少預后模型中的基因數量。作者使用glmnet包進行l(wèi)asso Cox回歸分析，每個自變量的變化軌跡如圖所示。隨著lambda逐漸增加，自變量系數的數量趨于0時逐漸增加（圖 7A）。接下來，作者使用10倍交叉檢驗來構建每個lambda下的模型和置信區(qū)間，如圖7B所示。當lambda = 0.52并且選擇兩個基因（CXCL13，IL26）來構建預后模型時，該模型是最佳的。 

作者計算了TCGA訓練集的風險評分并確定了風險評分分布，表明CXCL13和IL26 基因表達較低的患者的風險評分和死亡率越高（圖 7C-E）。中位風險評分標準化為0，樣本按照中位標準化分為高風險或低風險，高危組的預后較差（圖 7F）。 

圖7 構建遺傳預后模型

7.預后預后模型的驗證

為了確定預后模型的穩(wěn)健性，作者使用TCGA測試（圖 8A-C）和所有TCGA數據集（圖 8E-G）來計算風險評分，表明風險評分高的樣本明顯小于風險評分低的樣本。CXCL13和IL26的低表達被確定為風險因素。最后，圖 8D、H 中所示的 KM 曲線的結果揭示了低風險組和高風險組之間的顯著差異（p < 0.05）。 

圖 8預后模型的驗證

8.OC組織中特征基因的表達

此外，為了驗證雙基因特征的準確性，作者通過qPCR（圖 9A、B）和IHC（圖 9C、D）分析檢查了OC患者臨床樣本中特征基因（CXCL13 和 IL26）的表達，表明OC組織中CXCL13和IL26的表達較低。

圖9 qPCR和IHC結果顯示