基于整合生物信息學分析識別骨肉瘤干性相關基因

符 晴,王思喬,缐述源,黃潤之,張 杰

(1.同濟大學醫學院,上海 200092; 2.上海長海醫院燒傷外科,上海 200433)

具有高度自我更新和分化能力的癌癥干細胞(cancer stem cells,CSCs)被認為是癌癥侵襲和轉移的根源,但其潛在機制尚未被完全闡明[6]。此外,患者治療反應和生存時間的差異提示了OS腫瘤細胞的異質性。惡性腫瘤的異質性是腫瘤診治的關鍵和難點,不同基因型及表型的惡性腫瘤細胞表現出不同的生物學行為。最新的研究發現OS中癌癥干細胞與患者耐藥性和腫瘤發生有關[7]?；趍RNA測定的干細胞指數(mRNA expression-based stemness index,mRNAsi)是衡量腫瘤細胞與干細胞之間相似性的指標,可以被用來量化腫瘤干細胞,亦稱為腫瘤干性指數。大量研究表明,mRNAsi在多種腫瘤中被作為評價患者生存預后和腫瘤分型的有效指標,如結直腸癌、乳腺癌、頭頸部鱗狀細胞癌等[8-10]。然而目前仍然沒有研究基于mRNAsi來分析OS的腫瘤干性。近年來,單細胞測序(single cell RNA sequencing,scRNA-seq)技術的發展為腫瘤細胞異質性的分析提供了新方向[11]。單細胞測序可從基因組學、轉錄組學、表觀遺傳學以及多組學對腫瘤微環境進行分析,能進一步揭示腫瘤內部的細胞異質性,研究腫瘤的發生、發展及轉移過程,闡明腫瘤患者耐藥機制,繪制腫瘤微環境的單細胞圖譜[12-13]。本研究采用整合生物信息學分析的方法,識別轉移性骨肉瘤中干性相關基因(stemness-related genes,SRGs),在單細胞水平深入地探究OS腫瘤細胞的異質性,為實現轉移性OS分子治療的轉化醫學研究提供理論依據。

1 材料與方法

1.1 數據來源及其預處理

OS患者的RNA測序(RNA sequencing,RNA-seq)數據和臨床信息(年齡、性別、種族、存活時間和轉移狀態)來源于TARGET數據庫(https:∥ocg.cancer.gov/programs/TARGET)。在GEO數據平臺上(https:∥www.ncbi.nlm.nih.gov/gds)下載OS單細胞測序數據文件,GPL24676(平臺注釋文件)與矩陣文件GSE152048(包含11個OS患者樣本)。刪除臨床信息不完整以及不符合納入、排除標準的OS患者數據,最終篩選出87例OS患者,其中65例為原發性OS患者,其余22例為轉移性OS患者。數據截至2020年3月。為了減少測序誤差和批次效應,用微陣列分析線性模型軟件包來校正已篩選的RNA-seq數據[14]。在R3.5.1環境下,將OS樣本的轉錄組測序結果合并為1個矩陣,使用“gen-code.v22.annotation.gene.probeMap”文件注釋轉錄組測序數據,將基因名從Ensemble ID轉換為基因符號,“edgeR”包對轉錄組測序的Counts數據進行標準化。

1.2 獲取基于mRNA表達的干性指數

干性指數mRNAsi與腫瘤去分化程度(腫瘤干性)密切相關。并且mRNAsi值越高,患者腫瘤干細胞活性越高,發生轉移的可能性也越大。為了估計OS患者的腫瘤干性,基于OS患者的mRNA表達水平,使用一類邏輯回歸(one-class logistic regre-ssion,OCLR)機器學習算法計算上述87個OS患者的mRNAsi[15]。

1.3 干性指數關聯性分析

皮爾遜關聯性分析被用于篩選出與mRNAsi明顯線性相關的基因[16]。mRNAsi(記為X)與基因(記為Y)之間的皮爾遜相關系數r定義為二者之間的協方差和標準差的商。二者協方差的計算公式如下:

其中E(X)和E(Y)分別代表二者的期望,也就是平均值。方差的計算公式如下:

最后通過計算相關系數r來分析二者之間的關聯性及關聯強度。計算公式如下:

相關系數r>0.5并且具有統計學意義(P<0.05)的基因被定義為干性相關基因(Cor-sig gene),其高表達水平可能與腫瘤干細胞的自我更新能力與惡性生物學行為密切相關。

