基于生物信息學篩選三陰性乳腺癌潛在預后相關基因

鄭鈞元,王奕然,王 慧,秦一帆,田玉靜,崔榮軍,趙福陽,楊旭芳

(牡丹江醫學院 1.病理生理學教研室;2.第一臨床醫學院;3.生物化學與分子生物學教研室;4.公共衛生學院,黑龍江 牡丹江 157011)

乳腺癌(breast cancer,BC)是最常見的女性惡性腫瘤,每年有40多萬女性患者死于乳腺癌,其發病率常年居高不下[1]。其各分子類型間的預后,惡性程度,發病年齡,臨床特點等均有不同。乳腺癌在分子層面上諸多不同的特點已經越來越受到人們的關注,其分子分型也已被越來越多的臨床醫生們作為進行乳腺癌治療的主要參考標準,在臨床工作中,絕大多數學者根據孕激素受體(progesterone receptor,PR)、雌激素受體(estrogen receptor,ER)、人類表皮生長因子受體2(human epidermal growth factor receptor-2,HER-2)以及Ki-67的結果,將乳腺癌劃分為4種亞型,包括Luminal A型、Luminal B型、HER2過表達型和三陰性乳腺癌(triple-negative breast cancer,TNBC)[2]。TNBC是指雌激素受體、孕激素受體和人表皮生長因子受體2的表達均為陰性的一類乳腺癌,其對內分泌治療與曲妥珠單抗靶向治療均不敏感。TNBC具有臨床分期晚、侵襲度高、易發生轉移以及組織學分級高等特點[3]。因為缺乏有效的治療靶點,所以除手術治療外傳統化療仍然是TNBC患者術后重要治療手段。因此,通過尋找TNBC預后相關基因,進而找到更為可靠精確的靶向生物標志物必將成為目前研究的熱點。

生物信息學發源于上世紀70年代,隨著近年來的分子生物學研究的深入和計算機技術的進步,生物信息學迅猛發展,各類組學概念也應運而生。隨著步入21世紀,微陣列技術的興起以及高通量測序的發展,在各類檢測中具有更高分辨率和更廣的檢測范圍。檢測技術日益進步所帶來成本下降,使生物信息學得到更多應用,從而推動了各類疾病診療技術的進步[4]?，F常用于腫瘤研究的數據庫有:GEO數據庫、TCGA數據庫等,利用這些數據庫中上傳的樣本測序信息,能夠從多種組學的角度闡明乳腺癌的發生發展機制,明確乳腺癌預后、治療相關的分子標志物,為乳腺癌的臨床診斷、臨床個性化治療和對患者預后的綜合評估提供新的思路[5]。

隨著臨床醫學的進步以及分子生物學的發展,醫學研究已逐步進入分子時代,不同分型的乳腺癌類型間,其發展進程、治療方式及癌組織、癌細胞對治療的反應以及預后情況都不盡相同,對于惡性程度高的TNBC尋找更為可靠的預后相關基因,并發掘潛在的生物學標志物,從而將免疫分型和乳腺癌臨床病理學分類更好更全面的綜合起來,臨床醫生才能以更科學的視角為乳腺癌病人量身定制有效的個性化治療方案[6]。

截至2022年,在PubMed上以triple-negative breast cancer,bioinformatics和prognosic genes三個關鍵詞檢索,共能檢索到自2008年至2022年共有127篇文獻,大多數是從某種特定生物過程或多組學角度進行研究的,但目前少有學者就某一關鍵基因進行深入研究,因此,我們就這一思路,通過生物信息學手段進行分析,對三陰性乳腺癌和非三陰性乳腺癌樣本數據進行分組,篩選出三陰性乳腺癌預后相關基因,并聯合KEGG、GO途徑富集分析,免疫浸潤分析,ScRNA-seq等結果探討導致TNBC預后差,難治療的潛在因素。

1 材料與方法

