圖1?利用Monocle推斷類風(fēng)濕性關(guān)節(jié)炎PBMCs B細(xì)胞的分化軌跡[3]

2020年Nature Communication發(fā)表的Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy[4]，利用RNA速度表征了膠質(zhì)母細(xì)胞瘤的分化過程，證實(shí)了膠質(zhì)母細(xì)胞瘤層次結(jié)構(gòu)頂點(diǎn)是腫瘤祖細(xì)胞，星形細(xì)胞、間充質(zhì)、少突膠質(zhì)細(xì)胞和神經(jīng)元癌細(xì)胞比腫瘤祖細(xì)胞的分化程度更高，并且在每個(gè)患者中均可見特定細(xì)胞分化方向（圖2）。

圖2?基于RNA?velocity分析膠質(zhì)母細(xì)胞瘤細(xì)胞分化路徑[4]

單細(xì)胞測序數(shù)據(jù)分析之2、細(xì)胞通訊分析–解析細(xì)胞間的信號(hào)通訊關(guān)系

多細(xì)胞生物是由很多不同類型細(xì)胞組成的開放而復(fù)雜體系，配體受體復(fù)合物介導(dǎo)的細(xì)胞間通訊對(duì)協(xié)調(diào)發(fā)育、分化和炎癥等多種生物學(xué)過程至關(guān)重要。細(xì)胞通訊分析，又稱細(xì)胞受體-配體互作分析，是以細(xì)胞亞群的基因表達(dá)量數(shù)據(jù)為研究對(duì)象，通過獲得細(xì)胞中配體及受體基因的表達(dá)信息，比較細(xì)胞類型之間的配體與受體基因表達(dá)差異，分析得到細(xì)胞亞群間的信號(hào)通訊關(guān)系，在闡明生物學(xué)過程中細(xì)胞間通訊的復(fù)雜性、多樣性和動(dòng)態(tài)性方面有重要意義。

2020年Nature Communications發(fā)表的Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma[5]，對(duì)8例膀胱癌腫瘤樣本和3例配對(duì)癌旁粘膜樣本進(jìn)行了單細(xì)胞轉(zhuǎn)錄組測序分析，利用CellphoneDB[6]軟件進(jìn)行細(xì)胞通訊分析，分析結(jié)果發(fā)現(xiàn)，炎性癌癥相關(guān)成纖維細(xì)胞iCAFs與其他細(xì)胞類型的相互作用最多，且與內(nèi)皮細(xì)胞ECs的相互作用最強(qiáng)（圖3a）；iCAF表達(dá)更高水平的CXCL12，其受體包括 DPP4、CXCR3、CXCR4 和 ACKR3 (CXCR7)，由于CXCR4和CXCR3 在免疫細(xì)胞上廣泛表達(dá)，猜測iCAFs分泌 CXCL12是膀胱癌免疫浸潤狀態(tài)的原因（圖3b）；iCAFs產(chǎn)生VEGF（包括VEGFA和VEGFB），與ECs的VEGF受體(FLT1、KDR、MET和FLT4)結(jié)合，促進(jìn)血管生成，與腫瘤的惡性進(jìn)展有關(guān)（圖3c）。

圖3?利用CellPhoneDB分析膀胱癌不同細(xì)胞類型間配體-受體相互作用[6]

單細(xì)胞測序數(shù)據(jù)分析之3、GSEA/GSVA分析–基于功能基因集的富集分析策略

GSEA( Gene Set Enrichment Analysis)，是2005年由Broad Institute研究開發(fā)的一種基于基因集的富集分析方法[7]，用來評(píng)估一個(gè)預(yù)先定義的基因集的基因在與表型相關(guān)度排序的基因表中的分布趨勢，從而判斷其對(duì)表型的貢獻(xiàn)。GSEA是從所有基因的表達(dá)豐度出發(fā)，分析在不同的通路中的基因的整體表達(dá)影響，這也是區(qū)別于GO/KEGG富集分析的地方，GSEA不需要設(shè)定差異閾值篩選目標(biāo)基因集，理論上更容易囊括細(xì)微但協(xié)調(diào)性的變化對(duì)生物通路的影響。