1.4 Cox回歸分析

Cox回歸分析是一種半參數分析方法,用于研究各種因素對患者生存期長短的影響。本研究使用R軟件包中的“Survival”功能,基于OS患者的表達譜,應用單變量Cox回歸分析估計每個基因的表達水平高低對患者生存時間長短的影響。風險函數用h(t)表示,計算公式如下:

表示在t時刻OS患者死亡的瞬時發生率。隨后進一步構建Cox比例風險模型,計算公式如下:

lnh(t)=lnh0(t)+β1X1+…+βiXpi+…+βpXp

其中Xi是協變量,反映基因i的表達水平,βi是回歸系數,h0是基準風險,即當所有Xi=0時的事件風險。隨后分析特定基因的風險比(hazard ratio,HR),計算公式如下:

HR=eβi

如果值βi>0,即風險比>1,表示隨著基因表達水平的增加,患者死亡的風險也逐漸增加,于是生存期減少。HR>1表示該基因的表達水平高是危險因素,可能會減少患者的生存時間。單變量Cox回歸分析中HR>1并且生存時間差異具有統計學意義(P<0.05)的基因被進一步納入多變量Cox回歸分析。分析的變量包括HR,性別、年齡、原發部位、轉移狀態。多變量Cox回歸分析中預后價值顯著(P<0.05)的基因被定義為Cox風險相關基因(Cox-sig gene)。

1.5 Kaplan-Meier(KM)生存分析

KM生存分析方法即乘積極限法,由Kaplan和Meier于1958年提出,因此被命名為Kaplan-Meier法或簡稱為KM法[17]。作為當前生存分析最常用的非參數方法,KM法根據臨床數據對患者的生存時間進行分析和推斷,研究生存時間和終點事件(如死亡事件)與諸多影響因素之間的關聯程度。臨床數據有兩種類型。(1) 完全數據:所有觀察對象的生存時間或發生終點事件的具體時間都被完整記錄;(2) 刪失數據:觀察對象在研究結束前發生了研究之外的其他事件或生存結局(失訪、退出、終止等),無法明確記錄到發生終點事件的生存時間。

為研究基因的表達水平與OS患者的生存相關性,根據患者特定基因表達水平將OS患者分為高表達組和低表達組,在每個OS患者發生死亡事件的時間點上,進行生存概率的計算(觀察終點:12年)。其中生存概率指某時間段(tk)開始時存活的OS患者至該時間段結束時仍然存活的可能性大小,以p(tk)表示,計算公式如下:

p(tk)=tk結束時OS患者存活例數/

同組tk初始OS患者例數

由于存在刪失數據,本研究假定刪失事件在觀察時間內各個時間點等機會發生,分母改用校正觀察例數,計算公式如下:

校正觀察例數=觀察初始OS患者例數-觀察期內刪失OS患者例數/2

隨后計算生存率(survival rate),以S(tk)表示,指OS患者于tk結束后仍存活的概率,計算公式如下:

S(tk)=P(T>tk)=p(1)×p(2)×…×p(tk)

本研究的KM生存曲線以生存時間為X軸,以生存率為Y軸,使用R軟件包中的“Survival”功能繪制。KM生存曲線展示了生存率隨時間變化的趨勢,其斜率越高表示患者生存率下降越快。最后以對數秩檢驗法分析兩組曲線之間的差異。兩組總體生存曲線差異有統計學意義(P<0.05)的基因被定義為KM生存相關基因(KM-sig gene)。

1.6 繪制韋恩圖(Venn diagram)

基于上述分析結果,Cor-sig gene、Cox-sig gene、KM-sig gene三者交集內的基因被定義為關鍵基因。在線繪圖工具jvenn被用以對結果進行可視化(http:∥www.bioinformatics.com.cn/)[18]。

1.7 scRNA-seq數據標準化

將scRNA-seq數據進行前期的整理和校正,對重復基因的表達水平取平均值。采用“Seurat”R軟件包中的“NormalizeData”算法對基因表達矩陣進行標準化。使用“harmony”R軟件包去除數據的批次效應。使用“Percentage Featureset”函數來計算線粒體基因的百分比,并以此為標準來過濾數據(線粒體基因的百分比<5%,線粒體基因數量>50)。提取細胞間變異系數較大的基因進行后續分析。

