郭鵬飛,賀平安

(浙江理工大學(xué)理學(xué)院，杭州 310018)

?

乳腺癌癌癥干細(xì)胞的特異基因識別

郭鵬飛,賀平安

(浙江理工大學(xué)理學(xué)院，杭州 310018)

乳腺癌是一種嚴(yán)重威脅女性健康的惡性腫瘤，癌癥干細(xì)胞假說的提出為乳腺癌的起因以及治療提供了新的模型。對746個乳腺癌樣本中的18409個基因和1035個miRNAs，通過生物信息學(xué)方法構(gòu)建共表達(dá)網(wǎng)絡(luò)，將其劃分到不同的共表達(dá)模塊中；利用胚胎干細(xì)胞和間充質(zhì)干細(xì)胞的特性進(jìn)一步篩選模塊，得到兩個大小分別為2019和859且與上述兩類干細(xì)胞相關(guān)的基因集；最后通過構(gòu)建這兩個基因集的調(diào)控網(wǎng)絡(luò)，篩選出兩個胚胎干細(xì)胞的特異性關(guān)鍵基因TPX2和MCM10，以及間充質(zhì)干細(xì)胞的特異基因COL5A2。這些基因可以作為癌癥干細(xì)胞的候選特異性標(biāo)志物，有望成為潛在的乳腺癌治療標(biāo)靶。

乳腺癌；胚胎干細(xì)胞；間充質(zhì)干細(xì)胞；基因調(diào)控網(wǎng)絡(luò)；關(guān)鍵基因

0 引 言

乳腺癌是全世界女性最常見的一種惡性腫瘤，其發(fā)病率占女性全身其他惡性腫瘤的10%左右。據(jù)統(tǒng)計，每年約有120萬新發(fā)乳腺癌病例，乳腺癌已成為女性發(fā)病率最高的癌癥[1]。隨著醫(yī)學(xué)的發(fā)展和治療手段的進(jìn)步，乳腺癌患者的生存率已經(jīng)得到很大改善，但它的耐藥性以及預(yù)后復(fù)發(fā)等問題依然困擾著許多研究者。

癌癥干細(xì)胞假說認(rèn)為癌癥是由一小群癌癥干細(xì)胞造成的，即癌癥干細(xì)胞是癌癥異常增殖、侵襲、轉(zhuǎn)移、耐藥以及復(fù)發(fā)等的根源。Wicha等[2]利用不同的細(xì)胞表面標(biāo)記物對細(xì)胞進(jìn)行標(biāo)記，驗(yàn)證了乳腺癌中存在著癌癥干細(xì)胞，這為癌癥干細(xì)胞的研究提供了堅實(shí)基礎(chǔ)。此后，Takebe等[3]研究發(fā)現(xiàn)Notch、Hedgehog (HH) 和Wnt等細(xì)胞通路在癌癥干細(xì)胞的致癌方面起著至關(guān)重要的作用，并對通過抑制這些通路從而控制干細(xì)胞復(fù)制、存活和分化這一新治療策略進(jìn)行了研究。

此外，許多研究者從基因?qū)用嫜芯堪┌Y干細(xì)胞涉及到的相關(guān)生物過程。Xu等[4]研究發(fā)現(xiàn)miR-214在調(diào)控卵巢癌干細(xì)胞性質(zhì)方面起著至關(guān)重要的作用，并且miR-214可作為治療卵巢癌的潛在的治療靶標(biāo)。Li等[5]研究發(fā)現(xiàn)長鏈非編碼RNA(long noncoding RNAs, lncRNAs)能夠抑制其標(biāo)靶mRNA在膠質(zhì)母細(xì)胞瘤干細(xì)胞(glioblastoma stem cells, GSCs)中分化，并利用這種方法識別出一些可用來治愈GSCs的候補(bǔ)lncRNAs。此外，Kalamohan等[6]應(yīng)用加權(quán)的基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)和富集度分析的方法分析胃癌的mRNA數(shù)據(jù)，發(fā)現(xiàn)胃癌的兩種亞型分別與不同的干細(xì)胞特征相關(guān)。