上述膀胱癌的單細(xì)胞文獻(xiàn)中[5]，利用GSEA對(duì)iCAFs和肌癌相關(guān)成纖維細(xì)胞mCAF兩種細(xì)胞類型的基因進(jìn)行了GSEA富集分析，結(jié)果顯示，iCAFs與胞外基質(zhì)降解基因集顯著相關(guān)，暗示iCAFs可能在腫瘤轉(zhuǎn)移中發(fā)揮作用，同時(shí)iCAFs與細(xì)胞因子-細(xì)胞因子受體相互作用基因集也顯著相關(guān)（圖4，左圖）；mCAFs則與肌肉收縮、PGC1A通路顯著相關(guān)，這一結(jié)果與之前的體外研究結(jié)果一致（圖5，右圖）。

圖4?利用GSEA分析膀胱癌腫瘤相關(guān)成纖維細(xì)胞顯著富集的通路[5]

GSVA分析（Gene Set Variation Analysis）[8]與GSEA（Gene Set Enrichment Analysis）都是一種基于基因集的富集分析方法，能夠有效彌補(bǔ)傳統(tǒng)富集分析對(duì)微效基因的有效信息挖據(jù)不足等問題，更為全面地對(duì)某一功能單位的調(diào)節(jié)作用進(jìn)行解釋，Broad研究所在GSEA發(fā)布8年之后，開發(fā)了GSVA算法來拓展基因集分析的應(yīng)用。GSEA分析主要用于兩兩組間比較的方案設(shè)計(jì)中，對(duì)于分組比較復(fù)雜的方案設(shè)計(jì)則比較適合GSVA分析，GSVA不需要預(yù)先進(jìn)行樣本之間的差異分析，依據(jù)表達(dá)矩陣就可以計(jì)算每個(gè)樣本中特定基因集的變異分?jǐn)?shù)。

2020年Cell Research發(fā)表的Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma[9]，利用GSVA分析了鼻咽癌NPC不同細(xì)胞亞型間的通路活性差異，結(jié)果顯示，參與干擾素 (IFN)-α和IFN-γ通路在巨噬細(xì)胞中相對(duì)上調(diào)，而血管生成、NF-κB 介導(dǎo)的 TNF-α 信號(hào)傳導(dǎo)和缺氧相關(guān)信號(hào)通路在單核細(xì)胞中上調(diào)（圖5，左圖）；生發(fā)中心（GC）B細(xì)胞表現(xiàn)出丙酮酸代謝和MYC途徑的激活，而漿細(xì)胞中的糖酵解途徑被上調(diào)，在漿細(xì)胞中也檢測到巨噬細(xì)胞趨化途徑的高活性，表明這些途徑在促進(jìn)巨噬細(xì)胞募集到 NPC 中的 TME 中的潛在作用（圖5，右圖）。

圖5?利用GSVA分析鼻咽癌不同髓系細(xì)胞亞型間的通路活性差異[9]

單細(xì)胞測序數(shù)據(jù)分析之4、腫瘤拷貝數(shù)變異分析–揭示惡性細(xì)胞表型

拷貝數(shù)變異(Copy number variation, CNV)是由基因組發(fā)生重排而導(dǎo)致的，基因組大片段的拷貝數(shù)增加或者減少，基因組結(jié)構(gòu)變異(Structural variation, SV) 的重要組成部分，也是人類疾病的重要致病因素之一。與正常細(xì)胞相比，腫瘤基因組部分區(qū)域呈現(xiàn)過表達(dá)或低表達(dá)狀態(tài)，通過與一組參考的“正常”細(xì)胞相比，比較不同樣本間或不同細(xì)胞類型之間的CNV基因表達(dá)差異，探索腫瘤基因組位置上基因的表達(dá)強(qiáng)度，最終反映基因大片段區(qū)域的CNV事件，鑒定體細(xì)胞整個(gè)染色體或大片段染色體的擴(kuò)增或缺失。

