劉懿天，莊 然，丁 勇

(1空軍軍醫(yī)大學唐都醫(yī)院骨科，陜西 西安 710038；2空軍軍醫(yī)大學基礎醫(yī)學院免疫學教研室，陜西 西安 710032)

骨肉瘤(osteosarcoma，OS)是一種進展程度快、易發(fā)生早期轉移的惡性骨腫瘤，好發(fā)于青春期以及65歲以上的老年人[1]。OS最常見的病變部位為股骨遠端及脛骨近端的干骺端[2]，疼痛和腫脹是其早期最常見的癥狀。由于好發(fā)于青春期人群，癥狀常被誤認為是生長痛或劇烈運動導致的疼痛而耽誤早期診斷。作為一種惡性程度極高的骨腫瘤，OS具有親肺轉移的特性，在發(fā)生轉移的患者中，有高達82%的患者存在肺轉移[3]，并且大約20%的患者在確診OS時便已發(fā)生肺轉移[1]。盡管自上世紀80年代以來通過手術切除、化療等各種治療方式使未發(fā)生肺轉移的局灶性患者長期存活率已經(jīng)保持在60%以上，但對于已經(jīng)發(fā)生肺轉移的患者來說，其長期存活率卻一直沒有改善，仍然低于40%[4]。因此，亟需了解OS發(fā)生發(fā)展的分子機制，并尋找新的生物標志物用于OS的早期診斷和預后評估。

雖然OS具有高度侵襲性，但其發(fā)病率卻非常低，在0～19歲的人群中每年每百萬人只有不到6人確診[5]。這導致OS難以像結直腸癌、乳腺癌那樣擁有完整的大樣本量的轉錄組圖譜，癌癥基因組圖譜(The Cancer Genome Atlas，TCGA)數(shù)據(jù)庫中31種癌癥分類中甚至沒有OS，只有肉瘤。盡管既往已有研究采用生物信息學的方法對OS的預后標志物進行了初步篩選，但不是樣本量過小就是數(shù)據(jù)來源單一，只在特定數(shù)據(jù)集中有效[6]。

本研究從多個不同數(shù)據(jù)庫選取與OS相關的多種不同類型的基因表達譜數(shù)據(jù)集，分析了OS患者腫瘤組織與正常骨組織的差異表達基因(differentially expressed genes，DEGs)，通過功能富集分析、腫瘤微環(huán)境(tumor microenvironment，TME)分析以及單因素和多因素Cox回歸分析篩選影響患者預后的相關基因，構建預后預測模型并進行了蛋白水平的初步驗證，為OS患者的診斷和預后評估提供臨床新思路。

1 材料與方法

1.1 材料

