代佳佳，肖 何，張 琴，李松霖，陳 川，王 閣

0 引 言

對于局部晚期直腸癌（locally advanced rectal cancer，LARC）患者，術(shù)前新輔助化放療（neoadjuvant radiochemotherapy，NCRT）后進行全直腸系膜切除術(shù)（total mesorectal excision，TME）是目前的標準治療方式［1-3］。相對于術(shù)后輔助放化療，NCRT可降低術(shù)前分期、提高手術(shù)切除率和保肛率，并顯著降低局部復(fù)發(fā)率［4-5］。然而在臨床實踐中，僅10%～30%的患者顯示病理完全緩解［6-8］。因此除術(shù)后TNM分期以外，術(shù)后病理結(jié)果顯示的腫瘤退縮分級（tumor regression grading，TRG）也NCRT效果評價的一項重要參考指標。

美國國立綜合癌癥網(wǎng)絡(luò)（National Comprehensice Cancer Netword，NCCN）指南推薦評估腫瘤治療反應(yīng)的病理學(xué)分級系統(tǒng)包括完全反應(yīng)（TRG4）：無活的癌細胞殘留；中度反應(yīng)（TRG3）：單個或小簇癌細胞殘留；輕度反應(yīng)（TRG2）：殘留癌灶；反應(yīng)不良（TRG1）：僅少數(shù)或未見癌細胞消退［7］。研究表明，術(shù)前NCRT后腫瘤完全退縮（TRG4）和部分退縮（TRG3）與患者無復(fù)發(fā)生存和總生存密切相關(guān)，腫瘤達到完全退縮的患者較退縮情況較差患者具有明顯的生存獲益［9-10］。目前，第7版美國癌癥聯(lián)合委員會（American Joint Commettee on Cancer，AJCC）分期已提出對直腸癌標本檢查時應(yīng)評價新輔助治療后的治療反應(yīng)［11］。因此，對于術(shù)前放化療后能到達TRG 3/4級人群的預(yù)測有助于篩選出對NCRT最大獲益人群和避免對放化療不敏感患者的毒副反應(yīng)，對促進直腸癌的個體化治療有重要的臨床意義。

迄今已有大量文獻報道了各種臨床和分子生物學(xué)指標對NCRT完全病理學(xué)緩解的影響因素及預(yù)測模型的構(gòu)建［12-18］。其中包括治療前腫瘤大小和腸系膜淋巴結(jié)轉(zhuǎn)移狀態(tài)等臨床特征［12］、核磁共振影像組學(xué)特征［13］、腫瘤芽殖等病理形態(tài)特征［14］以及放化療敏感性相關(guān)基因標簽等［15-18］；而基于基因標簽的預(yù)測模型除具有潛在的臨床應(yīng)用價值外，還因與腫瘤放療抵抗機制密切相關(guān)而具有放射生物學(xué)探索的價值。從樣本量考慮，合并多個數(shù)據(jù)集擴大訓(xùn)練集樣本量對于構(gòu)建穩(wěn)定的預(yù)測模型具有重要作用［19］。因此，本研究試圖綜合分析4個局部晚期直腸癌NCRT基因表達譜數(shù)據(jù)集和腫瘤退縮等級，構(gòu)建預(yù)測模型以進一步判斷該類模型的臨床應(yīng)用價值和潛在的直腸癌放化療抵抗機制。

1 資料與方法

1.1 數(shù)據(jù)集的獲取和前處理以“preoperative chemoradiotherapy”和“rectal cancer”或“neoadjuvant chemoradi otherapy”和“rectal cancer”作為檢索詞搜索GEO數(shù)據(jù)庫，選擇可獲取基因微陣列表達值和腫瘤退縮分級的數(shù)據(jù)集，排除RNAseq、非編碼RNA數(shù)據(jù)集和沒有標示退縮等級的數(shù)據(jù)集。共有4個數(shù)據(jù)集符合入組標準，數(shù)據(jù)集名稱和相關(guān)信息見表1。采用R軟件包“annotate”將各微陣列探針I(yè)D號注釋到基因符號。對于多個探針對應(yīng)1個基因，取4分位距最大探針代表該基因表達值。以基因符號合并 GSE35452、GSE46862、GSE68204 3個數(shù)據(jù)集共174個樣本，其中79例達到TRG 3/4，緩解率為44.1%。利用“virtualArray”包函數(shù)“virtualArrayCom-Bat”對合并數(shù)據(jù)集的批次效應(yīng)進行校正［20］。δ統(tǒng)計量和主成分分析判斷批次效應(yīng)校正結(jié)果［21］。

