一種面向癌癥個(gè)性化的代謝分析方法

馮昌奇 鄭浩然

0 引言

腫瘤異質(zhì)性是癌癥治療的主要障礙之一，它可以發(fā)生在多個(gè)層面，包括基因、分子和組織學(xué)等[1-3]。由于基因的改變，癌細(xì)胞之間在疾病過程中出現(xiàn)了表型和功能上的差異[4]。同樣，腫瘤細(xì)胞與其微環(huán)境的相互作用或血管生成和缺氧的局部變化在腫瘤中也是不均勻的。這種巨大的生物、細(xì)胞和組織異質(zhì)性存在于瘤內(nèi)水平(一個(gè)腫瘤內(nèi)的分子差異)、患者內(nèi)水平(一個(gè)患者病變間的腫瘤特征差異)和患者間水平(患者間腫瘤特征的差異)。這種異質(zhì)性影響了腫瘤的侵襲性和耐藥性，是有效治療策略設(shè)計(jì)的一個(gè)重大挑戰(zhàn)[5]。

細(xì)胞代謝在癌癥[6]和其他各種疾病[7]的發(fā)病機(jī)制中起著重要的作用。癌細(xì)胞重新規(guī)劃其能量代謝，以滿足增殖和遷移的能量需求。侵襲和轉(zhuǎn)移的機(jī)制是復(fù)雜的，死亡率主要是由癌癥進(jìn)展到轉(zhuǎn)移狀態(tài)[8]引起的。了解正常細(xì)胞和癌細(xì)胞之間的代謝差異可以闡明疾病的發(fā)展機(jī)制，也可能有助于確定治療性代謝靶點(diǎn)。

在系統(tǒng)生物學(xué)領(lǐng)域，基因組尺度代謝模型(genome-scale metabolic model,GEM)整合了可用的組學(xué)數(shù)據(jù)和基因組序列，以提供對(duì)有機(jī)體細(xì)胞內(nèi)代謝的更好的機(jī)制理解。將轉(zhuǎn)錄組學(xué)數(shù)據(jù)整合到基因組規(guī)模的代謝網(wǎng)絡(luò)中，以執(zhí)行網(wǎng)絡(luò)級(jí)別的模擬，這是一個(gè)有用的步驟，可以將數(shù)據(jù)規(guī)范化，并試圖從轉(zhuǎn)錄組中表達(dá)的酶的綜合證據(jù)推斷代謝狀態(tài)。GEM可以作為一種強(qiáng)大的工具來對(duì)患者進(jìn)行分層，或確定癌癥代謝中的關(guān)鍵中介基因，使生物學(xué)途徑相互交織，以及提供癌細(xì)胞能量和生物合成需求的代謝物[9-10]。此外，GEM還可用于識(shí)別新的靶點(diǎn)和藥物，并在藥物發(fā)現(xiàn)或藥物重新定位的背景下創(chuàng)建和/或測(cè)試假設(shè)[11-12]。

目前已有大量的研究關(guān)注癌癥代謝上的異質(zhì)性。例如，Berndt等[13]應(yīng)用系統(tǒng)生物學(xué)方法對(duì)10例人類肝細(xì)胞癌(hepatocellular carcinoma,HCC)的代謝變化進(jìn)行了分析，發(fā)現(xiàn)最大的腫瘤間異質(zhì)性是乳酸的攝取和釋放以及細(xì)胞糖原含量的大小。Baloni等[14]對(duì)來自癌癥基因組圖譜(The Cancer Genome Atlas,TCGA)的1 156個(gè)乳腺正常和腫瘤樣本的基因表達(dá)數(shù)據(jù)進(jìn)行了散度分析，構(gòu)建了個(gè)性化的特異性代謝網(wǎng)絡(luò)，根據(jù)它們的不同代謝情況將樣本劃分成了4個(gè)類群，并識(shí)別出了兩種最具抗腫瘤活性的藥物。但目前有關(guān)腫瘤代謝異質(zhì)性的研究大都關(guān)注部分通路的代謝變化或腫瘤代謝系統(tǒng)的異常代謝，很少關(guān)注腫瘤細(xì)胞在整個(gè)代謝系統(tǒng)上發(fā)生的代謝變化之間的異質(zhì)性。

為了探究不同癌癥患者在代謝系統(tǒng)改變上的異質(zhì)性，本文基于679例患者的癌細(xì)胞和正常細(xì)胞的基因表達(dá)數(shù)據(jù)，構(gòu)造了679個(gè)個(gè)性化的環(huán)境特異性代謝網(wǎng)絡(luò)模型。通過對(duì)這些模型進(jìn)行研究，課題組希望分析不同個(gè)體在代謝變化上的獨(dú)特性和普遍性，以便根據(jù)需要設(shè)計(jì)藥物治療方案。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)源