1.8 主成分分析(principal component analysis,PCA)

使用“Seurat”R軟件包中的“FindVariable-Features”算法識別高度可變基因?；诟叨瓤勺兓虻谋磉_水平,使用“Seurat”R軟件包中的“RunPCA”算法進行PCA,選擇特征根>1的前15個主成分繪圖,其中每個主成分包含15個基因。

1.9 聚類分析

用R語言進行t分布隨機鄰接嵌入(t-distributed stochastic neighbor embedding,t-SNE)聚類分析,獲得不同細胞亞群的聚類。優化標準模塊化,進行細胞聚類的可視化。識別差異表達的特征性基因,確定每個聚類的marker基因。

1.10 細胞亞群類型注釋

采用“Seurat”R軟件包識別各個細胞亞群中的特異性上調基因,作為對應細胞亞群的潛在標記基因。參照CellMarker數據庫(http:∥xteam.xbio.top/CellMarker/)中的人類細胞標志基因或已有文獻報道的細胞標志物,結合已識別的潛在標志基因,排除算法自動注釋的偏倚,對細胞亞群進行注釋。

1.11 細胞周期與細胞通信分析

為了預測OS細胞的細胞周期狀態,基于既往文獻報道的43個G1/S周期和54個G2/M周期相關的關鍵基因,本研究采用“CellCycleScore”函數對細胞周期進行判定。采用ggplot2(https:∥github.com/tidyverse/ggplot2)軟件對細胞周期分析結果進行可視化?；诨虮磉_矩陣,結合已有的配體-受體信息,采用“iTALK”軟件包量化配體-受體相互作用的強度,推測細胞間的互作關系。

分析流程見圖1。

圖1 本研究的分析流程圖

2 結 果

2.1 干性指數的計算

本研究使用一類線性回歸算法(one class linear regression,OCLR)獲得每個OS樣本的干細胞指數(mRNA expression-based stemness index,mRNAsi),實現對87個OS樣本腫瘤干性的量化。87個OS患者的mRNAsi與臨床指標之間(存活狀態、性別、種族、是否轉移、原發灶位置等)的整體關系見圖2。不同存活狀態(存活vs死亡)、不同性別(男性vs女性)、不同種族(白人vs亞洲人vs黑人vs未知人種)、不同轉移狀態(轉移vs無轉移),以及不同原發部位(下肢vs骨盆vs上肢)OS樣本的mRNAsi之間差異均無統計學意義(P>0.05)。

圖2 骨肉瘤患者基線特征及mRNAsi水平

腫瘤干細胞(cancer stem cells,CSCs)是指擁有與正常干細胞相關特征的癌細胞,能在特定癌癥樣本中產生所有細胞類型。通常腫瘤干細胞被認為有形成腫瘤的潛力,特別是隨著癌癥的轉移,能產生新的癌癥病灶。mRNAsi值越高,患者腫瘤干細胞活性越高,發生轉移的可能性也越大。

2.2 關鍵基因

韋恩圖顯示,皮爾遜相關性分析篩選出5 131個干性相關基因;Cox回歸分析篩選出1 640個風險相關基因;KM生存分析篩選出7 066個生存相關基因。其中658個基因在3個基因集中都有統計學意義,被定義為關鍵基因,納入后續分析,見圖3。

圖3 關鍵基因的韋恩圖

此外,根據骨肉瘤干性基因與mRNAsi之間的相關系數R,將R>0.6或R<-0.5作為篩選標準,進一步篩選與腫瘤干性高度相關的關鍵基因,最終獲得了57個關鍵基因,其與mRNAsi之間的相關系數見圖4。

圖4 腫瘤干性相關基因與mRNAsi之間的相關系數

在關鍵基因中,NADH脫氫酶(泛醌)1β亞復合物9(NADH:Ubiquinone Oxidoreductase Subunit B9,NDUFB9)與mRNAsi之間的正相關性最強(相關系數r=0.57,P<0.001),并且NDUFB9表達水平高的OS患者預后更差(P=0.007),見圖5。

圖5 NDUFB9與mRNAsi的相關性及預后價值

相反,關鍵基因EHD家族蛋白2(EH domain containing 2,EHD2)與mRNAsi之間的負相關性最強(相關系數r=-0.43,P<0.001),并且EHD2表達水平低的OS患者預后更差(P=0.002),見圖6。