1.2 基于LASSO的Logistic回歸和SVM法篩選候選基因和模型構(gòu)建R軟件包“sample”函數(shù)隨機將174例合并數(shù)據(jù)集按7∶3比例分成訓(xùn)練集（n=121）和驗證集（n=53）。對訓(xùn)練集采用單因素Logistic回歸篩選與TRG 3/4顯著相關(guān)的基因，依據(jù)回歸系數(shù)顯著性P值對基因進行排秩，P＜0.05作為納入最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator，LOSSA）算法的候選基因。LASSO算法構(gòu)建預(yù)測模型，取λ對應(yīng)內(nèi)部交叉驗證錯配率最小值對應(yīng)變量作為模型變量，并在驗證集中判斷該模型對TRG 3/4的診斷準確性、特異性和靈敏度。在訓(xùn)練集中對每個基因利用t檢驗分析TRG3/4與TRG0/1/2組間表達值差異，并根據(jù)顯著性P值對基因進行排秩，取前50個基因利用支持向量機（support vector machine，SVM）構(gòu)建預(yù)測模型，在驗證集中判斷該模型對TRG 3/4的診斷準確性、特異性和靈敏度，SVM算法各參數(shù)使用缺失值。以上2種分析分別反復(fù)隨機取樣500次，以判斷模型預(yù)測作用的穩(wěn)定性；并對基因排秩取平均值和標準差以確定排秩前100基因的可重復(fù)性。LASSO算法和SVM分別采用R軟件包“glmnet”和“penalizedSVM”完成。

LASSO算法納入模型最多的前30個基因與獨立驗證集GSE53781有重疊的基因共21個，再次在174例合并數(shù)據(jù)集中驗證該21個基因?qū)Σ±韺W(xué)緩解預(yù)測效能。同時在GSE53781中采用這21個基因單因素Logistic回歸系數(shù)與表達值乘積和作為化放療敏感指數(shù)（sensitivity index，SI）以判斷該指數(shù)對病理學(xué)緩解的預(yù)測作用。

1.3 基于差異表達基因和轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)分析和基因功能注釋在174例合并數(shù)據(jù)集中。利用“l(fā)imma”包分析TRG3/4與TRG0/1/2組間差異表達基因。將多重校正前P＜0.05基因納入轉(zhuǎn)錄網(wǎng)絡(luò)調(diào)控分析，判斷差異表達基因潛在調(diào)控機制［22］?！皉eactome”包用于候選基因的功能注釋。

2 結(jié) 果

2.1 在 GSE35452、GSE46862、GSE68204構(gòu)成的數(shù)據(jù)集中構(gòu)建和驗證預(yù)測模型合并3個數(shù)據(jù)集后共12 803個基因納入分析。批次效應(yīng)校正前后δ統(tǒng)計量分別為2.350（P=0.008）和0.0257（P=1.000），主成分分析亦表明批次效應(yīng)校正后3個數(shù)據(jù)集來源病例在前2個主成分中無聚集性分布，見圖1。單因素Logistic回歸表明在500次隨機抽樣生成的訓(xùn)練集中，回歸系數(shù)顯著性P值最小的前100個基因P值均較大，僅3個基因P＜0.01，且有較大的變異幅度，見圖2。同樣，基因排秩靠前的100個基因平均秩也較大，并存在較大的變異幅度。LOSSA模型在訓(xùn)練集中的診斷準確性、特異性和靈敏度分別為0.998（95%CI：0.991～1.000）、1.000（95%CI：1.000～1.000）、0.998（95%CI：0.981～1.000）。但在驗證集中卻得到較弱的診斷效能，準確性、特異性和靈敏度分別為 0.523（95%CI：0.396～0.642）、0.578（95%CI：0.373～0.762）、0.464（95%CI：0.258～0.700）。同樣，SVM法在驗證集中的準確性、特異性和靈敏度分別為 0.504（95%CI：0.377～0.623）、0.596（95%CI：0.393～0.830）、0.405（95%CI：0.182～0.650）。見圖3。

表1 所用數(shù)據(jù)庫名稱和微陣列芯片Table 1 Series of GEO accession and platforms used for the analysis

圖1 GSE35452、GSE46862、GSE68204數(shù)據(jù)集合并后批次效應(yīng)校正前后主成分分析Figure 1 Principal component analysis before and after batch effect correction after merging three data sets

圖2 單因素Logistic回歸系數(shù)顯著性最高的前100個基因平均P值及排秩Figure 2 Average P values and ranks of the first 100 genes with the highest significance of coefficients from univariate Logistic regression in the training set