2021年Cell death discovery發(fā)表的Single-cell RNA transcriptome reveals the intra-tumoral heterogeneity and regulators underlying tumor progression in metastatic pancreatic ductal adenocarcinoma[10]，利用inferCNV[11]對(duì)胰腺導(dǎo)管腺癌轉(zhuǎn)移灶的5種導(dǎo)管細(xì)胞亞型進(jìn)行inferCNV分析，鑒別惡性細(xì)胞，分析結(jié)果發(fā)現(xiàn)，五種導(dǎo)管細(xì)胞亞型都經(jīng)歷了顯著的CNV事件，均為轉(zhuǎn)移灶中的惡性細(xì)胞（圖6E）；計(jì)算CNV scores并繪制小提琴，發(fā)現(xiàn)1型導(dǎo)管細(xì)胞的CNV水平顯著高于其他類型的導(dǎo)管細(xì)胞，提示1型導(dǎo)管細(xì)胞惡性程度更高（圖6F）。

圖6?利用inferCNV揭示胰腺導(dǎo)管腺癌惡性細(xì)胞類型[11]

單細(xì)胞測序數(shù)據(jù)分析之5、轉(zhuǎn)錄因子活性分析–挖掘關(guān)鍵調(diào)控轉(zhuǎn)錄因子

?單細(xì)胞研究通常會(huì)涉及到一個(gè)核心關(guān)鍵問題：細(xì)胞的異質(zhì)性以及這種異質(zhì)性是如何發(fā)展和維持的。這種細(xì)胞異質(zhì)性很大程度上是由潛在的基因調(diào)控網(wǎng)絡(luò)決定的，特定轉(zhuǎn)錄因子（transcription factor，TF）集合的協(xié)同表達(dá)驅(qū)動(dòng)各自靶標(biāo)基因的表達(dá)，從而建立特定的基因表達(dá)譜。因此，單細(xì)胞的基因調(diào)控網(wǎng)絡(luò)對(duì)于深入挖掘細(xì)胞異質(zhì)性背后的生物學(xué)意義是至關(guān)重要的。SCENIC Suite[12]是一套基于SCENIC（Single-Cell Regulatory Network Inference and Clustering）來研究和破譯基因調(diào)控的工具，能從單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)中推斷TF、基因調(diào)控網(wǎng)絡(luò)和細(xì)胞類型，基本原理是基于共表達(dá)和DNA調(diào)控保守序列（motif）分析推斷基因調(diào)控網(wǎng)絡(luò)，然后在每個(gè)細(xì)胞中分析網(wǎng)絡(luò)活性以鑒定細(xì)胞狀態(tài)。

上述鼻咽癌的單細(xì)胞文獻(xiàn)[9]，利用SCENIC分析髓系細(xì)胞的轉(zhuǎn)錄因子活性，發(fā)現(xiàn)BNR1H3 和 TFEC 轉(zhuǎn)錄因子活性在巨噬細(xì)胞中顯著上調(diào)（圖7e）；進(jìn)一步利用生存分析揭示了NR1H3和TFEC的高表達(dá)水平與 NPC 患者預(yù)后改善之間的關(guān)聯(lián)（圖k, n），暗示NR1H3和TFEC在鼻咽癌發(fā)展中的重要調(diào)控作用。

圖7??基于SCENIC分析鼻咽癌關(guān)鍵調(diào)控轉(zhuǎn)錄因子[9]

單細(xì)胞測序數(shù)據(jù)分析之6、加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析–篩選表型相關(guān)核心調(diào)控網(wǎng)絡(luò)

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（Weighted Gene Co-expression Network Analysis, WGCNA）[13]是構(gòu)建基因共表達(dá)網(wǎng)絡(luò)的常用方法，可以探索模塊與特定表型或疾病的關(guān)聯(lián)關(guān)系,最終達(dá)到鑒定基因網(wǎng)絡(luò)的目的；單細(xì)胞測序技術(shù)可以揭示特定腫瘤組織中的細(xì)胞特異性，對(duì)細(xì)胞進(jìn)行分類，并且識(shí)別特定的標(biāo)志物，但其檢測的細(xì)胞數(shù)量和病例來源都是有限的。利用WGCNA分析單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)，可以提供一套有別于高変基因、差異分析的方法，不依賴于數(shù)據(jù)庫直接用表達(dá)量的相關(guān)性值預(yù)測調(diào)控關(guān)系，篩選某些細(xì)胞亞群中有關(guān)聯(lián)作用的基因集（稱為模塊），可以從成千上萬的基因中挑選出高度相關(guān)的基因的模塊，并將模塊與表型進(jìn)行關(guān)聯(lián)，尋找marker gene或治療靶點(diǎn)。

