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