以上研究表明，癌癥干細(xì)胞假說能更好地解釋癌癥的起源，并指導(dǎo)癌癥治療。本文基于以上研究的成果，利用多種生物信息學(xué)方法分析現(xiàn)有數(shù)據(jù)庫中乳腺癌的基因表達(dá)譜數(shù)據(jù)，期望可以辨別出某些與癌癥干細(xì)胞相關(guān)的關(guān)鍵基因。該研究結(jié)果有助于人們更好地理解乳腺癌的發(fā)生、發(fā)展機(jī)制。

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

1.1 數(shù)據(jù)及數(shù)據(jù)預(yù)處理

本文研究的乳腺癌的miRNASeq、RNASeqV2數(shù)據(jù)以及乳腺癌患者的臨床數(shù)據(jù)，均下載于The Cancer Genome Atlas (TCGA)數(shù)據(jù)庫，其中RNASeqV2和miRNASeq數(shù)據(jù)是分別通過RNA測序和miRNA測序得到的樣本的基因表達(dá)和miRNA表達(dá)數(shù)據(jù)。

由于某些miRNAs可以通過調(diào)控基因表達(dá)進(jìn)而參與調(diào)控細(xì)胞分化和癌癥生成等重要生命過程，故本文合并兩種數(shù)據(jù)用于研究分析。首先，如果某個基因或miRNA在超過50%的樣本中其原始數(shù)據(jù)缺失，則將該基因或miRNA刪除；然后使用LIMMA package[7]將處理過的基因和miRNAs原始數(shù)據(jù)轉(zhuǎn)化為基因表達(dá)值，并將二者合并；最后使用反分位數(shù)歸一化(inversely normalize)的方法處理合并的基因表達(dá)數(shù)據(jù)，使其處于同一水平進(jìn)行后續(xù)分析。經(jīng)過預(yù)處理后得到一個由746個乳腺癌樣本中的18409個基因和1035個miRNAs的表達(dá)值構(gòu)成的數(shù)據(jù)集D=(dij)19444×746，其中dij表示第i個基因在第j個樣本中的表達(dá)值。

另外，從The Gene Expression Omnibus (GEO)數(shù)據(jù)庫下載了胚胎干細(xì)胞(GSE29625)和間充質(zhì)干細(xì)胞(GSE28974)的mRNA表達(dá)譜數(shù)據(jù)。

1.2 構(gòu)建基因共表達(dá)網(wǎng)絡(luò)

WGCNA算法[8]是一種從表達(dá)譜數(shù)據(jù)中挖掘模塊(module)信息的算法。在該算法中模塊被定義為一組具有類似表達(dá)譜的基因，即如果某些基因在一個生理過程或不同組織中總是具有相類似的表達(dá)變化，則將其定義為一個模塊。

(1)

(2)

(3)

另外考慮到基因i可以通過基因μ與基因j相互作用，故將鄰接矩陣被轉(zhuǎn)換成拓?fù)渚仃嚘?(wij)，其中：

(4)

這里lij=∑uaiuauj，表示與基因i、j都相鄰的基因μ之間的鄰接系數(shù)乘積和；ki=∑uaiu為基因i單獨(dú)連接的節(jié)點(diǎn)的鄰接系數(shù)的和。

理論上，胚胎干細(xì)胞和間充質(zhì)干細(xì)胞在原發(fā)性腫瘤，以及癌細(xì)胞的循環(huán)過程和轉(zhuǎn)移器官中都可以檢測到。因此本文選擇這兩種干細(xì)胞為代表通過基因集的富集度分析[9]，查找與癌癥干細(xì)胞相關(guān)的基因。為此，從MSigDB(molecular signatures database)數(shù)據(jù)庫[10]下載了15個與胚胎干細(xì)胞特性相關(guān)的基因集以及10個跟間充質(zhì)干細(xì)胞特征相關(guān)的基因集，作為背景基因進(jìn)行富集度分析。

基因集的富集度分析基于超幾何分布，服從超幾何分布(k-1,K,N-K,n)的概率p可通過式(5)來計算：

(5)

其中：n表示模塊中基因的個數(shù)；K表示與胚胎干細(xì)胞特性相關(guān)的基因集或者跟間充質(zhì)干細(xì)胞特征相關(guān)的基因集中基因的個數(shù)；k和N分別為上述模塊和基因集的交集和并集中基因的個數(shù)。錯誤發(fā)現(xiàn)率[11](false discovery rate,FDR)用于評價基因模塊是否富集于與兩種干細(xì)胞相關(guān)的基因集。