2021年Nature communications發(fā)表的Cellular and molecular landscape of mammalian sinoatrial node revealed by single-cell RNA sequencing[14]，利用WGCNA分析成年小鼠竇房結(jié)SAN細(xì)胞基因表達(dá)網(wǎng)絡(luò)，確定了十個(gè)共表達(dá)模塊（圖8a）；功能富集分析結(jié)果顯示有六個(gè)模板與離子通道的調(diào)節(jié)有關(guān)，其中紅色模塊基因主要集中在K+和Ca2+通道上（圖8c）。

圖8?利用WGCNA分析小鼠竇房結(jié)細(xì)胞基因共表達(dá)調(diào)控網(wǎng)絡(luò)[14]

?看完上述內(nèi)容，單細(xì)胞測序分析思路與方法get起來，屬于你的單細(xì)胞故事將由你自己講述~

百邁客擁有專業(yè)且強(qiáng)大的單細(xì)胞生信分析團(tuán)隊(duì)，在項(xiàng)目實(shí)戰(zhàn)中積累了大量的數(shù)據(jù)挖掘經(jīng)驗(yàn)，能提供但不限于以上內(nèi)容的高級(jí)分析，歡迎點(diǎn)擊下方按鈕聯(lián)系我們。

圖9 ?百邁客部分分析實(shí)例展示

 

參考文獻(xiàn)：

[1]?http://cole-trapnell-lab.github.io/monocle-release/docs/

[2]?La Manno, Gioele et al. “RNA velocity of single cells.” Nature vol. 560,7719 (2018): 494-498. doi:10.1038/s41586-018-0414-6

[3]?Wu, Xunyao et al. “Single-cell sequencing of immune cells from anticitrullinated peptide antibody positive and negative rheumatoid arthritis.” Nature communications vol. 12,1 4977. 17 Aug. 2021, doi:10.1038/s41467-021-25246-7

[4]?Bergen, Volker et al. “Generalizing RNA velocity to transient cell states through dynamical modeling.” Nature biotechnology vol. 38,12 (2020): 1408-1414. doi:10.1038/s41587-020-0591-3

[5]?Chen, Zhaohui et al. “Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma.” Nature communications vol. 11,1 5077. 8 Oct. 2020, doi:10.1038/s41467-020-18916-5

[6]?Subramanian, Aravind et al. “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.” Proceedings of the National Academy of Sciences of the United States of America vol. 102,43 (2005): 15545-50. doi:10.1073/pnas.0506580102

[7]?Chen, Zhaohui et al. “Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma.” Nature communications vol. 11,1 5077. 8 Oct. 2020, doi:10.1038/s41467-020-18916-5

[8]?H?nzelmann, Sonja et al. “GSVA: gene set variation analysis for microarray and RNA-seq data.” BMC bioinformatics vol. 14 7. 16 Jan. 2013, doi:10.1186/1471-2105-14-7

[9]?Chen, Yu-Pei et al. “Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma.” Cell research vol. 30,11 (2020): 1024-1042. doi:10.1038/s41422-020-0374-x

[10]?Xu, Qianhui et al. “Single-cell RNA transcriptome reveals the intra-tumoral heterogeneity and regulators underlying tumor progression in metastatic pancreatic ductal adenocarcinoma.” Cell death discovery vol. 7,1 331. 3 Nov. 2021, doi:10.1038/s41420-021-00663-1

[11]?https://github.com/broadinstitute/inferCNV

[12]?https://scenic.aertslab.org/

[13]?Langfelder, Peter, and Steve Horvath. “WGCNA: an R package for weighted correlation network analysis.” BMC bioinformatics vol. 9 559. 29 Dec. 2008, doi:10.1186/1471-2105-9-559

[14]?Liang, Dandan et al. “Cellular and molecular landscape of mammalian sinoatrial node revealed by single-cell RNA sequencing.” Nature communications vol. 12,1 287. 12 Jan. 2021, doi:10.1038/s41467-020-20448-x