本文從基因組數(shù)據(jù)共享(GDC)門戶網(wǎng)站(https://gdc.cancer.gov/)下載了TCGA中癌癥樣本的RNA-seq數(shù)據(jù)，其中包含了來自33種癌癥類型的60 483個(gè)基因的預(yù)處理基因表達(dá)數(shù)據(jù)。本文將研究對(duì)象定在其中的17種癌癥類型上，因?yàn)樗鼈儼俗銐驍?shù)量的正常細(xì)胞樣本，利于比較癌癥細(xì)胞與原發(fā)組織細(xì)胞之間的代謝差異。為了方便進(jìn)行后續(xù)的差異分析，本文使用TCGA數(shù)據(jù)庫level3層面中的read count數(shù)據(jù)。

TCGA的每個(gè)樣本都有自己的編號(hào)作為唯一標(biāo)識(shí)，其中包含了項(xiàng)目名稱、組織來源、參與者編號(hào)和樣本類型等信息。為了對(duì)不同的患者進(jìn)行個(gè)性化的分析，本文根據(jù)參與者編號(hào)對(duì)樣本進(jìn)行匹配，篩選出同時(shí)具有癌癥細(xì)胞和正常細(xì)胞基因表達(dá)數(shù)據(jù)的樣本。最終，本文篩選出了679對(duì)樣本，樣本數(shù)的詳細(xì)分布見表1。

表1 樣本在17種癌癥類型中的數(shù)量分布Table 1 The number distribution of samples in 17 cancer types

1.2 差異表達(dá)分析

差異表達(dá)分析以具有生物學(xué)意義的方式計(jì)算基因表達(dá)，然后通過對(duì)基因表達(dá)量進(jìn)行統(tǒng)計(jì)學(xué)的分析，找到具有顯著差異的基因。在本文的分析中，對(duì)每一對(duì)來自同一患者的癌癥和正常細(xì)胞樣本，使用Limma[15-16]R包對(duì)其進(jìn)行差異表達(dá)分析，以log2FC > 1.5和P值<0.01作為顯著差異的閾值，篩選出在癌癥細(xì)胞中表達(dá)水平顯著上調(diào)的基因。這些基因代表了癌細(xì)胞相對(duì)其原發(fā)組織顯著活躍的部分。

1.3 生成環(huán)境特異性代謝網(wǎng)絡(luò)

在本研究中，課題組使用當(dāng)前最新的人類基因組規(guī)模代謝網(wǎng)絡(luò)模型Recon3D[17]進(jìn)行癌癥患者的個(gè)性化環(huán)境特異性代謝網(wǎng)絡(luò)模型的生成。Recon3D中包含了2 248個(gè)代謝基因和10 600個(gè)代謝反應(yīng)，并且提供了基因-蛋白-反應(yīng)(gene-protein-reaction，GPR)關(guān)聯(lián)用于整合基因和反應(yīng)之間的關(guān)系。

在這一步中，本文首先從差異表達(dá)分析中獲得的差異上調(diào)基因中篩選出Recon3D中所包含的基因，之后結(jié)合GPR將這些基因映射到代謝反應(yīng)上，作為生成特異性代謝網(wǎng)絡(luò)模型的核心反應(yīng)集合，代表了癌癥細(xì)胞相對(duì)其原發(fā)組織顯著活躍的那部分反應(yīng)。特別地，為了保證生成的特異性代謝網(wǎng)絡(luò)能夠進(jìn)行正常的生物質(zhì)合成，本文中將生物質(zhì)合成反應(yīng)添加到核心反應(yīng)集合中。

之后，本文使用COBRA Toolbox v3.0[18]實(shí)現(xiàn)的FASTCORE[19]算法，基于Recon3D生成環(huán)境特異性代謝網(wǎng)絡(luò)模型。根據(jù)FASTCORE算法的特性，生成的模型能夠在包含核心反應(yīng)集合的條件下，滿足最小通量一致性。生成的代謝網(wǎng)絡(luò)模型代表了特定患者的癌細(xì)胞相對(duì)其正常細(xì)胞發(fā)生異變的那部分代謝系統(tǒng)，不同模型之間的差異充分體現(xiàn)出了癌癥患者之間代謝異常的獨(dú)特性。

1.4 鑒別特異性代謝網(wǎng)絡(luò)的潛在靶點(diǎn)

在本研究中，潛在靶點(diǎn)的鑒別分為兩步：基因敲除和毒性測(cè)試。

在基因敲除階段，本課題組使用RAVEN[20]工具箱提供的removeGenes函數(shù)模擬基因缺失后的網(wǎng)絡(luò)。對(duì)特異性網(wǎng)絡(luò)中的基因，每次模擬刪除其中的一個(gè)，然后計(jì)算敲除后模型的生物質(zhì)合成速率。如果一個(gè)基因被敲除后導(dǎo)致模型的生物質(zhì)合成速率變成0或幾乎為0，那么認(rèn)為這個(gè)基因?qū)λ诘奶禺愋跃W(wǎng)絡(luò)模型來說是必要的，其缺失會(huì)導(dǎo)致模型不能合成生長(zhǎng)所需要的物質(zhì)，進(jìn)而無法維持生存。