1.4 優(yōu)化與癌癥干細(xì)胞相關(guān)的基因集

通過基因集的富集度分析得到兩個分別與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集，為了優(yōu)化這兩個基因集，本文應(yīng)用多尺度嵌合基因共表達(dá)網(wǎng)絡(luò)分析(multiscale embedded gene co-expression network analysis, MEGENA)算法[12]對二者重新進(jìn)行精確分類。

MEGENA算法首先計算任意兩個基因之間的相關(guān)性，并依據(jù)相關(guān)性的大小對基因?qū)ε判?；接著通過平面最大過濾圖算法(planar maximally filtered graph , PMFG)將其嵌入拓?fù)渚W(wǎng)絡(luò)，從而構(gòu)建平面濾波網(wǎng)絡(luò)(planar filtered networks, PFNs)；然后通過最短路徑距離、本地路徑索引和整體模塊性三個標(biāo)準(zhǔn)對最初的PFNs進(jìn)行多次迭代處理，得到更精確的分類。

1.5 基因表達(dá)的調(diào)控網(wǎng)絡(luò)的建立

隨機(jī)森林算法[13]是一種基于決策樹模型的算法，它主要通過一個重要性評分矩陣來推斷調(diào)控網(wǎng)絡(luò)。相比于其他構(gòu)建調(diào)控網(wǎng)絡(luò)的方法，隨機(jī)森林算法可以得到一個有向的調(diào)控網(wǎng)絡(luò)，使得基因間的調(diào)控關(guān)系更加明確。故本文利用隨機(jī)森林算法在與癌癥干細(xì)胞相關(guān)的基因集中構(gòu)建基因調(diào)控網(wǎng)絡(luò)。

隨機(jī)森林算法將預(yù)測n個基因間的調(diào)控網(wǎng)絡(luò)的問題轉(zhuǎn)化為求解n個不同的回歸問題。首先選取一個基因作為靶基因(因變量)，其余n-1個基因作為輸入基因(自變量)，做回歸分析預(yù)測靶基因。每個輸入基因在預(yù)測靶基因過程中計算相應(yīng)的變量重要性評分(variable importance measure，VIM)，并以此作為推定基因間調(diào)控關(guān)系的指標(biāo)。將得到的所有靶基因與輸入基因之間的調(diào)控關(guān)系依據(jù)其大小排序，從而構(gòu)造調(diào)控網(wǎng)絡(luò)。本文用R語言中的randomForest package[14]構(gòu)建有向的基因調(diào)控網(wǎng)絡(luò)，同時用Cytoscape軟件[15]實(shí)現(xiàn)基因調(diào)控網(wǎng)絡(luò)的可視化。

在有向的基因調(diào)控網(wǎng)絡(luò)中，基因的頂點(diǎn)出度是以該點(diǎn)為起點(diǎn)的邊的個數(shù)。本文根據(jù)基因的頂點(diǎn)出度的大小篩選與癌癥干細(xì)胞相關(guān)的關(guān)鍵基因。

1.6Kaplan-Meier生存分析

生存分析是將事件的結(jié)果和出現(xiàn)這一結(jié)果所經(jīng)歷的時間結(jié)合起來分析的一種統(tǒng)計分析方法[16]。Kaplan-Meier生存分析將乘積極限法應(yīng)用于臨床數(shù)據(jù)中樣本生存或死亡這兩種狀態(tài)所對應(yīng)的生存時間，從而計算出樣本的生存率及其標(biāo)準(zhǔn)誤差。然后利用log-rank檢驗(yàn)來比較兩組或多組生存率，并通過p-value來評價不同組的生存率是否相同。

對得到的關(guān)鍵基因，本文通過構(gòu)建Kaplan-Meier生存曲線來驗(yàn)證它們對乳腺癌的重要性。

2 結(jié) 果

2.1 表達(dá)數(shù)據(jù)的聚類分析