2.2 GSE53781作為外部驗證集驗證候選基因在500次訓(xùn)練集中引入LASSO算法模型次數(shù)最多的前30個基因中，有21個與GSE53781重疊。取該21個基因作為獨立驗證的候選基因。LASSO算法表明這21個基因在合并數(shù)據(jù)集中對病理學(xué)完全緩解和接近病理學(xué)完全緩解病例的診斷準確性、特異性和靈敏度分別為 0.816、0.842、0.785。其調(diào)諧參數(shù)λ與對應(yīng)的錯配率見圖4。敏感指標病理學(xué)緩解診斷的AUC為0.863（95%CI：0.811～0.912）。21個候選基因符號和功能注釋以及回歸系數(shù)見表2。GSE53781數(shù)據(jù)集中有26例，其中10例達到病理學(xué)緩解，該組人群中計算得到的敏感指標1.674～54.985，中位值為7.133）。敏感指標有較高診斷效能（AUC=0.925，95%CI：0.817～1.000），見圖4。但是，對GSE53781隨機取出21個基因重復(fù)分析1000次得到對病理學(xué)緩解診斷的AUC 95%CI：0.769～0.969，21基因?qū)?yīng)的AUC處于該CI內(nèi)。

圖3 不同方法對新輔助化放療緩解診斷的特異性和靈敏度散點圖Figure 3 Scatter plots showing the sensitivity and specificity of the predictive model for response to neoadjuvant chemoradiotherapy in validation sets in 500 times of random sampling

圖4 合并數(shù)據(jù)集中21基因LASSO算法λ與內(nèi)部交叉驗證錯配率關(guān)系圖Figure 4 Relationship between internal cross-validation mismatch rate and lambda in LASSO algorithm with 21 candidate genes in the 174 merged datasets

2.3 LARC化放療抵抗基因參與調(diào)控網(wǎng)絡(luò)在174例合并數(shù)據(jù)集12803個基因中，差異表達基因共505個。RNEA分析表明，這些差異表達基因主要受 NF-κB、CTNNB1、IL6/STAT3、E2F1、GLI2 以及MYC轉(zhuǎn)錄因子調(diào)控基因，見圖5。

表2 21個候選基因在合并數(shù)據(jù)集中Logistic回歸系數(shù)和部分功能注釋Table 2 Logistic regression coefficients and functional annotations of 21 candidate genes in the 174 merged datasets

圖 5 GSE35452、GSE46862、GSE68204 3個合并數(shù)據(jù)集中RNEA揭示緩解組與未緩解組間505個差異表達基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)圖Figure 5 Transcription regulatory network of the 505 differentially expressed genes in the response and non-response groups in the 174 merged datasets

3 討 論

目前已有大量基于高通量基因組學(xué)和生物信息學(xué)分析的腫瘤放射內(nèi)在敏感性研究試圖建立個體化放療指標［23］。對于局部晚期直腸癌，已有若干研究篩選出了可構(gòu)建有效的NCRT病理學(xué)緩解或完全緩解預(yù)測模型的候選基因［15-18］。但這些基于統(tǒng)計量排秩而選擇“重要”基因納入模型構(gòu)建的方法易受訓(xùn)練集樣本量的影響。足夠樣本量的訓(xùn)練集才能得到穩(wěn)定的標簽基因和預(yù)測效能。合并多個數(shù)據(jù)集的綜合分析可能具有提高構(gòu)建預(yù)測模型穩(wěn)定性的優(yōu)勢［19，24］。本研究合并了 GEO 數(shù)據(jù)庫中的GSE35452、GSE46862、GSE68204 3個數(shù)據(jù)集，共納入了12 803個基因數(shù)據(jù)值，通過在緩解組和未緩解組間采用Logistic回歸和t檢驗，得出了平均排秩最小的100個基因，但這100個基因在121例訓(xùn)練集中并不穩(wěn)定；而且在不同次取樣構(gòu)成的訓(xùn)練集中，同一基因的排秩存在較大偏差，結(jié)合在各次取樣中對53例內(nèi)部驗證集較弱的預(yù)測效能。該結(jié)果提示預(yù)測模型具有一定的預(yù)測能力，但驗證效能不足，尚需進一步深入研究。有文獻報道其他二分類臨床結(jié)局預(yù)測模型僅使用了 120 例樣本已足夠［19，24］；但本研究單次用174例作為訓(xùn)練集也不足以得到穩(wěn)定的標簽基因，且在不同的研究中得到的標簽基因種類很少重疊［15-18］。提示了參與直腸癌NRCT抵抗或敏感基因種類在人群間存在較大異質(zhì)性，其機制仍需進一步探討。研究表明，在結(jié)直腸癌分子分型基礎(chǔ)上進行預(yù)后標志物篩選，或者聯(lián)合腫瘤細胞內(nèi)在基因與腫瘤微環(huán)境組成細胞如腫瘤相關(guān)成纖維細胞、各種腫瘤浸潤免疫細胞亞群基因表達構(gòu)建療效預(yù)測模型是克服腫瘤異質(zhì)性的有效手段［25-26］。

