莊黎明,丁雅婷,蔡成福,葉輝,蔡耿明

(1.福建醫(yī)科大學附屬泉州市第一醫(yī)院 耳鼻咽喉頭頸外科,福建 泉州 362000;2.泉州市醫(yī)學高等?？茖W校附屬人民醫(yī)院 耳鼻咽喉科,福建 泉州 362000;3.廈門醫(yī)學院附屬海滄醫(yī)院 耳鼻咽喉頭頸外科,福建 廈門 361026;4.福建醫(yī)科大學臨床醫(yī)學部,福建 福州 350122)

頭頸癌(head neck cancer,HNC)是世界上第8種最常見的癌癥。根據2020年全球癌癥統(tǒng)計,HNC的發(fā)病人數約為870 000例,占所有惡性腫瘤的4.5%。HNC的發(fā)病機制是一個多步驟的過程,涉及分子變化的漸進積累。進一步闡明HNC發(fā)展中涉及的分子事件,可能有助于識別潛在和有效的生物標志物,并為靶向治療提供新的方向[1]。頭頸部鱗狀細胞癌(head neck squamous cell carcinoma,HNSCC)是其中發(fā)生率最高的一類。標準療法如保留器官的微創(chuàng)手術方法、放療和多模式治愈性療法等的改進已經可以更好地保留患者的頭頸部器官的功能,并提高患者的5年生存率[2-4],但是復發(fā)性或轉移性患者的療效及預后仍差強人意。基于T細胞的免疫治療,已被證明能提高復發(fā)性或轉移性HNSCC患者的總生存率[5],但只有很低比例的患者對此療法有響應,還存在免疫耐受和免疫逃逸等問題,因此,轉移性、復發(fā)性和晚期HNSCC患者的預后仍然不理想。對HNSCC進展和轉移形成的機制已經有深入的研究,但仍需要開發(fā)新的治療策略[2]。

RNA修飾,例如N6-甲基腺苷(N6-methyladenosine,m6A),在表觀遺傳基因調控和細胞功能中發(fā)揮著不可忽視的作用,與許多人類疾病密切相關。m6A是在真核細胞中記錄的最常見的轉錄后RNA修飾,這些修飾大多發(fā)生在5’-和3’-非翻譯區(qū)(UTRs)以及內部長外顯子的信使核糖核酸終止密碼子附近,在許多生理和病理過程中起著重要作用,包括產生免疫反應、微小RNA的編輯和各種癌癥的進展[6]。目前已有研究表明,m6A修飾失調與膠質母細胞瘤、乳腺癌癥、胃癌、結直腸癌及甲狀腺乳頭狀癌[7]等腫瘤的發(fā)生發(fā)展密切相關。有研究表明[8],喉鱗狀細胞癌(laryngeal squamous cell carcinoma,LSCC)患者的預后,可以通過m6A相關的LncRNA表達量構建的預后模型預測。然而,m6A修飾在HNSCC進展中的作用還需要進一步研究[9]。

1 材料與方法

1.1 分析流程

分析過程流程見圖1。

圖1 數據收集分析過程 注:HNSCC(頭頸部鱗狀細胞癌);GO(基因本體);KEGG(京都基因與基因組);PPI(蛋白相互作用);GSVA(基因集變異分析)。

1.2 數據處理

1.2.1 RNA提取、文庫制備和測序 RNA提取、文庫制備和高通量測序的數據分析均由中國有限公司賽格健康科技有限公司進行。按照Chomczynski等[10]的方法,使用TRIzol(InvitrogenTM,Cat.No.15596018),從臨床收集的3對癌及癌旁組織樣品中提取總RNA。在DNaseI(Thermo ScientificTM, Cat.No.EN0521)提取RNA后進行DNA消化。RNA質量通過用NanodropTMOneC分光光度計(Thermo ScientificTM, Cat.No.ND-ONEC-W)檢測A260/A280來測定。RNA經證實完整性及合格后通過QubitTM3熒光計(InvitrogenTM, Cat.No.Q33216)與QubitTM RNA HS測定試劑盒(InvitrogenTM, Cat.No.Q32855)進行定量。對應于250～500 bp的文庫產物進行富集、定量,并最終在NovaSeq(Illumina?)上測序。