利用WGCNA算法，輸入數(shù)據(jù)集D中的數(shù)據(jù)，首先計算任意兩個基因之間的皮爾森相關(guān)性系數(shù)得到相關(guān)性矩陣。接著通過無尺度網(wǎng)絡(luò)原則確定β。如圖1所示，當(dāng)β=5時R2= 0.8220，因此本文選擇β=5作為加權(quán)系數(shù)將相關(guān)性矩陣轉(zhuǎn)化為鄰接矩陣。最后利用節(jié)點(diǎn)的相異程度進(jìn)行分層聚類，結(jié)果746個乳腺癌樣本中的18409個基因和1035個miRNAs被聚類到47個不同的模塊中。47個模塊分別被記作M1—M47，并且每個模塊的大小從33到2598數(shù)目不等。由于這些共表達(dá)基因傾向于功能相關(guān)的，故這種聚類方式也意味著乳腺癌的轉(zhuǎn)錄組包括47個不同或相關(guān)的生物過程。而這些生物過程有助于研究乳腺癌中的分子機(jī)制和關(guān)鍵性驅(qū)動因子，因此這些模塊值得深入研究[6]。

圖1 加權(quán)系數(shù) β的選取

2.2 識別與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的模塊

將WGCNA算法得到的模塊和MSigDB數(shù)據(jù)庫中下載的與兩種干細(xì)胞特性相關(guān)的基因集作為輸入數(shù)據(jù)進(jìn)行富集度分析，本文以FDR<0.05為標(biāo)準(zhǔn)來確定結(jié)果。

在與胚胎干細(xì)胞特性相關(guān)的基因集的富集度分析結(jié)果中，有8個基因模塊與胚胎干細(xì)胞特性相關(guān)，它們分別是M1—M8，具體結(jié)果見表1。另一方面，有6個模塊富集于間充質(zhì)干細(xì)胞，分別為M2、M5、M6、M8、M9和M10，其結(jié)果如表2所示。

表1 與胚胎干細(xì)胞性質(zhì)相關(guān)的模塊的富集度分析結(jié)果

表2 與間充質(zhì)干細(xì)胞特征相關(guān)的模塊的富集度分析結(jié)果

在表1和表2中，第一列和第二列分別表示模塊及其大小，第四列是每個模塊與從MSigDB數(shù)據(jù)庫中得到的任一基因集進(jìn)行一次富集度分析得到的FDR值，第三列是根據(jù)FDR的值從小到大的排序，第五列是每個基因模塊中富集于MSigDB數(shù)據(jù)庫中的基因集。

為了進(jìn)一步分析這些基因，本文將與胚胎干細(xì)胞性質(zhì)相關(guān)的8個模塊合并成一個基因集E1，內(nèi)含5751個基因；同時將富集于間充質(zhì)干細(xì)胞的6個模塊合并成一個基因集E2，內(nèi)含3441個基因。

2.3 分別確定與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集

在WGCNA算法的結(jié)果中，有4個模塊M2、M5、M6、M8同時富集于胚胎干細(xì)胞和間充質(zhì)干細(xì)胞，這使得基因集E1和E2中包含許多相同的基因。為了優(yōu)化上述兩個基因集，本文對這些基因模塊作如下處理：

首先，合并兩個基因集E1和E2得到新的基因集E3，它包含6102個基因。對于E3中的基因，根據(jù)由乳腺癌樣本中基因和miRNAs的表達(dá)值構(gòu)成的數(shù)據(jù)集D構(gòu)造它的一個子矩陣D1=(dij)6102×746。使用MEGENA算法對該數(shù)據(jù)重新分類得到新的模塊，并對新的模塊進(jìn)行基因集富集度分析，重新篩選出與兩種干細(xì)胞相關(guān)的模塊，最后合并模塊得到新的與兩種干細(xì)胞相關(guān)的基因集。在這一過程中，本文得到兩個基因個數(shù)分別為2009和572的且與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集F1和F2。

其次，對基因集E1中的基因，重復(fù)上述過程，得到分別由1824和802個基因構(gòu)成的與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集F3和F4，并且F3和F4無交集。

然后，對基因集E2中的基因，重復(fù)上述過程，得到兩個分別與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的無交集基因集F5和F6，其大小分別為57和386。