本研究篩選出的21個預(yù)測基因可能與直腸癌放療和化療敏感性有密切聯(lián)系，涉及糖代謝、絲氨酸/甘氨酸合成與代謝、促凝血機制、轉(zhuǎn)錄調(diào)控和細胞周期調(diào)控等。如ST8SIA5基因?qū)儆谔腔D(zhuǎn)移酶家族29，參與合成多種神經(jīng)節(jié)糖苷。同屬一個家族的另一個酶ST8SIA1則具有維持三陰乳腺癌細胞干性特征，并通過FAK-AKT-mTOR促進轉(zhuǎn)移的作用［27］。AMT基因編碼氨基甲基轉(zhuǎn)移酶是構(gòu)成甘氨酸解理系統(tǒng)蛋白之一。線粒體絲氨酸/甘氨酸代謝系統(tǒng)對快速增殖腫瘤細胞所需的嘌呤核苷合成具有重要作用。而甘氨酸解理系統(tǒng)可清除該過程中產(chǎn)生的過量甘氨酸，否則會在胞內(nèi)會產(chǎn)生氨基酮和丙酮醛等細胞毒性物質(zhì)［28］。而甘氨酸自身也可通過其受體抑制內(nèi)皮細胞血管生成和結(jié)直腸癌細胞轉(zhuǎn)移［29］。快速增殖的腫瘤細胞很容易造成乏氧和瘤體中心區(qū)域的壞死。推測AMT可能對于維持直腸癌在快速增殖情況下改善乏氧并誘導(dǎo)腫瘤轉(zhuǎn)移有重要作用。因此，本研究結(jié)果提示AMT基因表達越高的直腸癌患者達到病理學(xué)緩解的概率越低。

此外，本研究發(fā)現(xiàn)凝血酶（F2）表達可能對放療抵抗有重要影響。凝血酶介導(dǎo)的纖維蛋白/血小板凝聚被認為是腫瘤細胞微轉(zhuǎn)移灶形成的機制之一。凝聚的血小板通過TGF-β1減弱自然殺傷細胞的功能，使微轉(zhuǎn)移腫瘤細胞逃避免疫監(jiān)視［30］。大量研究提示，上皮間質(zhì)轉(zhuǎn)化或干性轉(zhuǎn)化的腫瘤細胞對放療和化療藥物具有較強的抵抗作用［31-32］。而本研究中，對緩解與未緩解組間差異表達基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)分析也提示這些差異表達基因主要受IL-6/STAT3、NF-kB和MYC等與腫瘤轉(zhuǎn)移或干性表型等相關(guān)轉(zhuǎn)錄因子調(diào)控。如IL-6不但可促進結(jié)直腸癌細胞系HCT116和HT29上皮間葉轉(zhuǎn)化和血管擬態(tài)的形成，且與轉(zhuǎn)移密切相關(guān)，同時還可增強A549和H157非小細胞肺癌細胞系CD133+干細胞對放療誘導(dǎo)雙鏈斷裂的修復(fù)能力和抵抗凋亡的能力［31-32］。臨床研究也表明，浸潤能力和淋巴結(jié)轉(zhuǎn)移也是直腸癌病理學(xué)完全緩解的不良預(yù)測因素［12，14］。因此，結(jié)合既往報道可以推測直腸癌中參與上皮間質(zhì)轉(zhuǎn)化和干性轉(zhuǎn)化機制并促進轉(zhuǎn)移的腫瘤微環(huán)境就是參與放化療抵抗的重要機制。

本研究納入分析的4個GEO數(shù)據(jù)集采用了不同的病理學(xué)緩解分級標準、所用雜交芯片檢測平臺不盡相同，微整列芯片重疊基因存在有限性以及未分細胞亞群分析等因素均成為本研究的局限性。盡管本研究并未建立明確的預(yù)測模型，分析結(jié)果提示直腸癌參與放化療抵抗的機制可能與其浸潤、轉(zhuǎn)移機制相關(guān)，并且強調(diào)了參與放化療抵抗基因在人群間的異質(zhì)性和基于腫瘤及其微環(huán)境多種細胞亞群研究的重要性。