將一個(gè)基因作為潛在的靶點(diǎn)，除了要滿足其缺失會(huì)影響到特異性網(wǎng)絡(luò)生存的條件外，還應(yīng)該對(duì)正常細(xì)胞是無害的。所以，在毒性測(cè)試階段，本文從代謝圖譜(metabolic atlas)[21]下載了83個(gè)正常細(xì)胞的代謝網(wǎng)絡(luò)模型，在這些網(wǎng)絡(luò)上模擬必要基因的缺失。通過檢查缺失后的模型能否完成包括能量與氧化還原、內(nèi)轉(zhuǎn)化、底物利用和產(chǎn)物生物合成等在內(nèi)的56個(gè)代謝任務(wù)[22]來確定該癌癥必要基因能否作為潛在的治療靶點(diǎn)。

2 結(jié)果

2.1 差異表達(dá)基因

本文中對(duì)每個(gè)成對(duì)樣本分別進(jìn)行了差異表達(dá)分析，以獲得特定患者的癌細(xì)胞中顯著上調(diào)的基因。679例患者的差異表達(dá)基因數(shù)量分布如圖1(a)所示?？梢钥闯?，不同患者癌細(xì)胞相對(duì)正常細(xì)胞的差異表達(dá)基因數(shù)量具有很大的差異，無論是上調(diào)還是下調(diào)。經(jīng)過課題組的統(tǒng)計(jì)，發(fā)現(xiàn)在多數(shù)樣本的結(jié)果中，差異下調(diào)基因的數(shù)量多于上調(diào)基因，這可能意味著癌細(xì)胞普遍會(huì)更多地抑制正常細(xì)胞基因的表達(dá)。特別地，沒有任何一個(gè)基因出現(xiàn)在全部樣本的上調(diào)基因或下調(diào)基因集合中，并且平均每對(duì)樣本有2.5個(gè)下調(diào)基因和1.5個(gè)上調(diào)基因是獨(dú)一無二的，不在其他任何樣本的差異表達(dá)基因中出現(xiàn)，這些信息充分體現(xiàn)了癌癥在不同患者體內(nèi)發(fā)生改變的獨(dú)特性，這種獨(dú)特性是研究普適性治療方案的主要障礙，有必要得到生物醫(yī)學(xué)研究人員的重視。

圖1 差異表達(dá)基因數(shù)量Figure 1 Number of differentially expressed genes

為了更好地分析不同癌癥類型之間在差異表達(dá)基因上的不同，本文計(jì)算了每種癌癥類型所有患者的平均差異表達(dá)基因數(shù)量，結(jié)果如圖1(b)所示?？梢钥闯觯琇IHC、PRAD、THCA 3種類型的差異表達(dá)基因數(shù)量明顯比其他癌癥要少，這意味著這3種癌癥相對(duì)其他癌癥來說，發(fā)生改變的程度更小。而相反的，LUSC、STAD和UCEC發(fā)生的變化更加明顯。

2.2 個(gè)性化特異性代謝網(wǎng)絡(luò)模型

本文分別將679對(duì)樣本進(jìn)行差異表達(dá)分析得到的差異上調(diào)基因與人類基因組規(guī)模代謝網(wǎng)絡(luò)模型Recon3D結(jié)合，獲取了679個(gè)個(gè)性化的環(huán)境特異性代謝網(wǎng)絡(luò)模型，這些模型充分反映了每位患者癌細(xì)胞相對(duì)其原發(fā)組織細(xì)胞發(fā)生變化的代謝系統(tǒng)。679個(gè)代謝網(wǎng)絡(luò)模型的基因和反應(yīng)規(guī)模如圖2所示。

圖2 679個(gè)個(gè)性化代謝網(wǎng)絡(luò)的規(guī)模Figure 2 The scale of 679 personalized metabolic networks