圖6 EHD2與mRNAsi的相關性及預后價值

2.3 OS腫瘤組織內的主要細胞亞群及其分子特征

經過嚴格質控,本研究最終共獲得123 307個細胞的高質量scRNA-seq數據?；赗語言,進行PCA并繪制熱圖,識別了15類特征根>1的主成分。各基因分組及表達差異顯著,提示OS腫瘤細胞表現出顯著分離趨勢。在11個OS患者的腫瘤組織中,使用tSNE聚類分析將15個主成分中的細胞分為10種細胞亞群,見圖7A,分別為成軟骨細胞(chondroblasts)、內皮細胞(endothelial cells)、成纖維細胞(fibroblasts)、淋巴細胞(lymphocytes)、間充質干細胞(mesenchymal stem cells,MSCs)、髓系細胞(myeloid cells)、成肌細胞(myoblasts)、成骨細胞(osteoblasts)、破骨細胞(osteoclasts)和周細胞(pericytes)。病理學分析發現,大部分單細胞來源于傳統骨肉瘤分型,見圖7B。

上述10種細胞亞群在不同OS患者腫瘤組織中的數目和比例具有明顯差異,見圖7C。細胞特征圖顯示分化簇24(cluster of differentiation 24,CD24)主要表達于成骨細胞、間充質干細胞及成肌細胞;分化簇44(cluster of differentiation 44,CD44)主要表達于成纖維細胞和破骨細胞;MKI67主要表達于成骨細胞和間充質干細胞;分化簇133(cluster of differe-ntiation 133,CD133或Prominin 1,PROM1)主要表達于成骨細胞和周細胞;NDUFB9在多種OS腫瘤細胞亞群中表達水平上調,包括間充質干細胞和成骨細胞等;而EHD2在各腫瘤細胞亞群的表達量都較低,見圖7D。CD24是一種小分子量、高度糖基化的細胞膜蛋白質,通過參與腫瘤發生發展相關信號通路(如Wnt信號通路和MAPK信號通路)調控腫瘤細胞的生長、增殖和轉移,與腫瘤的侵襲轉移密切相關。CD44是一種非激酶跨膜糖蛋白,可促進腫瘤細胞增殖和侵襲,在大多數肉瘤中高表達,大量研究證實其對腫瘤進展、轉移和耐藥性有直接影響。此外,CD44已被證實是骨肉瘤診斷和預后評價的有效標志物[19]。CD44也是檢測和分離化療抗性CSCs的特異性生物標志物[20]。MKI67表達水平反映細胞增殖狀態,在癌癥轉移、侵襲和進展中起著至關重要的作用,是癌癥的預后預測因子[21-22]。PROM1編碼五跨膜糖蛋白,通常在癌細胞中過表達,通過抑制分化來維持腫瘤細胞的干細胞特性[23]。PROM1在各類腫瘤干細胞的鑒定和分選中被廣泛應用。通過對OS腫瘤組織進行單細胞測序發現,本研究鑒定的關鍵基因NDUFB9在具有干細胞特征的細胞亞群中表達上調,進一步明確了上述基因與OS腫瘤干性特征密切相關,可作為潛在預后標志物。

為了進一步解析OS腫瘤細胞的異質性,基于CellMarker數據庫中既往文獻報道的細胞標志物,本研究識別并注釋了10個細胞亞群,包括成骨細胞(MMP9、ACP5、CTSK)、破骨細胞(MT1X、ALPL、CYR61)、成纖維細胞(MMP13、IBSP、SPP1)、間充質干細胞(SFRP2、COMP、PLA2G2A)、髓系細胞(HLA-DRA、CCL3、CCL3L1)、周細胞(FGFBP2、COL11A1、HAPLN1)、內皮細胞(SPARCL1、RAMP2、PLVAP)、成肌細胞(SFTPC、MYLPF、PTN)、淋巴細胞(CCL5、CD69、IL32)、成軟骨細胞(SSPN、COL 12A1、LRP1),見圖7E。差異表達分析從10個亞群中鑒定出59 578個標志基因,各細胞亞群基因表達差異明顯,見圖7F。

2.4 細胞通信和細胞周期分析