1.2.2 RNA序列數據分析 原始測序數據首先通過SOAPnuke(版本1.6.0)進行過濾,去除低質量及被污染的讀數。Clean Reads首先根據唯一分子標識符(unique molecular identifier,UMI)序列進行聚類,相同UMI序列分組到同一聚類。生成所有子聚類后,進行多序列比對,得到每個子聚類的一個一致序列。

1.2.3 新基因和轉錄物的預處理、比對和分析 使用Guppy軟件(版本5.0.16)[11]對原始讀取進行基礎調用,評估讀取質量,Guppy還對適配器序列進行了修剪。Nanofilt(2.7.1版)過濾低質量(Q值<7)和短長度讀數(<50bp)[12]。使用Fclmr2(版本0.1.2)和短長度RNA-Seq數據校正剩余的干凈讀數[13]。然后,使用Minimap2(2.17-r941版)將每個文庫的干凈讀數與其基因組進行比對[14]。使用Samtools(1.10版)計算干凈讀數與參考基因的比對率[15]。為了減少結果的冗余,合并僅在5’端外顯子差異的比對,并使用StringTie軟件(版本:2.1.4;參數:-保守性-L-R)獲得非冗余轉錄序列。Gffcompare軟件(版本:0.12.1;參數:-R-C-K-M)用于將轉錄物與基因組的已知轉錄物進行比較,并發(fā)現(xiàn)新的轉錄物和新的基因。為了獲得新轉錄物的全面功能信息,根據序列相似性和基序相似性對7個數據庫的轉錄物進行了注釋,包括NR、Pfam、UniProt、京都基因與基因組(kyoto encyclopedia of gengs and genomes KEGG)、基因本體(gene ontology,GO)、KOG/COG和pathway。

1.2.4 基因轉錄物的結構分析 將轉錄本與已知的基因組轉錄本進行比較,以校正邊界。

1.2.5 測序質量控制 對牛津納米孔測序技術(oxford nanopore technologies, ONT)[16]原始測序數據進行質量控制。

1.2.6 數據過濾 根據二代及三代數據過濾標準,對數據進行過濾,從而得到有效數據。