1.1.1 生物信息學數(shù)據(jù) 從TARGET數(shù)據(jù)庫(https://ocg.cancer.gov/programs/target)檢索下載TARGET-OS項目的RNA-seq表達矩陣，其表達量單位為TPM。從GEO數(shù)據(jù)庫(http://www.ncbi. nlm.nih.gov/geo)中檢索下載OS相關的RNA-seq或Microarray表達矩陣及對應注釋信息，其單位為FPKM。從GEO數(shù)據(jù)庫獲取單細胞RNA測序(single cell RNA-seq，scRNA-seq)數(shù)據(jù)，表達量單位為Count。從CCLE數(shù)據(jù)庫(https://sites.broadinstitute.org/ccle)檢索下載人OS細胞系的表達矩陣，表達量單位為TPM。31種癌癥分類目標基因的TPM表達譜文件獲得自GEPIA網(wǎng)站(http://gepia.cancer-pku.cn/)[7]的TCGA數(shù)據(jù)庫。部分需要進行數(shù)據(jù)處理的數(shù)據(jù)集及相關信息見表1。

表1 數(shù)據(jù)集基本信息

1.1.2 主要試劑和儀器 細胞培養(yǎng)箱(Thermo公司，美國)；細胞系FOB1.19、MG63、143B、HOS(ATCC中心，美國)、K7M2(普諾賽生命科技有限公司，中國)；高糖DMEM培養(yǎng)基、MEM培養(yǎng)基、胎牛血清(Gibco公司，美國)；青霉素-鏈霉素、胰蛋白酶、40 mL/L多聚甲醛溶液(米鼠生物科技有限公司，中國)；RIPA裂解液(新賽美生物科技有限公司，中國)；聚偏二氟乙烯膜(Kodak公司，美國)；50 mL/L牛血清白蛋白溶液、HE染色套裝(塞維爾生物科技有限公司，中國)；兔抗NEDD4L(博士德生物工程有限公司，中國)；兔抗β-actin、鼠抗兔二抗(三鷹生物技術有限公司，中國)；化學發(fā)光成像系統(tǒng)(Bio-Rad公司，美國)；IX73 型倒置相差顯微鏡(Olympus公司，日本)；BALB/c小鼠(空軍軍醫(yī)大學實驗動物中心，中國)。

1.2 方法

1.2.1 原始數(shù)據(jù)處理與差異表達基因篩選 使用R語言dplyr包將從GEO數(shù)據(jù)庫下載的FPKM表達矩陣轉換為TPM表達矩陣，并進行標準化處理[log2(TPM+1)]。使用sva包的Combat函數(shù)進行批次矯正并整合。使用stats包進行主成分分析(principal components analysis，PCA)并用ggplot2包對分析結果進行可視化展示。使用limma包進行DEGs分析，以經(jīng)本亞明-霍赫貝格法多重檢驗校正后的P值和差異倍數(shù)(fold change，F(xiàn)C)作為判定DEG的標準，|log2(FC)|>2的基因在本文被定義為DEGs，P<0.05表示差異有統(tǒng)計學意義。

1.2.2 功能富集分析 使用GSEA軟件進行基因集富集分析(gene set enrichment analysis，GSEA)[8]，以目標基因表達量排序，中位數(shù)為閾值，將患者分為高表達組和低表達組，使用從MSigDB(http://www.gsea-msigdb.org/gsea/downloads.jsp)下載含有京都基因和基因組百科全書[9]通路信息的c2.cp.kegg.v7.4.symbols. gmt參考子集合，用于評估高、低表達組基因集的變化，P<0.05和FDR<0.25的通路在本文被定義為顯著富集的通路。

1.2.3 構建單因素和多因素預后模型 使用survival包和survminer包對上述方法1.2.1合并后數(shù)據(jù)集的患者年齡、性別以及所有DEGs進行單因素Cox回歸分析，篩選P<0.05的變量進行多因素Cox回歸分析。使用rms包[10]繪制帶有多因素Cox回歸分析結果P值的列線圖。使用Predict函數(shù)對多因素進行風險評分，對評分進行排序并繪制相應基因的熱圖。

1.2.4 免疫浸潤分析 使用ESTIMATE算法[11]對數(shù)據(jù)集進行分析，并獲得每個樣本的Stromal評分、Immune評分和腫瘤純度(tumor purity，TP)。按照TP排序，以中位值為界，將樣本分為高純度組和低純度組，使用Wilcoxon秩和檢驗對兩組進行統(tǒng)計學檢驗并用ggplot2包繪制小提琴圖。使用ggpubr包計算TP與目標基因的Spearman相關性系數(shù)并繪制相關性散點圖。P<0.05代表TP與目標基因的表達量存在顯著相關性。

1.2.5 單細胞測序分析 使用seurat包[12]對數(shù)據(jù)集進行數(shù)據(jù)標準化處理。使用統(tǒng)一流形逼近與投影(uniform manifold approximation and projection，UMAP)算法[13]對數(shù)據(jù)進行降維分析，使用Louvain算法對數(shù)據(jù)進行聚類。使用UMAPPlot和FeaturePlot函數(shù)分別繪制UMAP降維聚類散點圖和目標基因分布特征圖。

1.2.6 預后模型的評估 使用生存曲線和受試者工作特征(receiver operating characteristic，ROC)曲線下面積(area under the curve，AUC)評估預后模型的效能。使用survival包和survminer包繪制總生存曲線，并采用Log rank檢驗評估不同分組之間的差異顯著性。使用pROC包繪制單因素或多因素在1、3、5年3個時間點的ROC曲線并計算相應AUC。AUC>0.7代表模型可以有效預測預后。

1.2.7 細胞培養(yǎng) 將凍存細胞置于37 ℃恒溫水浴鍋中復溫至細胞液完全解凍。吸取復溫后細胞液至15 mL離心管，加入10 mL含100 mL/L胎牛血清的完全培養(yǎng)基(FOB1.19、143B和K7M2使用高糖DMEM培養(yǎng)基；MG63和HOS使用MEM培養(yǎng)基)稀釋，1 000 r/min離心5 min，棄上清，加入對應完全培養(yǎng)基重懸，接種于6 cm培養(yǎng)皿，置于37 ℃、50 mL/L CO2培養(yǎng)箱中，5 h后換一次液。待細胞生長至80%以上融合時，消化并擴增至10 cm培養(yǎng)皿，以后常規(guī)傳代培養(yǎng)。取第3～5代對數(shù)生長期細胞進行后續(xù)實驗。

1.2.8 蛋白質(zhì)印跡實驗 將處于對數(shù)生長期的3種OS細胞系(MG63、HOS和143B)以及人成骨細胞FOB1.19使用PBS洗2次后提取蛋白，BCA法蛋白定量，加入loading buffer后，100 ℃水浴15 min，4℃ 15 000g離心5 min。利用100 g/L SDS-PAGE凝膠進行蛋白電泳。蛋白樣品按10 μg/孔上樣電泳，恒流濕轉法轉膜。轉膜結束后，50 mL/L牛血清白蛋白溶液室溫封閉30 min，孵育一抗：兔抗NEDD4L、β-actin，一抗?jié)舛染鶠?∶1 000。4 ℃搖床過夜，之后加入山羊抗兔二抗(1∶5 000)室溫孵育1 h，最后進行化學發(fā)光顯影并拍照。使用Image Lab軟件進行灰度分析。

1.2.9 小鼠OS肺轉移模型的建立 將處于對數(shù)生長期的小鼠高轉移OS細胞系K7M2消化、離心后用PBS重懸，調(diào)整濃度至5×109個/L，按200 μL/只的量尾靜脈注射至4周齡雄性BALB/c小鼠，21 d后肺組織取材，40 mL/L多聚甲醛溶液室溫固定，石蠟包埋。

1.2.10 切片染色 將小鼠肺轉移灶石蠟標本按厚度4 μm切片，脫蠟后進行HE和組化染色。對于HE染色，依次使用蘇木精和伊紅染色，脫水封片，可使細胞核呈藍色，細胞質(zhì)呈紅色。對于組化染色，抗原修復后孵育一抗：兔抗NEDD4L(1∶50)，4 ℃過夜。之后用50 mL/L牛血清白蛋白溶液封閉2.5 h，并用二氨基聯(lián)苯胺顯色，得到棕色產(chǎn)物。使用IX73型倒置相差顯微鏡進行拍照記錄。

1.2.11 統(tǒng)計學分析 生物信息學分析方法各部分已述。蛋白質(zhì)印記實驗獨立重復3次，結果使用GraphPad Prism 9.0軟件作圖并進行非配對t檢驗分析。P<0.05表示差異有統(tǒng)計學意義。

2 結果

2.1 DEG分析與單因素Cox回歸分析

使用stat包對數(shù)據(jù)集GSE126209(OS患者12例，健康骨組織11例)進行PCA分析并繪制降維散點圖(圖1A)。圖中點之間的距離越近代表點對應的樣本組成越相似，在整體上OS患者與健康骨組織存在顯著差異，OS患者組內(nèi)部也存在較為明顯的異質(zhì)性。

使用limma包對數(shù)據(jù)集GSE126209進行DEG分析，共得到2 806個DEGs，其中上調(diào)基因2 284個，下調(diào)基因522個，使用ggplot2包繪制DEGs火山圖(圖1B)。圖中橫坐標代表log2處理后的FC，F(xiàn)C越大，橫坐標離0越遠；縱坐標代表-log10矯正后的P值，差異越顯著，縱坐標越大；藍點代表在OS患者中表達量下調(diào)的基因，紅點代表上調(diào)基因，灰點代表差異不顯著的基因。

使用survival包對GSE16091、GSE21257和GSE39055合并后的數(shù)據(jù)集(OS患者123例)患者年齡、性別以及所有DEGs進行單因素Cox回歸分析并繪制森林圖(圖1C)。篩選到9個基因P<0.05，其中TFPI2和CREB3L1表達水平下調(diào)，其余為上調(diào)基因，年齡和性別不是患者生存預后相關的獨立因素。

A：OS患者與健康骨組織PCA降維散點圖；B：OS患者與健康骨組織DEG火山圖；C：DEGs與基礎信息單因素Cox回歸分析結果森林圖。OS:骨肉瘤；DEGs：差異表達基因；FC：差異倍數(shù)。

2.2 多因素Cox回歸預后模型的構建與評估

對2.1得到的9個基因進行多因素Cox回歸分析以構建預后模型，評估這9個基因在123個樣本中的預后顯著性。采用rms包繪制所有多因素分析P<0.05的基因的列線圖(圖2A)。結果顯示，只有神經(jīng)前體細胞表達發(fā)育下調(diào)樣基因4(neural precursor cell-expressed developmentally down-regulated 4 like，NEDD4L)、TFBI2和WASF3共同對預后存在顯著影響(P<0.05)，其中NEDD4L對總分貢獻度最大，P值最小(P<0.000 1)。以這3個基因構成的回歸模型一致性指數(shù)為0.76，提示模型與實際結果較為一致。

使用Predict函數(shù)計算多因素風險評分，對評分進行排序并繪制這3個基因的風險評分熱圖(圖2B)。圖中上方柱狀圖縱坐標為患者的風險評分，橫坐標按照風險評分由低到高的順序從左向右排列，中、下兩圖橫坐標均與之對應；中間散點圖縱坐標為患者生存時間；下方熱圖為3個基因Z-score歸一化后的表達量。結果顯示，NEDD4L表達量越低、WASF3和TFPI2表達量越高，患者風險評分越高、預后越差。

按風險評分排序，以中位數(shù)(-0.045 4)為閾值，將患者分為高風險組(62例)和低風險組(61例)，繪制總生存曲線(圖2C)與ROC曲線(圖2D)。結果顯示，低風險組與高風險組在總生存率上存在顯著差異(P<0.000 1)，并且低風險組預后較好，其1、3和5年AUC分別為0.76、0.84、0.86。圖2C～D說明以NEDD4L為主的多因素風險評分可以有效預測OS患者的預后。

2.3 NEDD4L TME分析

使用estimate包對數(shù)據(jù)集TARGET-OS進行免疫浸潤分析，獲得88例OS患者的TP。按照TP排序，將88個樣本分為高純度組(44例)和低純度組(44例)，對NEDD4L的表達量進行Wilcoxon秩和檢驗，并繪制小提琴圖(圖3A)。結果顯示，NEDD4L在高純度組表達量顯著高于低純度組(P=0.005 7)。

使用Spearman相關性系數(shù)對TP與NEDD4L進行相關性分析并繪制相關性散點圖(圖3B)。結果顯示，TP與NEDD4L表達量顯著正相關(R=0.35，P=0.001 0)。圖3A～B共同說明NEDD4L在TME中主要表達于OS細胞。

通過對CCLE數(shù)據(jù)庫和GEPIA網(wǎng)站進行檢索獲得前文2.1部分篩選到的9個基因在不同OS細胞系以及TCGA腫瘤分類上的表達量并繪制熱圖(圖3C)。結果顯示，NEDD4L相較于其他8個基因高表達于多種OS細胞系以及多種其他腫瘤。

使用seurat包對數(shù)據(jù)集GSE162454進行Louvain聚類后最終獲得了5群細胞，通過UMAP算法對數(shù)據(jù)集進行降維分析并繪制NEDD4L的特征圖(圖3D)。結果顯示，NEDD4L在TME中主要分布在髓樣細胞和OS兩群細胞。

A：免疫浸潤分析TME中NEDD4L表達量小提琴圖；B：TME中NEDD4L與TP相關性散點圖；C：9個DEGs在OS患者、細胞系以及TCGA腫瘤分類的表達量熱圖；D：OS單細胞測序UMAP降維散點圖(上)以及NEDD4L表達特征圖(下)。OS:骨肉瘤；TME：腫瘤微環(huán)境；TP：腫瘤純度；DEGs：差異表達基因。

2.4 基于NEDD4L的單因素預測模型及驗證

對數(shù)據(jù)集GSE14359(原發(fā)灶10例，轉移灶8例)以上9個基因的表達量進行Wilcoxon秩和檢驗并繪制分邊小提琴圖(圖4A)，結果顯示NEDD4L在肺轉移灶的表達量顯著高于原發(fā)灶(P=0.005 7)。

對TARGET-OS進行GSEA富集分析，以NEDD4L表達量的中位數(shù)為閾值，將88位患者分為高表達組(44例)和低表達組(44例)，檢驗兩組之間基因集的差異并用ggplot2包繪制P<0.05的通路(圖4B)。圖4B分為上、中、下三部分，三部分的橫坐標相同：上面部分縱坐標為富集分數(shù)(enrichment score，ES)，每個基因都有ES，從左至右ES依次相加并連成折線，折線峰值處即為基因集的ES；中間部分每條豎線對應相應通路下的一個基因；下面部分為每個基因的Rank值分布圖，縱坐標Rank值即為基因log標準化后的FC；右側圖例為通路名稱及對應的ES及P值。結果顯示，低表達NEDD4L與包括子宮內(nèi)膜癌、甲狀腺癌、前列腺癌、結腸癌在內(nèi)的多種癌癥通路有關。

繪制NEDD4L的總生存曲線并用Log rank檢驗進行評估(圖4C)。結果顯示，高表達組生存預后顯著好于低表達組(P=0.002 3)。繪制NEDD4L的ROC曲線(圖4D)，其1、3和5年的AUC分別為0.78、0.72和0.72。圖4C～D說明NEDD4L可以有效預測OS患者的預后。

2.5 NEDD4L表達情況的初步驗證

利用蛋白質(zhì)印跡實驗檢測人OS細胞系MG63、143B和HOS以及人成骨細胞系FOB1.19 NEDD4L蛋白表達情況(圖5A)，結果顯示，NEDD4L在3種OS細胞系中的表達量均顯著低于FOB1.19細胞(MG63：P=0.029 4；143B：P=0.007 1；HOS：P=0.026 0；圖5B)。

利用尾靜脈注射的方法構建小鼠OS肺轉移模型，21 d后鼠肺取材固定并進行HE(圖6A)和組化染色(圖6B)，結果顯示，NEDD4L主要定位在肺組織，OS轉移灶未見明顯染色。蛋白質(zhì)印跡和組化結果共同說明NEDD4L在OS中的表達水平低于健康組織。

A：不同人OS細胞系NEDD4L表達情況蛋白質(zhì)印跡條帶結果(n=3)；B：不同人OS細胞系NEDD4L表達情況統(tǒng)計圖。 aP<0.05， bP<0.01。

A：小鼠OS肺轉移灶HE染色結果；B：小鼠OS肺轉移灶NEDD4L組化染色結果。比例尺為100 μm，黑色箭頭為OS肺轉移灶。

3 討論

本研究通過GEO數(shù)據(jù)庫、TARGET數(shù)據(jù)庫、CCLE數(shù)據(jù)庫以及TCGA數(shù)據(jù)庫獲得了多個不同測序類型、不同測序?qū)ο蟮腛S相關數(shù)據(jù)集。對RNA-seq數(shù)據(jù)集GSE126209進行差異表達分析，共獲得2 806個DEGs，其中上調(diào)基因2 284個，下調(diào)基因522個。在Microarray數(shù)據(jù)集中對DEGs進行單因素Cox回歸分析，共得到9個影響預后的基因，對這9個基因進行多因素Cox回歸分析并構建預后預測模型，只有NEDD4L、WASF3和TFPI2與預后顯著相關，其中NEDD4L對預后的貢獻最大。進一步對NEDD4L進行免疫浸潤分析、單細胞測序分析發(fā)現(xiàn)，在TME中NEDD4L主要表達于OS細胞，且轉移灶的表達量顯著高于原發(fā)灶。使用ROC曲線對NEDD4L的單因素及多因素Cox回歸分析預后預測模型進行評估，無論是1年、3年還是5年，其單因素和多因素預后預測模型AUC均大于0.7。這表明本研究所構建的模型在轉錄組水平可能對OS患者的預后進行有效預測。最后，本研究初步對3種人OS細胞系和人成骨細胞FOB1.19進行了蛋白水平的驗證，同時對小鼠OS肺轉移灶進行了組化染色，結果均提示與正常組織對比，NEDD4L在OS中的表達減少，與預后預測變化趨勢一致。