最后，取F1、F3和F5的并集，得到一個包含2019個與胚胎干細(xì)胞相關(guān)的基因集G1；取F2、F4和F6的并集，得到一個大小為859的與間充質(zhì)干細(xì)胞相關(guān)的基因集G2。而且新得到的基因集G1和G2沒有交集。表3—表4為上述三組數(shù)據(jù)利用MEGENA算法分類后，新的模塊進(jìn)行富集度分析的結(jié)果。

表3 MEGENA算法分類與胚胎干細(xì)胞性質(zhì)相關(guān)的模塊的富集度分析結(jié)果

表4 MEGENA算法分類與間充質(zhì)干細(xì)胞特征相關(guān)的模塊的富集度分析結(jié)果

2.4 驗(yàn)證基因集G1和G2

為了進(jìn)一步驗(yàn)證上述過程得到的兩個與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集。首先合并從GEO數(shù)據(jù)庫下載的癌癥胚胎干細(xì)胞(GSE29625)和間充質(zhì)干細(xì)胞(GSE28974)的mRNA表達(dá)譜數(shù)據(jù)，然后分別使用分位數(shù)歸一化[17]和反分位數(shù)歸一化的方法處理合并后的數(shù)據(jù)，最終得到一個由10195個基因構(gòu)成的基因表達(dá)數(shù)據(jù)T1=(tij)10195×24，這里tij為第i個基因在第j個樣本中的表達(dá)值。T1用于驗(yàn)證上述過程得到的兩個與胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的基因集G1和G2。

圖2(a)是由與胚胎干細(xì)胞相關(guān)的基因集G1在數(shù)據(jù)T1中的表達(dá)值構(gòu)成的熱圖，圖中每個小方格表示一個基因在樣本中的表達(dá)量，顏色表示表達(dá)量的大小。其中首字母為E的是GSE29625中的樣本，首字母為M的是GSE28974中的樣本。圖中結(jié)果表明該基因集的大部分基因在整合數(shù)據(jù)中的GSE29625樣本中具有顯著的高表達(dá)，在GSE28974樣本中具有顯著的低表達(dá)。類似地，圖2(b)是由與間充質(zhì)干細(xì)胞相關(guān)的基因集在整合數(shù)據(jù)T1中的表達(dá)值構(gòu)成的熱圖。該基因集中的大部分基因在整合數(shù)據(jù)中的GSE28974樣本中具有顯著的高表達(dá)。

圖2 與胚胎干細(xì)胞相關(guān)的基因集和與間充質(zhì)干細(xì)胞的基因集在兩種干細(xì)胞的整合數(shù)據(jù)中的表達(dá)模式

基因集G1和G2在整合數(shù)據(jù)T1中的不同表達(dá)模式表明二者與兩類干細(xì)胞具有明顯的相關(guān)性，進(jìn)一步說明使用本文方法得到的結(jié)果具有較強(qiáng)的可靠性。

2.5 構(gòu)建基因調(diào)控網(wǎng)絡(luò)

對與胚胎干細(xì)胞相關(guān)的基因集G1中的2019個基因，利用其在數(shù)據(jù)集D中的表達(dá)值，使用R語言中的randomForest package[14]構(gòu)建它們的調(diào)控網(wǎng)絡(luò)。隨機(jī)森林算法中最重要的參數(shù)有兩個：一個是建立決策樹的個數(shù)，本文取1000；另一個是每個節(jié)點(diǎn)可選擇的候選輸入基因個數(shù)，在本文中該參數(shù)取全部輸入基因數(shù)的平方根。

利用多重假設(shè)檢驗(yàn)，求出重要性評分矩陣的每一個值的FDR值，首先取FDR<0.01的邊來控制調(diào)控網(wǎng)絡(luò)的大小，得到一個包含53298條邊的調(diào)控網(wǎng)絡(luò)。圖3是由此網(wǎng)絡(luò)中全部基因之間的VIM值繪制的直方圖，從圖中可以看出大部分基因間的重要性評分值小于0.01。為了進(jìn)一步控制調(diào)控網(wǎng)絡(luò)的規(guī)模，本文僅選取基因之間的重要性評分值大于0.01的邊。最后得到了一個含有15720條邊的有向調(diào)控網(wǎng)絡(luò)。

圖3 53298個基因?qū)Φ?VIM值的頻率分布直方圖