克利夫蘭圖顯示,4種經典Marker在10種細胞亞群中表達水平有明顯差異。如前所述,上述10種細胞亞群在11個OS患者腫瘤組織中的比例具有顯著差異,見圖8A。使用“Seurat”包中的“Find-AllMarkers”函數識別不同細胞亞群的標記基因,并通過熱圖展示這些細胞亞群中前5個差異表達標志基因的表達水平,見圖8B。細胞通信分析顯示,成軟骨細胞和MSCs占據了細胞通信網絡的中心,主要通過配體:膠原蛋白Ⅰ型α 1鏈(collagen type Ⅰ Alpha 1 chain,COL1A1)和膠原蛋白Ⅰ型α 2鏈(collagen type Ⅰ Alpha 2 chain,COL1A2)發送信號,并通過受體低密度脂蛋白受體相關蛋白1(low-density lipoprotein receptor-related protein 1,LRP1)接收信號,見圖8C、D。細胞周期分析顯示淋巴細胞、MSCs和成骨細胞表現為細胞增殖較為活躍,主要處于G2/M期;而周細胞和成肌細胞主要處于相對靜息的S期,見圖8E。而NDUFB9在MSCs和成骨細胞等增殖較為活躍的細胞亞群中的表達水平較高,并且MSCs在OS腫瘤微環境細胞間通訊中占據著中心樞紐地位,這進一步驗證了NDUFB9高表達與腫瘤細胞干性之間的潛在關聯及其在促進腫瘤發生發展中的潛在作用。此外,EHD2在上述增殖活躍的腫瘤細胞亞群中表達量極低,與之前的結果一致,即EHD2與腫瘤細胞干性之間呈負相關關聯。

圖8 骨肉瘤細胞亞型識別、細胞通信與細胞周期分析

基于單細胞水平進一步分析先前識別的57個差異表達干性相關基因(R>0.6或R<-0.5)是否在不同細胞亞群間仍存在差異表達。結果顯示,與腫瘤干性呈負相關關系的20個基因中,共有18個基因表現出明顯差異表達的特征,并且這些基因在淋巴細胞、MSCs和成骨細胞等增殖活躍的細胞亞型中表達水平極低,見圖9。此外,與腫瘤干性呈正相關關系的37個基因中,共有25個基因表現出明顯差異表達的特征,并且這些基因在淋巴細胞、MSCs和成骨細胞等增殖活躍的細胞亞型中普遍高表達,且表現出一定的共表達趨勢,這進一步驗證了這些基因與OS腫瘤干性之間的潛在關聯,見圖10。

圖9 單細胞測序識別的差異表達干性相關基因(與mRNAsi呈負相關關系)的細胞特征圖

圖10 單細胞測序識別的差異表達干性相關基因(與mRNAsi呈正相關關系)的細胞特征圖

3 討 論

CSCs是參與腫瘤侵襲轉移及耐藥的關鍵細胞亞群[24]。腫瘤細胞分化程度越低,則其干細胞樣特征越顯著,腫瘤惡性程度也越高。腫瘤干性的本質是多基因表達異常介導下游信號通路失調,通過維持CSCs的自我更新或誘導腫瘤細胞去分化而獲得干細胞樣特征。然而,目前OS干性特征的發生機制尚未明確。本研究通過TARGET數據庫聯合機器學習描繪了OS干性特征,并挖掘了調控OS腫瘤干性特征的潛在機制。結果發現,mRNAsi與OS患者的原發部位、生存狀態和腫瘤轉移相關,是影響OS患者預后的不利因素。

基于OS患者RNA-seq數據,本研究共篩選出658個關鍵基因。這些基因不僅與腫瘤干性顯著相關,而且與OS患者預后及生存顯著相關,是OS干細胞靶向治療的潛在靶點。在關鍵基因中,NDUFB9和EHD2與mRNAsi之間的相關性最強,且均與OS患者的總生存期顯著相關。

NDUFB9編碼線粒體膜呼吸鏈NADH脫氫酶(復合物Ⅰ)的一個輔助亞單位[25]。NDUFB9突變可導致復合物Ⅰ缺乏[26]。復合物Ⅰ可以氧化NADH,減少泛醌,并通過線粒體內膜運輸質子,從而產生質子動力。眾所周知,線粒體和能量代謝的改變與癌變之間存在著緊密關聯。癌細胞會偏向以無氧糖酵解取代正常的有氧循環進行能量供給,所以癌細胞線粒體的生物氧化方式與正常細胞有顯著區別,這一現象被稱為瓦氏效應。瓦氏效應促進了腫瘤細胞適應低氧和快速增殖[27]。多項研究證明上調NDUFB9可促進腫瘤轉移[28-29]。最近有研究基于單細胞測序分析和轉錄組分析發現NDUFB9在子宮內膜癌組織中顯著上調,在子宮內膜癌的發生發展中發揮著重要作用[30]。此外,NDUFB9表達異常引起的復合物Ⅰ活性改變增強了癌細胞的侵襲性,導致患者的不良預后[29]。