NEDD4L，又叫NEDD4-2，是一種E3泛素連接酶，能夠與多種膜蛋白結合并調(diào)節(jié)其功能，對于維持細胞內(nèi)環(huán)境穩(wěn)態(tài)是必不可少的[14]。NEDD4L在多種癌癥的診斷與預后中已經(jīng)被證明存在保護作用，如蔣興旺等[15]利用實時熒光PCR、蛋白質(zhì)印跡等方法對胃癌組織與癌旁組織進行NEDD4L表達量的檢測，發(fā)現(xiàn)低表達NEDD4L可能與胃癌的發(fā)生密切相關，且可能作為判斷胃癌預后的潛在指標；LI等[16]利用生物信息學方法證明低表達NEDD4L促進肺腺癌的進展并進行了實驗驗證；GUARNIERI等[17]利用miRNA直接抑制NEDD4L，發(fā)現(xiàn)這種抑制可以促進乳腺癌的發(fā)生。在本研究中，蛋白質(zhì)印跡實驗表明、NEDD4L在OS細胞系中的表達量顯著低于成骨細胞對照，而GSEA富集分析結果也顯示低表達NEDD4L與多種腫瘤發(fā)生通路相關，這提示NEDD4L的作用機制在不同腫瘤之間存在共性，其保護作用可能極為重要。有趣的是，免疫浸潤分析和單細胞測序分析結果均表明在TME中NEDD4L主要表達在OS細胞上，其在免疫細胞或間質(zhì)細胞上的表達量反而更低。盡管NEDD4L在其他腫瘤中的作用已經(jīng)初步得到研究，但尚無其在OS中相關作用的直接報道。既往研究表明NEDD4L可以下調(diào)腫瘤的TGF-β和Wnt信號通路[14，18]，而這兩個通路在OS的進展中發(fā)揮著重要作用[19-20]，因此高表達NEDD4L很有可能對OS患者的預后具有保護作用，本研究通過生物信息學方法和初步實驗驗證推導出了這一點。

本研究具有以下優(yōu)勢：第一，在不同數(shù)據(jù)庫來源的傳統(tǒng)RNA-seq測序、Microarray以及單細胞測序數(shù)據(jù)集上使用多種生物信息學方法初步論證了NEDD4L對OS的預后具有預測價值；第二，首次構建以NEDD4L為主的OS預后預測模型，揭示NEDD4L對于OS可能具有保護作用。同時本研究存也存在局限性：第一，本研究主要結論由生物信息學方法分析得出，僅對NEDD4L的表達情況進行了初步驗證，缺少大樣本量的預后驗證；第二，部分臨床信息并非公開，無法納入分析，而諸如病理類型、分期、化療、手術等臨床自變量與OS預后密切相關，NEDD4L與這些臨床自變量之間可能存在不同程度的相關性。因此，在未來的研究中我們將進一步收集OS患者的臨床信息和標本，通過體內(nèi)以及體外實驗對NEDD4L的作用進行進一步驗證并對其在TME中表達情況相反的原因進行探討。

綜上所述，本研究基于生物信息學方法篩選目標基因，并對其進行大樣本分析，對預測OS患者的預后評估提供臨床新思路。