相比于包含2 248個(gè)代謝基因和10 600個(gè)代謝反應(yīng)的Recon3D而言，個(gè)性化的特異性代謝網(wǎng)絡(luò)平均只包含1 217個(gè)代謝基因和2 181個(gè)代謝反應(yīng)，在去除了多余信息的同時(shí)，充分反映了不同患者癌細(xì)胞發(fā)生改變的代謝系統(tǒng)的獨(dú)特性。

本文統(tǒng)計(jì)了每種癌癥類型的個(gè)性化代謝網(wǎng)絡(luò)模型的平均規(guī)模，發(fā)現(xiàn)LIHC、PRAD和THCA 3種類型的網(wǎng)絡(luò)平均規(guī)模最小，這與差異表達(dá)基因得到的結(jié)果是一致的，即這3種癌癥的癌細(xì)胞相對(duì)正常細(xì)胞發(fā)生的代謝變化程度最小。

2.3 潛在靶點(diǎn)基因

從代謝網(wǎng)絡(luò)模型中刪除一個(gè)基因可能會(huì)對(duì)整個(gè)代謝系統(tǒng)的生長(zhǎng)產(chǎn)生巨大的干擾，也可能幾乎沒有影響。為了分析代謝系統(tǒng)中可能引發(fā)干擾的基因，本文對(duì)679個(gè)個(gè)性化代謝網(wǎng)絡(luò)模型進(jìn)行了基因缺失的測(cè)試。如果一個(gè)基因被敲除后會(huì)影響其所在模型的生長(zhǎng)，那么認(rèn)為這個(gè)基因是該模型的必要基因。如果同時(shí)該基因的缺失又不會(huì)影響到正常組織細(xì)胞的正常代謝活動(dòng)，即認(rèn)為這個(gè)必要基因的缺失對(duì)正常細(xì)胞是無影響的，可以作為候選的治療靶點(diǎn)。

本文對(duì)679個(gè)模型分別進(jìn)行了基因敲除和毒性測(cè)試，并計(jì)算了每種癌癥類型下的模型平均必要基因和潛在靶點(diǎn)數(shù)量，如圖3所示。特別地，PRAD的平均必要基因和潛在靶點(diǎn)數(shù)量顯著高于其他癌癥類型。結(jié)合前面差異表達(dá)基因和代謝網(wǎng)絡(luò)模型的結(jié)果，本課題組考慮這是因?yàn)镻RAD相對(duì)其他癌癥來說，發(fā)生變化的那部分代謝網(wǎng)絡(luò)規(guī)模更小，因而沒有更多的同工酶和同工反應(yīng)來保證生物質(zhì)的合成，因此其網(wǎng)絡(luò)結(jié)構(gòu)在穩(wěn)定性上較差。

圖3 17種癌癥類型的平均必要基因和靶點(diǎn)數(shù)量Figure 3 The average numbers of necessary genes and targets for 17 cancer types

本文統(tǒng)計(jì)了679個(gè)模型的潛在靶點(diǎn)，得到227個(gè)不重復(fù)的基因，這些基因的缺失會(huì)影響一個(gè)到若干個(gè)模型中生物質(zhì)的合成速率。227個(gè)靶點(diǎn)基因在679個(gè)模型中出現(xiàn)次數(shù)的分布如圖4所示?？梢钥闯?，有93%(211/227)的基因在少于100個(gè)模型中被認(rèn)為是潛在的靶點(diǎn)，其中絕大多數(shù)(195/211)出現(xiàn)在679個(gè)模型中的出現(xiàn)次數(shù)少于30次。這意味著對(duì)一例患者可行的治療靶點(diǎn)對(duì)其他患者大概率是無效的，從而導(dǎo)致大多數(shù)的藥物治療方案可能只能對(duì)極少數(shù)的患者具有治療意義。

圖4 227個(gè)靶點(diǎn)基因在679個(gè)模型中出現(xiàn)次數(shù)的分布Figure 4 Distribution of the occurrence number of the 227 target genes in 679 models