1.2.7 功能注釋和通路富集分析 對上述對照組鑒定的差異表達基因(differentially expressed gene,DEG)進行功能注釋,GO富集分析使用了GO術語的注釋和可視化和元景觀。GO術語的不同表達基因列表之間的重疊通過富集分析圓圖進行。然后將DEG引入FunRich(功能富集分析工具)(http://www.funrich.org/)用于KEGG通路分析。GENEMANIA 用于構建DEG的基因-基因相互作用網絡,以評估這些基因的功能。

1.2.8 蛋白-蛋白相互作用(protein protein interaction,PPI)網絡的構建 交集的基因通過string V11.0進行基因之間的相互作用預測,得到它們之間的關系,篩選標準為combined_score>=0.4。并通過Cytoscape軟件(V3.5.1)進行網絡圖的構建可視化。 共表達的hub基因通過對PPI中的基因進行共表達分析,從而發(fā)現(xiàn)它們之間的表達關系,進而發(fā)現(xiàn)它們在此實驗中的調控關系。

1.2.9 基因集變異分析(gene set variation analysis,GSVA) 使用軟件R(版本3.6)的GSVA包分析m6A。

2 結果

2.1 轉錄組差異基因與差異m6A的差異基因的交集

首先通過癌與癌旁組織對比得到轉錄組差異基因篩選出來,標準為P<0.05且log2>1或者log2<-1。將m6A的差異基因篩選出,標準為P<0.05且log2>1或者log2<-1。通過兩組差異基因篩選出轉錄組上調和m6A上調的基因以及轉錄組下調和m6A下調的基因。我們得到96個上調基因(轉錄組上調m6A上調)和159個下調基因(轉錄組下調m6A下調)。

2.2 交集基因的富集分析

首先將交集基因按照上調和下調基因來做富集分析,從而得到它們各自的方向,具體見表1。

表1 交集基因的富集分析

對這部分的基因進行功能以及pathway的富集分析,對顯著性差異基因功能分別做了柱狀圖(圖2),顯著性篩選的標準為P<0.05。

圖2 部分基因功能及pathway富集分析中顯著性差異的基因 注:縱坐標為差異基因功能名稱, 橫坐標為P值的負對數(-LgP);圖中僅顯示按照-LgP排序最大的各前10項。collagen-containing extracellular matrix(含膠原蛋白的細胞外基質);secretory granule lumen(分泌顆粒腔);cytoplasmic vesicle lumen(細胞質囊腔);vesicle lumen(囊泡腔);blood microparticle(血液微粒子);tertiary granule lumen(三級顆粒腔);microvillus membrane(微濾膜);endoplasmic reticulum lumen(內質網腔);apical part of cell(細胞頂端);apical plasma membrane(頂端質膜);extracellular matrix structural constituent(細胞外基質結構成分);glycosaminoglycan binding(糖胺聚糖結合);heparin binding(肝素結合);extracellular matrix structural constituent conferring compression resistance(細胞外基質中具有抗壓性的結構成分);oxidoreductase activity(氧化還原酶活性);collagen binding(膠原蛋白結合);endopeptidase inhibitor activity(內肽酶抑制劑活性);peptidase inhibitor activity(肽酶抑制劑活性);endopeptidase regulator activity(內肽酶調節(jié)器活性);sulfur compound binding(硫化合物結合);complement activation, alternative pathway(補體激活,替代途徑);O-glycan processing(O型糖加工);iron ion homeostasis(鐵離子穩(wěn)態(tài));body fluid secretion(體液分泌);secretion by tissue(組織分泌);protein O-linked glycosylation(蛋白質 O 鏈接糖基化);transition metal ion homeostasis(過渡金屬離子穩(wěn)態(tài));iron ion transport(鐵離子運輸);regulation of ossification(骨化調節(jié));lung development(肺發(fā)育)。下同。

2.3 差異基因信號通路分析

根據差異基因通路中每一通路的顯著性水平(P值)和所包含的差異基因數據(表2)構建顯著性差異基因通路的分布圖。對部分顯著性差異基因通路做了氣泡圖(圖3)。從通路的結果中可以看出:在癌組織中PPAR通路下調,免疫相關通路下調,以及代謝通路的上調,從而影響了腫瘤的發(fā)展和增殖。

表2 差異基因數據分析

圖3 顯著性差異基因通路氣泡圖 注:縱坐標為差異基因通路名稱,橫坐標為富集程度(enrich factor);mucin type O-glycan biosynthesis(黏蛋白型 O-糖的生物合成);complement and coagulation cascades(補體和凝血級聯(lián));renin secretion(腎素分泌);PPAR signaling pathway (PPAR信號通路);glycosphingolipid biosynthesis-lacto and neolacto series(糖磷脂的生物合成-內酯和新內酯系列);staphylococcus aureus infection(金黃色葡萄球菌感染);ferroptosis(鐵死亡);fatty acid degradation(脂肪酸降解);other types of O-glycan biosynthesis(其他類型的 O-糖生物合成);phosphonate and phosphinate metabolism(膦酸鹽和膦酸鹽代謝);intestinal immune network for IgA production(產生 IgA 的腸道免疫網絡);malaria(瘧疾)。下同。

2.4 PPI的建立

通過string V11.0對上述交集的基因進行基因之間的相互作用預測,得到它們之間的關系,并通過Cytoscape軟件進行網絡圖的構建可視化構建,點用來表示基因,線用來表示基因之間的關系,其中顏色表示上下調的關系和程度,紅色表示上調基因,而綠色表示下調基因,顏色越深表示變化越顯著。

通過構建差異基因的作用網絡,得到基因之間的作用關系。基因的度排序結果見表3。通過癌與癌旁組織進行差異基因分析發(fā)現(xiàn),DCN、GLU、ARG2、C3均呈現(xiàn)下調趨勢,而CDK1、APOE、POSTN、RPL8、YWHAZ、MMP13均呈現(xiàn)上調的趨勢(表3)。

表3 基因的度與作用趨勢

2.5 m6A_GSVA

通過癌和癌旁組織的6個樣本的轉錄組測序技術(RNA sequencing,RNAseq),其每千個堿基的轉錄每百萬映射的Transcript(Transcript per Kilobase per Million mapped reads,TPM)值進行GSVA分析(GSVA:用來體現(xiàn)每個樣本的在某個通路的變化情況)。運用t檢驗方法篩選,通過基因的信號值預測每個樣本的富集得分值,計算兩組之間的富集差異性,得到有差異富集的通路。從圖4可見PI3K/AKT信號通路為在腫瘤中上調,PPAR通路在腫瘤中為抑制通路。

圖4 6個樣本在不同分型的樣本富集得分值 注:圖中顏色是由藍向紅來變化,表示富集分數數值的升高,表示通路被激活。

2.6 共表達的hub基因

通過對PPI中的基因進行共表達分析,從而發(fā)現(xiàn)它們之間的表達關系,進而發(fā)現(xiàn)它們在此實驗中的調控關系。通過計算PPI中的基因之間的共表達關系,具體見圖5。從圖5中PPI共表達關系,我們發(fā)現(xiàn)AGR2、CDK1在m6A調控下,與周圍的基因聯(lián)系很密切,具體結果見表4。

圖5 PPI中基因的共表達關系 注:深藍色表示基因的度越大,而淺色表示基因的度越小。

表4 共表達基因的度

3 討論

由于腫瘤的異質性、對包括免疫治療在內的多種治療效果不佳,以及耐藥性等問題的存在。HNSCC患者,特別是復發(fā)性和中晚期患者的治療效果和預后并不令人滿意,適當的個體化治療方案需要進一步去探索和建立,從而有效提高和改善HNSCC患者的治療效果和預后。既往的研究表明,m6A調節(jié)因子的異常表達可能與多種癌癥的發(fā)生發(fā)展有關,包括乳腺癌癥、膀胱癌癥、膠質瘤和結直腸癌等。與正常組織相比,某些m6A相關基因在HNSCC組織中的表達上調,與免疫細胞浸潤以及免疫檢查點抑制劑(immune checkpoint inhibitors,ICIs)治療的反應和結果有關[6]。從PPI共表達關系圖中,我們發(fā)現(xiàn)在m6A調控下,共表達基因與周圍的基因聯(lián)系很密切,它們處于網絡的中心位置。

細胞周期蛋白依賴性激酶1(cell cycle protein-dependent kinase 1, CDK1)是一種絲氨酸/蘇氨酸激酶,是細胞周期的主調節(jié)因子[17],它對細胞周期調節(jié)至關重要。CDK1在控制DNA復制、修復及mRNA轉錄等功能作用中均能起到一定的作用[18]。一般情況下,細胞周期控制的失調主要引起腫瘤生長[19],CDK1的上調與腫瘤的增殖和進展相關。一項研究表明[20],鐵死亡相關基因與CDK1密切相關。并且有研究已驗證了CDK1 可使LA 相關蛋白 1(LARP1)磷酸化[21]。幾項研究均顯示了LARP1在幾種惡性腫瘤中顯著上調[22-25]。與之相關的是,LARP1可作為CDK1的底物在部分癌癥中上調,從而提高異常翻譯存活率,導致癌細胞增殖增加。載脂蛋白E(apolipoprotein E,ApoE)是一種參與脂質轉運和脂蛋白代謝的多功能蛋白,介導組織和細胞中的脂質分布/再分配。它還可以調節(jié)炎癥和免疫功能,保持細胞骨架的穩(wěn)定性,改善神經組織功能。近年來,ApoE頻繁出現(xiàn)在腫瘤研究中,并逐漸成為腫瘤的生物標志物。研究發(fā)現(xiàn)ApoE在許多腫瘤細胞,包括神經膠質瘤、卵巢腫瘤、口腔鱗狀細胞癌癥、胰腺導管細胞癌等大多數實體瘤組織中高表達。ApoE可以調節(jié)巨噬細胞的極化變化,參與腫瘤免疫微環(huán)境的構建,調節(jié)腫瘤炎癥和免疫反應,在腫瘤侵襲和轉移中發(fā)揮作用。許多研究已經證實ApoE促進腫瘤的發(fā)展、增殖和轉移,但最近的研究也發(fā)現(xiàn)ApoE具有保護作用[26]。Periostin(POSTN)屬于基質細胞蛋白,分泌到細胞外環(huán)境中的非結構ECM蛋白,在大多數成年組織中以低水平表達。這些蛋白質與細胞表面受體相互作用,介導細胞和細胞外通訊。Periostin在各種實體上皮腫瘤中過表達,其與細胞表面受體整合素的相互作用調節(jié)細胞內信號通路,對癌癥有直接影響。它在轉移中上調,并可影響轉移病灶的大小和數量[27]。編碼核糖體蛋白(ribosomal protein L8,RPL8)是60S亞基的組成部分,該蛋白屬于核糖體蛋白的L2P家族,位于細胞質中。RPL8是一種鐵死亡相關基因,它的表達隨胰腺腫瘤分級而顯著改變,與胰腺癌癥預后相關,與膠質瘤癌癥進展有關[28-29]。血小板參與腫瘤血管生成和癌癥進展,是癌癥生物標記物的潛在來源。研究發(fā)現(xiàn),RPL8在癌癥患者的血小板蛋白質組中上調[30]。

簇蛋白是一種進化高度保守的糖蛋白,簇蛋白通過調節(jié)多種信號通路介導各種癌癥(如前列腺癌、乳腺癌、肺癌、肝癌、結腸癌、膀胱癌和癌癥)中的癌癥進展,在調節(jié)幾個基本生理過程中發(fā)揮了突出作用,包括程序性細胞死亡、轉移、侵襲、增殖和細胞生長[31]。簇蛋白介導的線粒體自噬可以抑制口腔癌癥的干性和自我更新能力[32]。它在神經保護和癌癥以及化療耐藥性方面具有重要意義,被認為既是一種生物標志物,也是一種非常可探索的新治療靶點,適用于多種疾病[33]。但簇蛋白在惡性腫瘤中的作用一直存在爭議,它可以促進體外的腫瘤發(fā)生和對化療藥物的耐藥性。然而, 簇蛋白在小鼠體內也可能會限制神經母細胞瘤和癌癥的發(fā)展,研究發(fā)現(xiàn),它在抑制細胞死亡、調節(jié)促生存信號傳導和增強腫瘤細胞化療耐藥性方面至關重要。簇蛋白高表達與膀胱癌患者生存率低明顯相關,被定義為膀胱癌中的一種癌基因。它通過抑制細胞凋亡來增強癌癥對順鉑的耐藥性[34]。蛋白聚糖由鐵死亡的細胞釋放,然后作為警報信號觸發(fā)先天和適應性免疫反應。鐵死亡期間蛋白聚糖的早期釋放是一個活躍的過程,涉及分泌性大自噬/自噬和溶酶體胞吐作用。一旦釋放,細胞外蛋白聚糖與巨噬細胞上的受體高級糖基化終產物特異性受體(AGER)結合,以NFKB/NF-κB依賴的方式觸發(fā)促炎細胞因子的產生[35]。YWHAZ(也稱為14-3-3ζ)是許多信號轉導途徑的中心樞紐蛋白,在腫瘤進展中起著重要作用。越來越多的證據表明,YWHAZ在多種類型的癌癥中經常上調,并在細胞生長、細胞周期、凋亡、遷移和侵襲等廣泛的細胞活動中起致癌基因的作用。此外,YWHAZ可能是幾種癌癥的診斷、預后和化療耐藥性的潛在生物標志物[36]。臨床上, YWHAZ高表達提示胃癌、膀胱尿路上皮癌、乳腺癌、前列腺癌、結腸癌、非小細胞肺癌等多種癌癥患者預后差有關[37-42]。經典的鐵死亡誘導劑erastin促進了軟骨細胞中基質金屬蛋白酶13的表達,同時抑制了Ⅱ型膠原的表達[43-44]。AGR2基因位于染色體7p21.3內,是一種蛋白質二硫鍵異構酶(protein disulfide isomerase,PDI),PDI家族成員作為氧化還原酶而發(fā)揮作用,負責催化二硫鍵的形成或重排,這是發(fā)生于細胞內質網中的一種關鍵蛋白質修飾。許多研究已證實,AGR2高表達于含有黏液分泌功能細胞的組織。腫瘤微環(huán)境(tumour microenvironment,TME)中成分復雜,包含腫瘤細胞、免疫細胞、腫瘤相關炎癥細胞和成纖維細胞等細胞,此外還有組織間質、血管、細胞因子和趨化因子等。頸癌研究中發(fā)現(xiàn),AGR2的表達于與ALDH1、Sox2和Oct4之間存在相關性,表明AGR2可能與頭頸癌有關[45]。此外,AGR2在各種上皮腫瘤中表達上調,通過激活包括p53激活途徑在內的多種途徑促進腫瘤的侵襲性生長、血管生成、侵襲、轉移和耐藥性等惡行生物學行為[46]。

有趣的是,與已有研究相反,測序結果發(fā)現(xiàn)CLU、AGR2在中晚期HNSCC癌組織中表達,我們認為:① 實驗所收集的標本為中晚期HNSCC標本,與上述研究的標本有所不同。②它們可以逃離這個位置,在細胞表面或細胞外基質中發(fā)揮作用。③轉化生長因子-β(transforming growth factor-β,TGF-β)的缺失可以抑制AGR2。因為除了通過關鍵信號系統(tǒng)誘導AGR2外,還可以通過TGF-β的介導調節(jié)AGR2,即TGF-β/AGR2軸。TGF-β的缺失可以抑制AGR2從而可以抑制p53腫瘤抑制活性[20]。④測序顯示 FOXA1低表達,而功能研究表明,FOXA1在前列腺癌中存在兩種相悖的功能:FOXA1通過雄激素受體(androgen receptor,AR)途徑誘導前列腺癌細胞繁殖,同時FOXA1也可以通過不依賴AR途徑抑制癌細胞遷移和上皮-間皮轉化(epithelialto-mesenchymal transition,EMT)。AR及其輔因子FOXA1與一個特定的內含子雄激素響應元件結合,通過MBOAT2調節(jié)前列腺癌的鐵死亡[47]。FOXA1是Yes相關蛋白(yes-associated protein,YAP)依賴性的轉錄因子,在轉錄因子CP2(transcription factor CP2,TFCP2)作為YAP調控轉錄因子的情況下,可以于YAP一起調節(jié)FOXA1,誘導鐵蛋白輕鏈(ferritin light chain,FLT)轉錄,從而抑制鐵死亡[35,48]。

我們測序的數據結果表明,m6A修飾與FOXA1、RPL8、DCN等鐵死亡相關基因緊密相關。鐵死亡是鐵依賴性的,以脂質過氧化為主要特征的一種區(qū)別于凋亡細胞死亡的一種獨特的程序性死亡類型。已在包括癌癥等的自身免疫疾病中被探索研究,并被作為新的治療策略[49-50]。因此我們推測,m6A修飾鐵死亡相關基因在中晚期HNSCC中起重要作用,這需要進一步的實驗深入研究證實。本研究探討了m6A介導的中晚期HNSCC表觀轉錄組學,對中晚期HNSCC的免疫逃逸進行初步的探討,為HNSCC的個體化精準治療提出新的觀點。