近些年,單細胞技術(shù)的興起為科研人員的精準研究提供了堅實的基礎(chǔ),單細胞轉(zhuǎn)錄組數(shù)據(jù)與公共大隊列數(shù)據(jù)的多組學整合已然成為癌癥治療研究的新熱點。與此同時,一個關(guān)鍵問題誕生了---如何更好地將兩種類型數(shù)據(jù)進行整合從而挖掘重要癌癥的發(fā)生發(fā)展機制呢?接下來,小編將給大家介紹兩種不同的單細胞轉(zhuǎn)錄組數(shù)據(jù)與大隊列bulk轉(zhuǎn)錄組數(shù)據(jù)整合的算法:BayesPrism和Scissor
單細胞轉(zhuǎn)錄組測序聯(lián)合bulk轉(zhuǎn)錄組測序(整合篇)
一.BayesPrism算法介紹
下面我就先來具體介紹BayesPrism這項整合方法。為什么選擇這個方法呢?因為我覺得這個方法很有潛力,作者通過比較當前主流的方法如CIBERSORT,證明BayesPrism方法從bulk RNA-seq數(shù)據(jù)集推斷單細胞細胞類型以及對整體基因表達變化的解析上具有比較明顯的優(yōu)勢。簡單來說,BayesPrism算法以scRNA-seq作為先驗信息,從bulk RNA-seq數(shù)據(jù)集中預(yù)測細胞類型的組成和基因表達。該文章主要對膠質(zhì)母細胞瘤(GBM)、頭頸部鱗狀細胞癌(HNSCC)和皮膚黑色素瘤(SKCM)進行了整合分析,以將細胞類型組成與不同腫瘤類型的臨床結(jié)果相關(guān)聯(lián),并探索惡性腫瘤細胞和非惡性細胞狀態(tài)的空間異質(zhì)性。
本篇文章發(fā)表在期刊: Nature Cancer 該期刊在最近一年的影響因子: 23.013。
1. BayesPrism算法研發(fā)背景
大量研究表明,腫瘤微環(huán)境(TME)中惡性細胞和不同類型的非惡性細胞之間的相互作用可以促進血管生成、癌癥轉(zhuǎn)移和免疫抑制等。非惡性細胞在不同的患者和腫瘤類型之間存在顯著差異,某些非惡性細胞群常被用作臨床生物標志物和治療靶點。單細胞轉(zhuǎn)錄組測序(scRNA-seq)技術(shù)的興起,使得在TME內(nèi)對單個細胞的轉(zhuǎn)錄組進行直接的全基因組測量和表征其異質(zhì)性成為可能。然而,scRNA-seq的成本和對高質(zhì)量樣本的要求限制了可檢測的樣本的數(shù)量。此外,scRNA-seq在細胞捕獲中容易受到技術(shù)偏差的影響,這妨礙了細胞類型組成的復(fù)原。為了盡量規(guī)避以上問題,研究人員擬從bulk RNA-seq數(shù)據(jù)中推斷單細胞類型組成與基因表達。然而,現(xiàn)有的反卷積方法沒有完全支持異質(zhì)腫瘤細胞群體中基因表達的預(yù)測。因此,現(xiàn)有的方法未能解決這些關(guān)鍵問題:在TME中,惡性細胞如何影響非惡性細胞的組成?哪些基因與這些相互作用相關(guān)?為了回答這些問題,研究者創(chuàng)建模型,它可以準確地表示每個bulk RNA-seq樣本中的細胞類型比例和細胞類型特異性基因表達譜。因此,該研究提出了BayesPrism,一個貝葉斯模型,使用scRNA-seq參考表達譜作為先驗信息,從bulk RNA-seq數(shù)據(jù)中推斷細胞類型組成和基因表達。
2. BayesPrism算法推斷細胞類型組成和細胞類型特異性基因表達
BayesPrism算法具體步驟:BayesPrism使用細胞狀態(tài)(cell state)、細胞類型、基因列表、bulk RNA-seq數(shù)據(jù)集樣本個數(shù)、bulk RNA-seq 表達矩陣以及單細胞表達矩陣作為輸入文件,先模擬出每個細胞狀態(tài)的基因表達矩陣和每個細胞狀態(tài)的比例,最后通過對每種細胞狀態(tài)的求和,估算出細胞類型的比例和各細胞類型特異的基因表達水平。
為了驗證BayesPrism估計細胞類型比例的效能,作者生成模擬數(shù)據(jù),并發(fā)現(xiàn)BayesPrism在估計惡性腫瘤細胞比例方面比CIBERSORTx的效果更好(圖1 c, d)。同時,以PBMC scRNA-seq數(shù)據(jù)為參考數(shù)據(jù)集,相對于經(jīng)典的方法如CIBERSORTx S mode和MuSiC, BayesPrism 對B細胞、髓系細胞、CD4+ T 細胞、CD8+ T細胞以及NK細胞進行了更準確的細胞類型估計(圖1e,f)??傊?,這些數(shù)據(jù)分析表明BayesPrism完善了現(xiàn)在的反卷積算法的性能。該研究在惡性細胞中使用BayesPrism算法,進行跨患者異質(zhì)性的基因表達的預(yù)測。首先使用了來自8個GBM的scRNA-seq作為參考數(shù)據(jù),估計了GBM bulk RNA-seq數(shù)據(jù)中的細胞類型和基因表達,并發(fā)現(xiàn)bulk RNA-seq數(shù)據(jù)中惡性細胞的基因表達的估計與已知的真實表達高度相似(圖1g)。并且研究發(fā)現(xiàn),當腫瘤純度大于50%,BayesPrism基因表達估計值與已知基因表達真實值的相關(guān)性要大于0.95(圖1h)。這表明BayesPrism評估的基因表達可以準確地從bulk樣本中復(fù)原惡性細胞中的基因表達,并且,使用BayesPrism對基因表達的估計要比使用CIBERSORTx或無反卷積更準確(圖1h)。
3. 浸潤免疫細胞類型和狀態(tài)對生存的影響
該工作分析了TCGA中GBM、HNSCC和SKCM三種腫瘤的1142個bulk樣本中細胞類型的比例。利用GBM、HNSCC和SKCM三種腫瘤的單細胞參考數(shù)據(jù)集,估計了6種GBM細胞類型,10種HNSCC細胞類型,8種SKCM細胞類型(圖2a)。接下來,使用兩個Cox比例風險模型檢查了細胞類型比例和生存率之間的關(guān)聯(lián),與以往報道一致,在SKCM中,發(fā)現(xiàn)CD8+ T細胞和巨噬細胞比例與生存有更強的相關(guān)性(圖2b,c)。并且,根據(jù)M1和M2兩個巨噬細胞亞群的標記基因, 發(fā)現(xiàn)而來自SKCM的M2型巨噬細胞評分最低,M1型巨噬細胞評分與來自HNSCC的巨噬細胞無顯著差異(圖2d)。此外,在SKCM中,巨噬細胞高M1極化低M2極化狀態(tài)與生存率有極強的相關(guān)性(圖2e)。這些發(fā)現(xiàn)表明了巨噬細胞比例以及巨噬細胞狀態(tài)對不同惡性腫瘤具有不同臨床結(jié)果的影響。
4. BayesPrism識別惡性細胞固有基因程序
作者在BayesPrism中開發(fā)了一個模塊,在剔除非惡性細胞類型的基因表達后,用來推斷惡性細胞固有基因程序,該基因程序可以更好地解釋bulk RNA-seq數(shù)據(jù)中的表達異質(zhì)性(圖3a)。BayesPrism模擬得到的基因程序與對惡性細胞進行因子分解得到的基因程序基本一致(圖3b)。BayesPrism識別到的每個基因程序的權(quán)重與每個腫瘤中分配給每個主要亞型的細胞比例相關(guān)(圖3c,d)。BayesPrism發(fā)現(xiàn)GBM中的幾個程序與以往的研究相似,包括program 3(classical和AC-like亞型),program 4(mesenchymal)和progam 5(proneural, OPC和NPC-like)(圖3e)。在HNSCC中,program 1富集了以往單細胞研究確定的Partial EMT亞型(圖3f)。在SKCM中,發(fā)現(xiàn)了多個生存相關(guān)的基因程序,以及一個T細胞排斥程序(圖3g-j)。
5. GBM基因程序和細胞類型的空間異質(zhì)性
作者將122個樣本的bulk RNA-seq GBM數(shù)據(jù)集解析為成5個空間結(jié)構(gòu): LE、IT、CT、MVP和PAN(圖4a)。值得注意的是,已知這些不同結(jié)構(gòu)的TME在血液供應(yīng)、氧水平和免疫應(yīng)激這幾個方面也呈現(xiàn)很大的不同,所有這些改變都可能影響細胞類型組成和惡性細胞狀態(tài)。利用GBM scRNA-seq數(shù)據(jù)進行了反卷積,檢查在解剖結(jié)構(gòu)中,哪些細胞類型和基因程序富集(圖4b,c)。根據(jù)相應(yīng)的解剖結(jié)構(gòu),正如預(yù)期的那樣,MVP區(qū)域在內(nèi)皮細胞和周細胞中高度富集,而LE和IT區(qū)域在少突膠質(zhì)細胞和神經(jīng)元中富集。值得注意的是,PAN區(qū)域在巨噬細胞和T細胞中富集。將解剖結(jié)構(gòu)與基因程序聯(lián)系起來,發(fā)現(xiàn)LE和IT區(qū)域在program 1和2中富集,CT區(qū)域在program 3中富集,PAN區(qū)域在program 4和program 5中富集,MVP區(qū)域在program 4和program 5中富集。并且發(fā)現(xiàn),CT和MVP區(qū)域具有高度的增殖能力,與它們在program 3和program 5中的富集一致,這些program在細胞增殖中富集。MVP和PAN都富含組織重塑和免疫相互作用(program 4),而MVP更具有血管生成性,PAN更具有炎癥性。IT和LE都是富集了呼吸鏈復(fù)合物組裝通路的,LE是具有最強的呼吸鏈復(fù)合物組裝能力、但是具有較低的增殖能力,這解釋了它們在program 1中的富集。IT也在促炎免疫過程中富集,尤其是干擾素反應(yīng)(圖4 d)。綜上所述,BayesPrism可以將基因程序與空間解剖結(jié)構(gòu)聯(lián)系起來。
6.小結(jié)
在這里,作者通過開發(fā)一個嚴格的統(tǒng)計模型來整合scRNA-seq和bulk RNA-seq數(shù)據(jù),進行綜合分析也為疾病進展提供了新的見解。以GBM為例,聯(lián)合分析bulk RNA-seq隊列和空間解剖數(shù)據(jù),提出了一個將惡性細胞狀態(tài)和非惡性細胞浸潤與腫瘤進展聯(lián)系起來的模型(圖4e)。當惡性細胞快速生長時,它們會消耗營養(yǎng)物質(zhì),也可能遇到免疫壓力,導(dǎo)致壞死(圖4e)。與此相一致,觀察到免疫細胞和間充質(zhì)program 4的富集,在PAN區(qū)域顯示更強的間充質(zhì)激活和更低的呼吸活性(圖4b,c)。惡性細胞可能激活這些組織重塑通路,促進M2巨噬細胞極化和血管生成。隨著微血管結(jié)構(gòu)的發(fā)展,惡性細胞迅速增殖(圖4e),MVP附近惡性細胞的細胞周期評分較高(圖4d)。增殖細胞侵入鄰近的正常腦組織,那里的氧氣供應(yīng)充足。在此過程中,它們的主要任務(wù)從快速增殖轉(zhuǎn)變?yōu)楹粑饔?,以產(chǎn)生合成必要分子機制所需的ATP(圖4e)。這一發(fā)現(xiàn)是基于LE和IT結(jié)構(gòu)中呼吸通路的富集(圖4d)。最后,隨著局部氧水平的降低,惡性細胞恢復(fù)快速增殖。該研究還表明,program 3可能反映了血液供應(yīng)充足的癌癥生長的早期階段。綜上所述,該模型說明了GBM細胞如何重塑和響應(yīng)局部微環(huán)境的變化以適應(yīng)自身生長。與以前的方法相比,BayesPrism更準確地將bulk RNA-seq解析為細胞類型的比例和基因表達,從而深入了解腫瘤與微環(huán)境的相互作用。
二.Scissor算法介紹
下面我再具體介紹一下Scissor這項整合方法。為什么選擇這個方法呢?因為我覺得這個方法很具有創(chuàng)新性,作者開發(fā)了Scissor算法根據(jù)bulk RNA-seq數(shù)據(jù)的表型信息從單細胞數(shù)據(jù)識別特定的細胞亞群。在肺癌scRNA-seq數(shù)據(jù)集中,Scissor確定了與生存率較差和TP53突變相關(guān)的細胞亞群。在黑素瘤中,Scissor發(fā)現(xiàn)一個與免疫治療應(yīng)答相關(guān)的低PDCD1/CTLA4和高TCF7表達的T細胞亞群。
本篇文章發(fā)表在期刊: Nature Biotechnology 該期刊在最近一年的影響因子: 54.908,水平很高。
1. Scissor算法研發(fā)背景
單細胞測序技術(shù)使復(fù)雜組織細胞的綜合表征成為可能,從而使生物醫(yī)學研究和臨床實踐發(fā)生了革命性的變化。與測量整個組織平均特性的bulk數(shù)據(jù)相比,scRNA-seq允許在異質(zhì)組織生態(tài)系統(tǒng)中識別不同細胞亞群的細胞類型和狀態(tài)。要從單細胞數(shù)據(jù)中識別關(guān)鍵亞群,標準方法是執(zhí)行無監(jiān)督聚類來定義細胞群,并評估已知細胞類型和通路中標記基因的富集情況,以評估每個細胞集群的重要性。然而,識別驅(qū)動表型(如疾病階段、腫瘤轉(zhuǎn)移、治療反應(yīng)和生存結(jié)果)的細胞亞群具有不可缺少的重要性,因為它將促進細胞類型靶向治療和預(yù)后生物標志物的發(fā)現(xiàn)。Scissor的新穎之處在于,它使用來自bulk數(shù)據(jù)的表型信息來識別與疾病高度相關(guān)的細胞亞群,進而揭示疾病機制,提高疾病的診斷和治療。
2. Scissor算法基本流程
Scissor算法具體步驟:Scissor使用單細胞表達矩陣、bulk表達矩陣和感興趣的表型(圖5a)。每個樣本的表型注釋可以是一個連續(xù)的因變量、二元分類向量或臨床生存數(shù)據(jù)。Scissor的關(guān)鍵步驟是量化單細胞數(shù)據(jù)和批量數(shù)據(jù)之間的相似性,通過測量,如每對細胞和批量樣本的皮爾遜相關(guān)性。在這之后,剪刀優(yōu)化與樣本表型相關(guān)矩陣的回歸模型(圖5b)?;貧w模型的選擇取決于輸入表型的類型,例如,連續(xù)變量的線性回歸,二分變量的logistic回歸和臨床生存數(shù)據(jù)的Cox回歸。根據(jù)估計回歸系數(shù)的符號,系數(shù)為非零的細胞可表示為剪刀陽性(Scissor+)細胞和剪刀陰性(Scissor?)細胞,它們分別與感興趣的表型呈正相關(guān)和負相關(guān)(圖5c)。此外,為了控制單細胞和bulk數(shù)據(jù)之間的虛假關(guān)聯(lián),設(shè)計了一個可靠性顯著性檢驗,以確定所選數(shù)據(jù)是否適合的表型-細胞關(guān)聯(lián)(圖5d)。最后,Scissor選擇的細胞將在下游分析中進一步表征,如特征基因和功能富集途徑的探索(圖5e)。
首先評估了Scissor在一系列模擬數(shù)據(jù)集上的性能,結(jié)果表明Scissor識別的與已知的表型相關(guān)的細胞亞群與真實結(jié)果在很大程度上是一致的(圖5f, g, h)。
3.識別正常和腫瘤表型相關(guān)的細胞亞群
使用577個TCGA-LUAD bulk樣本中的腫瘤和正常表型對肺癌單細胞數(shù)據(jù)進行Scissor分析。在29,888個來自不同細胞類型的細胞中(圖5i), Scissor選擇了361個Scissor+細胞和534個Scissor?細胞,它們與腫瘤和正常表型具有高置信關(guān)系(圖5j)。正如預(yù)期,超過98%的Scissor+細胞被證實為惡性細胞(圖5k)。至于Scissor?細胞,主要分布于非惡性細胞類型中(圖5k)。髓系細胞和肺泡細胞是兩種主要選擇的細胞類型,分別占Scissor?細胞總數(shù)的42.3%和36.9%。剪刀細胞中的所有細胞類型,尤其是肺泡細胞,都是正常肺組織中的重要細胞類型。因此,這些分析證明了Scissor可以從單細胞數(shù)據(jù)中準確地識別出表型相關(guān)的細胞。
4. 發(fā)現(xiàn)低氧亞群與較差的生存
使用Scissor,基于帶有生存信息的471個TCGA-LUAD bulk樣本為識別肺癌scRNA-seq數(shù)據(jù)集中具有侵襲性的癌細胞亞群。這些細胞被分離成12個Cluster(圖6a)。在205個Scissor選擇的細胞中,201個Scissor+細胞與較差的存活率相關(guān)(定義為Scissor_WS細胞),只有4個Scissor?細胞與良好的存活率相關(guān)(圖6b)。Scissor_WS細胞為主要來自Cluster 1和Cluster 3(圖6c)。差異基因分析顯示,與其他所有細胞相比,Scissor_WS細胞中分別有23個上調(diào)基因(有多個重要的缺氧相關(guān)基因)和205個下調(diào)基因差異表達(圖6d,e)。功能富集分析也證實了缺氧相關(guān)的通路,如糖酵解和糖代謝通路在Scissor_WS細胞中被激活(圖6f)。并發(fā)現(xiàn)在獨立的GEO數(shù)據(jù)集中,特征分數(shù)高的患者預(yù)后明顯差于簽名分數(shù)低的患者(圖6g)。并在單因素Cox生存分析中,發(fā)現(xiàn)只有病理分期和Scissor_WS特征與患者生存顯著相關(guān)(圖6h)。此外,在對多變量Cox生存分析中的腫瘤分期進行調(diào)整后,Scissor_WS特征在兩個數(shù)據(jù)集中仍然具有統(tǒng)計學意義(圖6i)??傊?,Scissor算法從LUAD scRNA-seq數(shù)據(jù)中發(fā)現(xiàn)了一個具有侵襲性的癌細胞亞群,該亞群與較差的生存結(jié)果相關(guān),其特征是缺氧相關(guān)基因的過度表達。高度缺氧信號可能會推動LUAD進展,從而給腫瘤中含有大量此類細胞的患者帶來不良結(jié)果。
5. 分析與TP53突變相關(guān)的細胞亞群
基于TCGA-LUADTP53突變狀態(tài)(突變型或野生型)作為表型特征來指導(dǎo)對肺癌單細胞數(shù)據(jù)的細胞亞群的鑒定。Scissor共鑒定了414個與TP53突變相關(guān)的Scissor+細胞和318個與野生型相關(guān)的Scissor?細胞(圖7a)。差異基因分析顯示Scissor+細胞上調(diào)337個基因包括E2F靶基因和細胞周期進展相關(guān)基因,如AURKA、CDK1、CCNB2和TOP2A(圖7b)。功能富集分析也證實了E2F等細胞周期相關(guān)的通路在Scissor+細胞中被激活(圖7c)。轉(zhuǎn)錄因子分析顯示,在Scissor+細胞中,E2F轉(zhuǎn)錄因子家族成員E2F1和E2F4的活性均顯著升高(圖7d),然而TP53在Scissor+細胞中是失活的(圖7d)。此外,通過關(guān)聯(lián)這337個上調(diào)基因(定義為TP53突變特征)與臨床結(jié)果,證明TP53突變評分較高的患者預(yù)后明顯較差(圖7e)。該研究還發(fā)現(xiàn)MHC類相關(guān)基因HLA-A、B2M和CD74在Scissor+細胞中均下調(diào)(圖7f)。并已有報道稱在免疫治療耐藥的癌癥患者中B2M的功能喪失突變。因此,Scissor的表型分析表明TP53突變可能是檢查點抑制劑治療耐藥性的一個機制。
6.鑒定與免疫治療相關(guān)的T細胞亞群
為了了解免疫檢查點阻斷治療(ICB)反應(yīng)的機制,對黑素瘤scRNA-seq數(shù)據(jù)集中的T細胞以及已知免疫治療反應(yīng)信息的黑素瘤患者的bulk RNA-seq數(shù)據(jù)進行了Scissor分析,以識別與ICB反應(yīng)相關(guān)的T細胞亞群。在scRNA-seq數(shù)據(jù)分析中,這些T細胞被聚集成6個Cluster(圖8a)。通過使用Scissor算法,確定了105個T細胞為Scissor+細胞,這些細胞與良好的免疫治療反應(yīng)相關(guān)(定義為Scissor_FR細胞)(圖8b),這105個Scissor_FR細胞主要分布在Cluster 2和Cluster 3中(圖8c)。差異基因分析表明Scissor_FR細胞增加了與T細胞記憶相關(guān)基因(CCR7, SELL和IL7R)的表達以及低表達抑制基因(HAVCR2, LAG3, PDCD1和CTLA4)和MHC II類基因(HLA-DRB5, HLA-DRB1, HLA-DPA1, HLA-DQB2和HLA-DRB6)(圖8d,e)。同時,Scissor_FR細胞也表現(xiàn)出轉(zhuǎn)錄因子TCF7的表達增強,這與ICB治療的良好結(jié)果相關(guān)(圖8e)。此外,富集分析顯示,Scissor_FR細胞TNF-α信號通路較高,CTLA4、PD1信號通路活性較低(圖8f)。以上與有效免疫治療應(yīng)答相關(guān)的差異表達基因(定義為免疫治療應(yīng)答特征)可以用于預(yù)測治療應(yīng)答率。該研究發(fā)現(xiàn)有ICB應(yīng)答者的簽名分數(shù)明顯高于無ICB應(yīng)答者(圖8g)。此外,免疫治療應(yīng)答信號中上調(diào)和下調(diào)的基因在應(yīng)答者和非應(yīng)答者中也顯著富集(圖8h)。進一步評估了五種不同分化狀態(tài)的腫瘤浸潤淋巴細胞的免疫治療反應(yīng)特征,發(fā)現(xiàn)LAG3-low/PD1-low效應(yīng)CD8 T細胞和記憶前體CD8 T細胞的特征評分最高(圖8i,j)。這一結(jié)果表明,Scissor_FR細胞更像PD1-low記憶前體細胞,具有較高的TCF7表達,并與良好的免疫治療反應(yīng)有關(guān)??傊?,對黑色素瘤scRNA-seq數(shù)據(jù)集的Scissor分析獨立揭示了PDCD1/CTLA4低和TCF7高T細胞亞群,其獨特的轉(zhuǎn)錄組特征對免疫治療的良好反應(yīng)至關(guān)重要。
三.小編總結(jié)
通過使用BayesPrism算法整合,既可以解析大樣本隊列的細胞類型的比例,又能獲得細胞特異的基因表達,通過使用Scissor算法整合,可以獲得與臨床表型緊密相關(guān)的細胞群體??偠灾?,單細胞轉(zhuǎn)錄組數(shù)據(jù)和bulk轉(zhuǎn)錄組數(shù)據(jù):優(yōu)勢互補,將單細胞轉(zhuǎn)錄組數(shù)據(jù)和bulk轉(zhuǎn)錄組數(shù)據(jù)進行有效結(jié)合,將以全新研究思路出發(fā),會發(fā)現(xiàn)更多未知且精細化結(jié)果。