特別地，69.6%(158/227)的基因只在少于10個(gè)模型中被認(rèn)為是潛在的靶點(diǎn)，基于這些基因設(shè)計(jì)的藥物治療方案將具有極高的特異性。其中，又有63個(gè)基因的缺失只在某一個(gè)特定的模型中影響生物質(zhì)的合成，這些基因被認(rèn)為是針對(duì)特定患者的特異性治療靶點(diǎn)。這些特異性治療靶點(diǎn)中的大多數(shù)已經(jīng)被各項(xiàng)研究確定為可能的癌癥藥物靶點(diǎn)或被證實(shí)與癌癥的發(fā)生、進(jìn)展和轉(zhuǎn)移有關(guān)聯(lián)。例如，B3GTN3在胰腺癌[23]和宮頸癌[24]的進(jìn)展和轉(zhuǎn)移中發(fā)揮了顯著的作用，與不良預(yù)后相關(guān)，被預(yù)測(cè)為這兩種癌癥的治療靶點(diǎn)。IDH2的突變?cè)谏窠?jīng)膠質(zhì)瘤中被發(fā)現(xiàn)，但這種突變并不常見，研究者們分析了87例神經(jīng)膠質(zhì)瘤病例，只有1例發(fā)生了IDH2突變[25]，間接體現(xiàn)了IDH2在癌癥治療中的特異性。ME1的過表達(dá)與乳腺癌腫瘤體積大、生存期差和化療耐藥性成正相關(guān)，被預(yù)測(cè)為乳腺癌潛在的預(yù)后標(biāo)志物和治療靶點(diǎn)[26]。有研究利用miR-885-5p抑制了胃癌中ME1的過表達(dá)，有效地抑制了癌細(xì)胞的增殖和發(fā)展，改善了預(yù)后不良的狀況[27]。NNMT的表達(dá)在卵巢癌[28]和胃癌[29]中顯著升高，并與術(shù)后預(yù)后不良相關(guān)，被認(rèn)為是卵巢癌治療的重要靶點(diǎn)[30]。PHGDH在乳腺癌和胰腺癌細(xì)胞的增殖、遷移和侵襲中發(fā)揮調(diào)節(jié)作用，可能是乳腺癌和胰腺癌重要的預(yù)后指標(biāo)和治療靶點(diǎn)[31-32]。有研究通過抑制劑下調(diào)PHGDH的表達(dá)，顯著抑制了膀胱癌細(xì)胞的增殖能力并誘導(dǎo)細(xì)胞的凋亡[33]。

63個(gè)特異性靶點(diǎn)基因中大部分已經(jīng)被先前的研究證實(shí)可以作為癌癥的潛在靶點(diǎn)或與癌癥密切相關(guān)，沒有被證實(shí)的部分則代表了日后進(jìn)一步開展癌癥個(gè)性化分析和治療的潛在研究對(duì)象。有趣的是，本文中發(fā)現(xiàn)有14個(gè)特異性靶點(diǎn)都屬于溶質(zhì)運(yùn)載蛋白家族，這可能意味著這一家族的基因在癌癥進(jìn)展中發(fā)揮了重要作用，其研究?jī)r(jià)值有待深入發(fā)掘。

3 討論與結(jié)論

在本文中，課題組使用來自TCGA的17種癌癥類型的679例患者的癌癥細(xì)胞樣本和正常組織細(xì)胞樣本，分析了不同癌癥患者的癌細(xì)胞相對(duì)正常組織細(xì)胞在代謝表型和代謝網(wǎng)絡(luò)層面上發(fā)生的變化。本文的結(jié)果表明，無論是在代謝表型上還是在更深層的代謝網(wǎng)絡(luò)上，不同的癌癥患者細(xì)胞癌變后發(fā)生的代謝異常都在很大程度上是不同的。在代謝表型上，不同患者的差異表達(dá)基因有很大的差異，且平均每位患者有2.5個(gè)下調(diào)基因和1.5個(gè)上調(diào)基因是獨(dú)有的，且沒有任何一個(gè)基因在所有患者的差異表達(dá)基因中出現(xiàn)。在代謝網(wǎng)絡(luò)層面上，不同患者癌細(xì)胞發(fā)生改變的代謝系統(tǒng)規(guī)模也不盡相同，這反映了發(fā)展癌癥個(gè)性化治療的必要性。

本文的實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)了若干個(gè)性化的靶點(diǎn)基因，這些基因中的一部分已經(jīng)被證實(shí)與癌癥相關(guān)，另一部分則為生物醫(yī)學(xué)研究提供了新的線索。但由于目前的實(shí)驗(yàn)分析都停留在理論層面，缺少相應(yīng)的臨床醫(yī)學(xué)實(shí)驗(yàn)驗(yàn)證，因此本研究仍然具有一定的局限性，得到的結(jié)果對(duì)于臨床腫瘤個(gè)性化治療的應(yīng)用價(jià)值還需要進(jìn)一步的研究。在未來的工作中，課題組將嘗試擴(kuò)展可用的樣本數(shù)量和覆蓋的癌癥類型，以期不斷提高實(shí)驗(yàn)結(jié)果的準(zhǔn)確性和可信度。