今天要給大家分享的這篇文章是今年9月底發(fā)表在Aging-US（IF=4.831）上的一篇純生信分析文章（“Glioblastoma cell differentiation trajectory predicts theimmunotherapy response and overall survival of patients”），該文章基于GEO的單細(xì)胞測(cè)序數(shù)據(jù)發(fā)現(xiàn)了具有不同分化特征的膠質(zhì)母細(xì)胞瘤（GBM）細(xì)胞，繼而進(jìn)行差異表達(dá)分析找到分化相關(guān)的基因（GDRG），最后根據(jù)這些基因分別構(gòu)建了分子分型與預(yù)后預(yù)測(cè)模型。

研究背景介紹

GBM是最常見同時(shí)也是死亡率最高的顱內(nèi)原發(fā)腫瘤，5年生存率僅5-6.8%，綜合治療手段，諸如手術(shù)、放療、化療、免疫治療等都未能有效延長(zhǎng)患者生存時(shí)間。腫瘤微環(huán)境中的多種因素會(huì)影響腫瘤細(xì)胞與腫瘤干細(xì)胞（CSC）的分化，從而導(dǎo)致細(xì)胞分化狀態(tài)的異質(zhì)性，而單細(xì)胞轉(zhuǎn)錄組學(xué)分析為我們提供了一個(gè)強(qiáng)大的方法，用于研究整個(gè)腫瘤樣品的基因組本質(zhì)，繼而提供了定義和明確細(xì)胞狀態(tài)及其分化狀態(tài)的機(jī)會(huì)。

圖1

結(jié)果

1.基于單細(xì)胞測(cè)序數(shù)據(jù)確定細(xì)胞亞群

從GSE84456獲取了來(lái)源于4個(gè)GBM標(biāo)本的3589個(gè)細(xì)胞測(cè)序結(jié)果，僅有2343個(gè)來(lái)源于腫瘤核心部位的細(xì)胞被納入進(jìn)一步研究。Seurat R包被應(yīng)用于單細(xì)胞數(shù)據(jù)分析。經(jīng)過(guò)質(zhì)控及標(biāo)準(zhǔn)化數(shù)據(jù)，194個(gè)低質(zhì)量細(xì)胞被排除（圖2A-B）?？偣布{入19752個(gè)基因，方差分析顯示了1500個(gè)高度變異的基因（圖2C）。主成分分析PCA未見GBM細(xì)胞顯著分離趨勢(shì)，我們最終選了20個(gè)主成分行進(jìn)一步分析（圖2D-E）。之后，tSNE分析將GBM細(xì)胞分為13個(gè)亞群，差異表達(dá)分析從所有13個(gè)亞群中鑒定出總共8025個(gè)標(biāo)記基因（圖2F-G）。根據(jù)標(biāo)記基因的表達(dá)量，我們使用singleR和CellMarker對(duì)各個(gè)亞群進(jìn)行注釋（圖3A）。

圖2

2.腫瘤分化軌跡分析并確定GBM分化相關(guān)基因（GDRG）

使用Monocle 2算法進(jìn)行時(shí)序及軌跡分析，發(fā)現(xiàn)腫瘤干細(xì)胞主要位于樹根部（root），而2個(gè)樹枝（branch）分布有不同的GBM細(xì)胞，branch I包含了434個(gè)GBM細(xì)胞，而branch II包含了444個(gè)細(xì)胞（圖3B）。差異表達(dá)分析得到265個(gè)I型GDRG和193個(gè)II型GDRG。基因集富集分析（GSEA）發(fā)現(xiàn)I型GDRG與免疫反應(yīng)調(diào)控相關(guān)通路顯著負(fù)相關(guān)，而II型GDRG與代謝相關(guān)通路顯著正相關(guān)（圖3C-D）。

圖3

3.TCGA中GDRG的分析

為探索能否在TCGA轉(zhuǎn)錄組測(cè)序結(jié)果中找到具有不同分化特征的GBM亞群，相關(guān)性分析表明我們可以通過(guò)傳統(tǒng)轉(zhuǎn)錄組測(cè)序數(shù)據(jù)中GDRG表達(dá)來(lái)鑒定I型和II型GBM細(xì)胞（圖4A-C），因?yàn)檫@些高度相關(guān)的基因可能表明共同的細(xì)胞起源。為了確定來(lái)自不同GBM細(xì)胞亞群的基因是否在功能上相關(guān)，我們利用metagene來(lái)代表相應(yīng)基因的整體表達(dá)模式，即I型和II型的metagene分別由I型和II型GDRG組成。圖4D-F顯示I型和II型的metagene在單細(xì)胞及傳統(tǒng)轉(zhuǎn)錄組數(shù)據(jù)中高度相關(guān)，表明I型和II型GDRG在功能上是高度相互關(guān)聯(lián)的。此外，90.8%的GDRG出現(xiàn)了體細(xì)胞突變。