EHD2編碼Eps15同源結構域,是一種腫瘤抑制基因,在多種腫瘤中表達水平明顯降低[31]。EHD2表達與組織學分級、淋巴結轉移和腫瘤大小有關。此外,EHD2的過度表達抑制癌細胞的遷移和侵襲,而敲除EHD2則促進癌細胞的遷移和侵襲[31]。最新研究表明,EHD2與多種癌癥顯著相關。與正常組織相比,食管鱗狀細胞癌中EHD2蛋白水平明顯降低,并且EHD2的低表達促進了食管鱗狀細胞癌細胞的遷移[32]。在卵巢癌、惡性黑色素瘤和肝細胞癌中也發現了EHD2的明顯低表達[33-35],而低表達EHD2的患者生存期更短,預后更差。

腫瘤異質性包括瘤內異質性和瘤間異質性,源于形態學和表觀遺傳學的可塑性,也可歸因于基因表達譜不同的腫瘤細胞克隆共存[36]。單細胞測序在探究腫瘤異質性方面具有獨特優勢,特別是在OS這種具有高度異質性、難以用整體測序方法進行描述的復雜腫瘤中。

本研究繪制了OS患者的轉錄組學單細胞圖譜,并探究了復雜的OS腫瘤微環境以及腫瘤異質性。通過對11個OS患者的scRNA-seq數據進行質控、降維后,識別并注釋了10種不同的細胞亞群,提示OS腫瘤細胞之間存在異質性。在通過質控的123 307個單細胞中,成骨細胞比例最高,其次是髓系細胞和成纖維細胞。在復雜的OS腫瘤微環境中,細胞之間存在著廣泛的細胞通信。本研究使用“iTALK”R軟件包來分析不同細胞亞群間的受體配體對,從而預測細胞間通訊模式。成軟骨細胞和MSCs占據了通訊網絡中心,主要通過配體COL1A1和COL1A2發送信號,并通過受體LRP1接收信號,在該通訊網絡中發揮重要作用。COL1A1與COL1A2編碼產物是Ⅰ型膠原蛋白的重要組成成分,可以形成Ⅰ型膠原蛋白的三螺旋結構[37]。Ⅰ型膠原蛋白是骨骼中膠原蛋白形成所必需的,并且它也參與細胞外基質(extra cellular matrix,ECM)的合成和細胞形態改變。研究表明,COL1A1和COL1A2在多種腫瘤組織中的表達比在正常組織中更高,并且都與腫瘤細胞的侵襲和增殖有關[38-39]。重要的是,本研究鑒定的關鍵基因NDUFB9在多種細胞亞群中表達水平都較高,包括間充質干細胞和成骨細胞等;而EHD2在各腫瘤細胞亞群的表達量極低或不表達。NDUFB9與多種特異性干性標志物存在共表達趨勢,進一步明確了NDUFB9與OS腫瘤細胞干性的密切聯系。相反,EHD2在各類OS細胞亞群中低表達或幾乎不表達。結合細胞周期分析結果發現,分裂增殖旺盛的細胞亞群(如間充質干細胞和成骨細胞)高表達NDUFB9,低表達EHD2?；趩渭毎麥y序,可更加深入地認識OS腫瘤細胞的分子特征,為剖析其腫瘤微環境、探索腫瘤異質性、尋找診斷和治療靶點提供理論依據。

總之,本研究通過整合生物信息學分析篩選了OS預后相關的關鍵干性基因,為靶向藥物研究提供了重要的分子基礎。此外,基于OS患者scRNA-seq數據,本研究描繪了OS單細胞轉錄組圖譜,初步揭示了OS腫瘤異質性,分析了潛在的細胞通信模式,從單細胞水平驗證了關鍵基因在OS腫瘤細胞中的表達水平。

致謝：特別感謝上海格致中學(奉賢校區)的湯慧婷同學在論文參考文獻整理中的貢獻。