在構(gòu)造的有向調(diào)控網(wǎng)絡(luò)中，如果一個基因同時調(diào)控多個基因，那么它肯定在某個生物過程中起重要作用，所以文本重點(diǎn)關(guān)注那些處于調(diào)控關(guān)系上游的基因。本文使用R語言中的igraph package[18]計算網(wǎng)絡(luò)中每個頂點(diǎn)的出度，并根據(jù)其大小進(jìn)行排序。在這個有向調(diào)控網(wǎng)絡(luò)中，TPX2、MCM10、CEP55、BUB1、NCAPG、NCAPH和BUB1B等基因有較大的頂點(diǎn)出度，具體結(jié)果如表5所示。特別是TPX2和MCM10這兩個基因調(diào)控下游基因的個數(shù)都超過110，所以本文認(rèn)為這兩個基因是與胚胎干細(xì)胞相關(guān)的關(guān)鍵基因。圖4(a)—(b)是以基因TPX2和MCM10為核心，以及它們調(diào)控的基因之間的調(diào)控關(guān)系構(gòu)成的調(diào)控子網(wǎng)絡(luò)。為了顯示清晰，該調(diào)控子網(wǎng)絡(luò)中僅畫出了VIM值大于0.02的邊。

表5 基因集G1調(diào)控網(wǎng)絡(luò)中Top12基因的頂點(diǎn)出度

對于包含859個與間充質(zhì)干細(xì)胞相關(guān)的基因集G2，本文重復(fù)同樣的過程。結(jié)果發(fā)現(xiàn)，在與間充質(zhì)干細(xì)胞相關(guān)的基因集的調(diào)控網(wǎng)絡(luò)中，COL5A2、FBN1和COL1A2等基因具有較大的頂點(diǎn)出度，具體結(jié)果如表6所示。其中COL5A2基因調(diào)控的基因數(shù)超過100。因此，基因COL5A2被看作是與間充質(zhì)干細(xì)胞相關(guān)的關(guān)鍵基因。它的調(diào)控子網(wǎng)絡(luò)見圖4(c)。

(a)TPX2

(b)MCM10

(c)COL5A2圖4 三個關(guān)鍵基因及其下游調(diào)控基因構(gòu)成的調(diào)控子網(wǎng)絡(luò)

基因頂點(diǎn)出度基因頂點(diǎn)出度COL5A2102COL6A371FBN187VCAN67COL1A287LUM67COL3A174BNC266THBS272CDH1165

由于上述3個基因在相應(yīng)的調(diào)控網(wǎng)絡(luò)中具有極高的頂點(diǎn)出度，說明TPX2、MCM10和COL5A2在癌癥的胚胎干細(xì)胞和間充質(zhì)干細(xì)胞的自我更新、分化過程中具有重要作用。故上述三個基因可作為辨別乳腺癌的癌癥干細(xì)胞的特征基因，以及治愈乳腺癌的潛在的生物靶基因。

2.6 關(guān)鍵基因的生物學(xué)分析

本文將3個關(guān)鍵基因的表達(dá)譜數(shù)據(jù)與TCGA數(shù)據(jù)庫中的臨床數(shù)據(jù)整合成一個由657個樣本構(gòu)成的新數(shù)據(jù)T2進(jìn)行Kaplan-Meier生存分析，從而進(jìn)一步研究關(guān)鍵基因的表達(dá)方式對乳腺癌患者的生存率的影響，結(jié)果如圖5所示。在圖5中，對應(yīng)的曲線分別為在基因TPX2和MCM10的高表達(dá)和低表達(dá)情況下乳腺癌患者的生存曲線，其中x軸表示乳腺癌患者的生存時間，y軸表示患者的生存率；event=1代表患者死亡。觀察圖5發(fā)現(xiàn)處于TPX2和MCM10高表達(dá)組的癌癥患者相比于低表達(dá)組的患者有明顯高的死亡率。而且假設(shè)檢驗(yàn)的p-value都小于0.05，也表明關(guān)鍵基因的不同表達(dá)方式對乳腺癌患者生存率的影響顯著不同。