圖4

4.基于GDRG的分子分型

通過(guò)非監(jiān)督一致性聚類分析，將TCGA患者分為2組，MC1和MC2（圖5A-C）。K-M生存曲線顯示MC1的OS顯著差于MC2（圖5D）。MC1的I型metagene表達(dá)量顯著高于MC2，而MC1的II型metagene表達(dá)量顯著低于MC2。此外，MC1患者的年齡顯著低于MC2患者，而其他臨床指標(biāo)在2組間未見顯著差異。同樣的結(jié)論在CGGA驗(yàn)證集中得到驗(yàn)證。

圖5

5.免疫檢查點(diǎn)分子（ICM）及免疫治療反應(yīng)性分析

PD1、PDL1、PDL2在I型GBM細(xì)胞亞群（單細(xì)胞測(cè)序數(shù)據(jù)）和MC1患者（傳統(tǒng)轉(zhuǎn)錄組測(cè)序）中相對(duì)高表達(dá)，而CTLA4、CD80、CD86在II型GBM細(xì)胞亞群和MC2患者中相對(duì)高表達(dá)（圖6A-C）。根據(jù)TIDE算法，MC2患者對(duì)免疫治療的反應(yīng)性優(yōu)于MC1。此外，SubMap分析發(fā)現(xiàn)MC1對(duì)PD1抑制劑反應(yīng)性更好，而MC2對(duì)CTLA4抑制劑反應(yīng)性更好（圖6D-E）。

圖6

6.構(gòu)建基于GDRG的預(yù)后風(fēng)險(xiǎn)評(píng)分模型

套路化的單因素COX回歸，再接著LASSO、多因素COX回歸，繼而構(gòu)建一個(gè)四基因風(fēng)險(xiǎn)評(píng)分模型。分析評(píng)分計(jì)算公式如下：

根據(jù)風(fēng)險(xiǎn)評(píng)分的中位數(shù)，將所有患者被分為2組，即高危組和低危組，K-M生存曲線提示高危組患者OS顯著差于低危組。ROC曲線及C指數(shù)均提示較好的預(yù)測(cè)價(jià)值。同樣的結(jié)論在CGGA驗(yàn)證集中得到驗(yàn)證。

圖7

7.構(gòu)建預(yù)后相關(guān)的列線圖Nomogram

套路化的單因素+多因素COX回歸，繼而構(gòu)建Nomogram，C指數(shù)、ROC曲線、校正曲線均提示較好的預(yù)后預(yù)測(cè)價(jià)值。同樣的結(jié)論在CGGA驗(yàn)證集中得到驗(yàn)證。

圖8

總結(jié)

該文章首先挖掘GEO數(shù)據(jù)庫(kù)里的單細(xì)胞數(shù)據(jù)，通過(guò)分化軌跡分析找到了與分化相關(guān)的基因，繼而通過(guò)TCGA數(shù)據(jù)庫(kù)以這些基因?yàn)榛A(chǔ)分別構(gòu)建了分子分型以及預(yù)后預(yù)測(cè)模型。該文章思路新穎，可謂生信分析中的一股清流，令編輯和審稿人耳目一新（僅用1個(gè)月即接收），大家不妨按照相似的套路模仿一下，想必不難發(fā)表高分文章。

總結(jié)來(lái)看，這篇文章主要分為兩個(gè)板塊，單細(xì)胞測(cè)序分析+模型構(gòu)建。對(duì)于這兩個(gè)板塊，近期我們就分別有針對(duì)性的課程，歡迎參 加！

[10.17-18 上海班] 單細(xì)胞數(shù)據(jù)挖掘?qū)崙?zhàn)班2020

[10月31-11月1日 在線直播]手把手全流程5分+多組學(xué)生信數(shù)據(jù)挖掘速成