1.1 數據收集與處理為了獲取三陰型乳腺癌組織中基因表達量和患者臨床信息,首先使用GEO數據庫(Gene Expression Omnibus,GEO)(https://www.ncbi.nlm.nih.gov/geo/),下載GSE25066(Hatzis C等人,2011年)、GSE62931(Boac BM等人,2019年)、GSE129559(Sinn BV等人,2019年)三個TNBC相關的基因表達微陣列。GSE25066數據集中包括27個三陰性乳腺癌組織樣本和39個非三陰性乳腺癌組織樣本。GSE62931數據集中包含47個三陰性乳腺癌組織樣本和53個非三陰性乳腺癌組織樣本。GSE129559數據集中包含17個三陰性乳腺癌組織樣本和99個非三陰乳腺癌組織樣本。由于數據來自在線數據庫,因此不需要倫理委員會的進一步批準。

通過GEO2R標準化基因微陣列(http://www.ncbi.nlm.nih.gov/geo/geo2r)篩選癌組織與非癌組織的差異表達基因,定義滿足|LogFC|>1和FDR<0.05的基因為差異表達基因(Differentially Expressed Genes,DEGs),移除沒有相應基因符號的探針集,在線篩選出TNBC相關的DEGs,并繪制火山圖及熱圖。

為了使篩選結果更加精準,在癌癥基因組數據庫(The Cancer Genome Altas,TCGA)(https://www.cancer.gov)中下載了TNBC基因表達譜,使用R語言篩選差異表達基因。利用Draw Venn Diagram在線網站(http://bioinformatics.psb.ugent.be/webtools/Venn/)將上述四個數據集的差異表達基因取交集,最終得到共表達差異基因。

1.2 差異基因的KEGG與GO富集分析為了進一步了解TNBC差異基因所富集的相關通路以及功能,我們通過DAVID在線數據庫(DAVID;https://david.ncifcrf.gov)進行GO和KEGG富集分析。GO富集分析包括三個基本內容:分子功能(Molecular Function,MF)、細胞組分(Cellular Component,CC)、生物過程(Biological Process,BP)。P<0.05表明該差異具有統計學意義。

1.3 TNBC預后相關差異表達基因的篩選與確定

1.3.1 篩選預后相關差異表達基因 TNBC患者的臨床信息數據從TCGA數據庫下載。為了得到預后相關的DEGs,我們使用單變量Cox回歸分析篩選預后相關基因。P值<0.05表明該差異具有統計學意義。并利用Draw Venn Diagram在線網站繪制Venn圖篩選預后相關差異表達基因。

1.3.2 生存曲線驗證預后差異基因 根據核心基因表達的中位值將TNBC患者分為高低兩組,使用R軟件繪制生存曲線驗證預后差異基因,P<0.05表明該差異具有統計學意義。

1.4 核心基因的ROC曲線分析去除只在單因素Cox回歸分析中或只在生存曲線中有統計學意義的基因,根據TCGA數據庫下載的TNBC基因表達譜中上述預后相關基因的表達量,分別進行ROC曲線分析,以進一步評估這些基因對TNBC診斷的特異性與敏感性。

1.5 免疫浸潤分析用Cibersort算法評估預后相關基因的表達與TNBC中的腫瘤浸潤性免疫細胞相關性。

1.6 單細胞RNA測序使用人類蛋白圖譜數據庫(the human protein atlas,HPA)(https://www.proteinatlas.org/)分析COL9A3在正常乳腺組織中的分布。

2 實驗結果

2.1 篩選三陰性乳腺癌中差異表達基因通過GEO數據庫篩選差異表達基因,火山圖和熱圖結果顯示,GSE25066數據集篩選出672個差異表達基因,其中包括331個上調基因和341個下調基因。GSE62931數據集篩選出2 318個差異表達基因,其中包括1 248個上調基因和1 070個下調基因。GSE129559數據集篩選出994個差異表達基因,其中包括456個上調基因和538個下調基因(圖1A～B)。將上述三個數據集獲得的DEGs與TCGA數據庫中的DEGs取交集得到了217個共表達DEGs(圖1C)。

圖1 三陰性乳腺癌基因差異表達基因的鑒定

2.2 KEGG和GO富集分析通過KEGG和GO富集分析了解TNBC差異基因的生物學功能, KEGG富集分析表明DEGs主要富集于雌激素信號通路(圖2A),GO結果顯示,生物過程富集在乳腺上皮發育,腺體發育,腺體形態發生等方面。細胞組分主要涉及微絨毛膜、細胞基底外側膜等。分子功能主要富集于跨膜受體蛋白酪氨酸激酶等(圖2B～D)。

圖2 217個差異表達基因的富集分析

2.3 預后生存相關核心基因的篩選及驗證

2.3.1 預后相關差異表達基因的篩選 通過單因素Cox回歸分析篩選出了971個預后相關基因,將這971個預后相關基因與217個篩選出的DEGs取交集,共得到了11個預后相關DEGs,其中的ALCAM、KITLG、MAST4、MATN3、MCCC2、UGCG的HR>1,提示其可能為風險因素,而TJP3、TTYH1、SOX10、IL12RB2、COL9A3的HR<1提示其可能為保護性因素(P<0.05)(圖3)。

圖3 預后相關差異表達基因的篩選

2.3.2 預后相關差異表達基因的驗證 根據TCGA數據庫下載的患者臨床信息分別繪制上述11個DEGs的生存曲線,結果顯示COL9A3,TTYH1,SOX10,MATN3,MAST4,ALCAM,IL12RB2七個基因對TNBC患者預后有顯著影響(圖4)。

圖4 生存曲線驗證預后相關差異表達基因

2.4 預后相關基因的ROC曲線分析同時在單因素Cox比例回歸分析和生存曲線分析結果中均有顯著影響的7個基因,分別繪制ROC曲線進行分析,數據來自于TCGA數據庫的TNBC基因表達譜,結果顯示COL9A3、TTYH1、MATN3、MAST4、IL12RB2的AUC>0.8(圖5)。

圖5 ROC曲線分析

2.5 癌組織中COL9A3的表達與腫瘤浸潤性免疫細胞的相關性分析目前關于三陰性乳腺癌中COL9A3基因的相關研究十分有限,因此我們使用CiberSort算法研究了COL9A3的表達是否與腫瘤浸潤性免疫細胞有關,結果顯示,TNBC中22個免疫細胞相對含量(圖6A)和22個免疫細胞之間的相關性(圖6B)。

圖6 癌組織中COL9A3的表達與腫瘤浸潤性免疫細胞的相關性分析

根據COL9A3表達量的中位數將樣本分為高、低兩組,計算兩組22種免疫細胞浸潤程度的差異,P<0.05為具有統計學意義。小提琴圖結果表明,在TNBC組織中COL9A3高表達組中初始B細胞,記憶B細胞,單核細胞浸潤程度高,而低表達組中靜息態樹突狀細胞,活化的肥大細胞浸潤程度高(圖6C)。

散點圖結果顯示,記憶B細胞(R=-0.22,P=0.049),單核細胞(R=-0.29,P=0.007 1)與COL9A3表達呈負相關,而初始B細胞(R=0.26,P=0.016),中性粒細胞(R=0.22,P=0.05)與COL9A3表達呈正相關(圖6D)。

2.6 COL9A3在正常乳腺組織中的表達單細胞RNA測序結果顯示,COL9A3在正常乳腺組織中主要分布在成纖維細胞和內皮細胞中(圖7),而結果顯示TNBC高表達COL9A3,提示TNBC細胞有發生EMT的傾向,更易侵襲、轉移。

圖7 COL9A3基因ScRNA-seq分析

3 討論

根據目前主流觀點認為,Luminal A型、Luminal B型屬于內分泌治療敏感的乳腺癌亞型,使用內分泌療法能有較為顯著的治療效果, HER2 陽性的乳腺癌病人對HER2靶向治療具有良好的反應性[7],然而由于 TNBC患者缺少 ER、PR和 HER2 的表達,因此TNBC病人往往不能從內分泌治療及抗HER-2的靶向治療中取得良好的療效,TNBC 復雜的生理生化特性導致其在臨床治療中困難重重[8],無法準確評估其預后,尋找確切的預后分子成為了目前亟待解決的問題[9]。

本實驗為了得到TNBC中的DEGs,首先搜集、整理了大量的數據集,并通過GEO2R和R語言對其進行深度挖掘,發現了217個DEGs。在此,我們通過GO和KEGG對這217個基因進行了富集分析,KEGG分析結果顯示這些基因主要在雌激素受體信號通路富集,推測其可能是通過調控此通路促進TNBC的發生發展,這也印證了TNBC對內分泌療法不敏感的特性,GO分析結果顯示這些基因在生物過程中主要富集在乳腺發育,在細胞成分中主要富集在細胞基底外側膜和細胞基膜,可能是其遠處轉移的機制之一[10],此外,值得注意的是這些基因在分子功能中主要富集在跨膜受體蛋白酪氨酸激酶活性,而此蛋白與多個乳腺癌相關信號通路密切相關[11],為三陰性乳腺癌的有效藥物篩選提供參考。為了進一步篩選出預后相關DEGs,使用單因素Cox回歸分析找到了971個乳腺癌預后相關基因,并將其與篩選出的217個差異表達基因取交集,將研究范圍縮小至11個預后相關的DEGs。對這11個基因進行生存分析,結果顯示COL9A3,TTYH1,SOX10,MATN3,MAST4,ALCAM,IL12RB2七個基因對TNBC患者預后有顯著影響。其中COL9A3的ROC曲線AUC>0.8,這提示COL9A3對TNBC預后評估有較好的特異性和敏感度。

在最近幾年里,免疫療法在癌癥治療的各類方法中備受矚目[12],因此,研究并討論了COL9A3的表達對腫瘤浸潤性免疫細胞的影響,結果顯示初始B細胞,記憶B細胞,單核細胞,中性粒細胞與COL9A3表達顯著相關。COL9A3高表達可能通過募集初始B細胞誘導體液免疫,發揮抗體介導的抗腫瘤效應,但與效應性T淋巴細胞的浸潤沒有明顯相關性,提示其對抗腫瘤免疫的影響并不明顯。

單細胞RNA測序可在單個細胞的水平上對基因組,轉錄組進行測序分析,其使用從臨床獲得的樣本在單細胞水平衡量正常組織細胞與癌組織細胞之間的差異性,對腫瘤微環境的研究有極其重要的價值,目前有學者提出EMT并非二元開關式的變化,而是漸進的過程[13],我們發現TNBC高表達COL9A3基因,提示其兼具上皮特性和間充質特性,推測其可能是一種自然存在的早期混合型EMT細胞,介于上皮細胞與間充質細胞兩種狀態之間,這使TNBC侵襲轉移能力增強[14],是其惡性程度高的有力證據之一,并且有研究提出細胞外基質中存在腫瘤相關成纖維細胞,通過多種機制誘導癌細胞侵襲轉移[15]。

雖然在大數據樣本中篩選并探討了COL9A3作為一種預后相關基因對診斷和治療TNBC的潛在價值,但是由于沒有濕實驗的驗證,目前還不能完全確定COL9A3基因是否可以明確指導臨床診斷、治療TNBC,在未來的研究中,我們計劃使用siRNA轉染沉默COL9A3基因,然后從細胞增殖、遷移侵襲等方面去驗證COL9A3基因的具體作用機制,進一步探討其在臨床的應用價值。