基因COL5A2沒有上述結(jié)論，但Weng等[19]通過研究血小板反應(yīng)蛋白2(thrombospondin2,THBS2)的表達(dá)模式在肺癌發(fā)展中的作用，發(fā)現(xiàn)COL5A2基因作為THBS2的一個共表達(dá)基因，它們的高表達(dá)使得肺癌患者具有較低存活率。此外，F(xiàn)ischer等[20]通過對比膠原蛋白的基因在結(jié)腸直腸癌患者和正常結(jié)腸上皮的組織樣品中的差異表達(dá)，發(fā)現(xiàn)基因COL5A2在基質(zhì)中的表達(dá)與結(jié)腸直腸癌相關(guān)。Zhang等[21]利用TCGA數(shù)據(jù)庫中的卵巢癌數(shù)據(jù)構(gòu)建了貝葉斯網(wǎng)絡(luò)，其中基因COL5A2同樣被發(fā)現(xiàn)是關(guān)鍵基因。

3 結(jié) 論

癌癥干細(xì)胞假說認(rèn)為癌癥很可能起源于干細(xì)胞的非正常分化，那么通過中斷癌癥干細(xì)胞的自我更新從而造成其自我更新障礙，完全可以成為一種治療癌癥的理想方式。本文基于這一思想，利用生物信息學(xué)中WGCNA算法和MEGENA算法分析乳腺癌的基因和miRNA混合表達(dá)數(shù)據(jù)，將其劃分為具有不同或相似生物功能的基因共表達(dá)模塊，同時利用MSigDB數(shù)據(jù)庫中與胚胎干細(xì)胞特性和間充質(zhì)干細(xì)胞特征相關(guān)的基因集進(jìn)行富集度分析，得到了兩個分別與乳腺癌胚胎干細(xì)胞和間充質(zhì)干細(xì)胞相關(guān)的且由2019個基因和859個基因組成的基因集。另外從GEO數(shù)據(jù)庫中下載了癌癥胚胎干細(xì)胞和間充質(zhì)干細(xì)胞的mRNA表達(dá)譜數(shù)據(jù)，并通過查看上述兩個基因集在mRNA表達(dá)譜數(shù)據(jù)中的表達(dá)模式，驗(yàn)證了本文得到的兩個基因集是可靠的。

進(jìn)一步利用隨機(jī)森林算法在上述兩個基因集中構(gòu)建有向的調(diào)控網(wǎng)絡(luò)。通過調(diào)控網(wǎng)路的分析，發(fā)現(xiàn)了3個在癌癥的胚胎干細(xì)胞和間充質(zhì)干細(xì)胞自我更新、分化過程中起重要作用的關(guān)鍵基因：TPX2、MCM10和COL5A2。生存分析和已有的結(jié)果進(jìn)一步說明這三個基因是與乳腺癌密切相關(guān)的，可以作為治療乳腺癌的潛在的治療靶點(diǎn)。

[1] BENSON J R, JATOI I, KEISEH M, et al. Early breast cancer [J]. Lancet,2009,373(9673):1463-1479

[2] WICHA M S, DONTU G, AL-HAJJ M, et al. Stem cells in normal breast development and breast cancer[J].Breast Cancer Research,2003,36(1):59-72.

[3] TAKEBE N, MIELE L, HARRIS P J, et al. Targeting Notch, Hedgehog, and Wnt pathways in cancer stem cells: clinical update [J].Nature Reviews Clinical Oncology,2015,12(8):445-464.

[4] XU C X, XU M, TAN L, et al. MicroRNA miR-214 regulates ovarian cancer cell stemness by targeting p53/Nanog [J].Journal of Biological Chemistry,2012,287(42):34970-34978.

[5] LI H. Differential long non-coding RNA and mRNA expression in differentiated human glioblastoma stem cells [J].Mol Med Rep,2016,14(3):2067-2076.

[6] KALAMOHAN K, PERIASAMY J, BHASKAR R D, et al. Transcriptional co-expression network reveals the involvement of varying stem cell features with different dysregulations in different gastric cancer subtypes [J].Molecular Oncology,2014,8(7):1306-1325.

[7] LAW C W, CHEN Y, SHI W, et al. Voom: precision weights unlock linear model analysis tools for RNA-seq read counts [J].Genome Biology,2014,15(2):1417.

[8] 宋長新,雷萍,王婷.基于WGCNA算法的基因共表達(dá)網(wǎng)絡(luò)構(gòu)建理論及其R軟件實(shí)現(xiàn)[J].基因組學(xué)與應(yīng)用生物學(xué),2013,32(1):135-141.

[9] RIVALS I, PERSONNAZ L, TAING L, et al. Enrichment or depletion of a GO category within a class of genes: which test ?[J].Bioinformatics,2007,23(4):401-407.

[10] LIBERZON A, SUBRAMANIAN A, PINCHBACK R, et al. Molecular signatures database (MSigDB) 3.0[J].Bioinformatics,2011,27(12):1739-1740.

[11] BENJAMINI Y. Discovering the false discovery rate[J].Journal of the Royal Statistical Society,2010,72(4):405-416.

[12] SONG W M, ZHANG B. Multiscale embedded gene co-expression network analysis [J].Plos Computational Biology,2015,11(11),e1004574.

[13] 侯艷,楊凱,李康.基于隨機(jī)森林回歸的網(wǎng)絡(luò)構(gòu)建方法及應(yīng)用[J].中國衛(wèi)生統(tǒng)計,2015,32(4):558-561.

[14] LIAW A, WIENER M. Classification and regression by randomforest [J].R News,2002,2(3):18-22.

[15] SHANNON P, MARKIEL A, OZIER O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks [J].Genome Res,2003,13(11):2498-2504.

[16] 孫振球. 醫(yī)學(xué)統(tǒng)計學(xué)[M].3版. 北京：人民衛(wèi)生出版社,2014:306-313.

[17] BOLSTAD B M. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias[J].Bioinformatics,2003,19(2):185-193.

[18] CSARDI G, NEPUSZ T. The igraph software package for complex network research[J].InterJournal Complex Systems,2006,1695(5):1-9.

[19] WENG T Y, WANG C Y, HUNG Y H, et al.Differential expression pattern of THBS1 and THBS2 in Lung Cancer: clinical outcome and a systematic-analysis of Microarray Databases [J].Plos One,2016,11(8):e0161007.

[20] FISCHER H, STENLING R, RUBIO C, et al. Colorectal carcinogenesis is associated with stromal expression of COL11A1 and COL5A2 [J].Carcinogenesis,2001,22(6):875-878.

[21] ZHANG Q, BURDETTE J E, WANG J P. Integrative network analysis of TCGA data for ovarian cancer[J].BMC Systems Biology,2014,8(1):1-18.

(責(zé)任編輯： 康 鋒)

Identification of Specific Genes of Cancer Stem Cells of Breast Cancer

GUOPengfei,HEPingan

(School of Sciences, Zhejiang Sci-Tech University, Hangzhou 310018, China)

Breast cancer is a kind of malignant tumor which seriously threats the health of global female. However, the hypothesis of cancer stem cell (CSC) provides a new model for breast cancer causes and treatment. In the paper, coexpression network was constructed with the bioinformatics method for 18409 genes and 1035 miRNA in 746 breast cancer samples, and they were divided into different coexpression modules. The characteristics of embryonic stem cells and mesenchymal stem cells were utilized to further screen the modules, and two gene sets related to the above two types of stem cells (size: 2019 and 859) were gained respectively. Finally, regulatory network for the two gene sets were constructed to screen specific hub genesTPX2 andMCM10 of two embryonic stem cells as well as specific geneCOL5A2 of mesenchymal stem cells. These genes can be considered as candidate specific biomarkers of CSC and potential therapeutic targets in the treatment of breast cancer.

breast cancer; embryonic stem cells; mesenchymal stem cells; gene regulatory network; hub genes

10.3969/j.issn.1673-3851.2017.05.023

2016-11-22 網(wǎng)絡(luò)出版日期： 2017-03-28

國家自然科學(xué)基金項(xiàng)目(61170110, 61272312)；浙江省自然科學(xué)基金項(xiàng)目(LY14F020049)

郭鵬飛(1990-)，男，山西忻州人，碩士研究生，主要從事生物信息學(xué)方面的研究。

賀平安，E-mail：pinganhe@zstu.edu.cn

Q612

A

1673- 3851 (2017) 03- 0451- 10