癌癥相關(guān)基因選擇性剪接進(jìn)化數(shù)據(jù)庫的構(gòu)建

徐敏杰，竇同海，徐佳熹，江建平,，劉 匆，付茂賓，高 原，陳誠文，張 亮，周 雁,

(1. 復(fù)旦大學(xué) 生命科學(xué)學(xué)院 遺傳工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200438；2. 國家人類基因組南方研究中心，上海市疾病與健康基因組學(xué)重點(diǎn)實(shí)驗(yàn)室，上海 201203 )

?

癌癥相關(guān)基因選擇性剪接進(jìn)化數(shù)據(jù)庫的構(gòu)建

徐敏杰1,2,*，竇同海1,*，徐佳熹1，江建平1,2，劉 匆1，付茂賓1，高 原1，陳誠文2，張 亮2，周 雁1,2

(1. 復(fù)旦大學(xué) 生命科學(xué)學(xué)院 遺傳工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200438；2. 國家人類基因組南方研究中心，上海市疾病與健康基因組學(xué)重點(diǎn)實(shí)驗(yàn)室，上海 201203 )

選擇性剪接是真核生物基因調(diào)控的基本調(diào)節(jié)機(jī)制之一，與各種類型的生理和病理活動相關(guān).癌癥相關(guān)基因的不正常剪接可能與多種癌癥的發(fā)生發(fā)展有關(guān).作為選擇性剪接進(jìn)化的主體，選擇性剪接外顯子展示了其在不同物種中多樣的進(jìn)化功能.本文系統(tǒng)整理了2989個癌癥相關(guān)基因的各項(xiàng)功能，通過比較基因組學(xué)的分析方法總結(jié)了癌癥相關(guān)基因的選擇性剪接外顯子的功能和進(jìn)化關(guān)系，建立了癌癥相關(guān)基因選擇性剪接進(jìn)化數(shù)據(jù)庫(ASeeDB數(shù)據(jù)庫).ASeeDB包含癌癥相關(guān)基因外顯子區(qū)域的進(jìn)化保守性、結(jié)構(gòu)域預(yù)測、Ka/Ks值、以及基因、轉(zhuǎn)錄本、外顯子區(qū)域3個層次的表達(dá)量統(tǒng)計(jì)等信息，結(jié)合這些信息用戶可以方便的檢索有研究意義的基因或者外顯子.外顯子區(qū)域蛋白質(zhì)結(jié)構(gòu)域的預(yù)測可以幫助了解其可能的功能，而外顯子是否存在選擇性剪接又可以推及包含該外顯子的轉(zhuǎn)錄本是否具有相似的功能以及推測剪接是否會對蛋白質(zhì)功能發(fā)生影響.數(shù)據(jù)庫提供的物種間的序列比對可以幫助用戶發(fā)現(xiàn)沒有注釋的外顯子區(qū)域或者是保留的失去功能的外顯子.相比于基因?qū)用娴腒a/Ks，數(shù)據(jù)庫提供的外顯子層面的Ka/Ks對于發(fā)現(xiàn)適應(yīng)性進(jìn)化事件具有更高的敏感性，可以更加方便地預(yù)示基因中未知功能的區(qū)域.

癌癥; 選擇性剪接; 數(shù)據(jù)庫; 外顯子

選擇性剪接(alternative splicing)，又稱可變剪接，是一種存在于真核生物細(xì)胞內(nèi)的重要的基因表達(dá)調(diào)控機(jī)制.選擇性剪接是指同一個mRNA前體通過不同的剪接方式選擇不同的剪接位點(diǎn)進(jìn)而產(chǎn)生不同的成熟的mRNA(轉(zhuǎn)錄本isoform)的過程.選擇性剪接最初在腺病毒(adenovirus)研究中發(fā)現(xiàn)，之后在1981年，真核生物的降鈣素(calcitonin)基因中也觀察到了該現(xiàn)象[1-2].大量的研究表明，選擇性剪接與細(xì)胞內(nèi)的多種生理以及病理活動相關(guān)[3].選擇性剪接外顯子，即在轉(zhuǎn)錄過程中同一個基因由于選擇性剪接出現(xiàn)在不同轉(zhuǎn)錄本中的外顯子.通過選擇性剪接機(jī)制，外顯子的不同組合方式在蛋白質(zhì)多樣性和功能分化中發(fā)揮重要的作用.以外顯子作為進(jìn)化分析的主要對象，可以明顯地提高進(jìn)化選擇壓力的計(jì)算敏感度，對于縮小研究范圍具有重要意義.

長期以來的研究已經(jīng)表明，可變剪接對于信號轉(zhuǎn)導(dǎo)、細(xì)胞凋亡、免疫應(yīng)答等有著廣泛的影響[4-7]，其紊亂會導(dǎo)致多種疾病產(chǎn)生，如中樞神經(jīng)系統(tǒng)病變、脂代謝紊亂、阿爾茨海默氏癥、多種類型的腫瘤等一系列的疾病[8].在腫瘤中可變剪接對其發(fā)生發(fā)展的影響是多方面的.參與可變剪接的調(diào)控因子(如hnRNP，SR蛋白)時空表達(dá)的改變會造成腫瘤相關(guān)基因(如CD44、Ron、S6K1等)的不同剪接體異常表達(dá)，從而誘發(fā)肺癌、直腸癌、卵巢癌.另一方面，腫瘤基因本身突變(如APC、BRCA1的突變)，以及含有外顯子剪接增強(qiáng)子(ESE)的側(cè)翼序列突變(如NF1、NF2)也會影響其正常剪接位點(diǎn)的選擇，從而引發(fā)腫瘤.目前，研究人員已經(jīng)明確剪接突變體Survivin 2B具有促凋亡作用，其表達(dá)的下調(diào)會導(dǎo)致乳腺癌或轉(zhuǎn)移性胃癌的發(fā)生[8].這些例子說明在研究探尋腫瘤相關(guān)基因的功能時，有必要考慮不同轉(zhuǎn)錄本的功能差異，否則將會產(chǎn)生難以解釋或相互矛盾的結(jié)果.

近年來真核生物基因組與基因表達(dá)數(shù)據(jù)不斷增多，目前僅Ensembl數(shù)據(jù)庫(www.ensembl.org)中可供橫向比較的脊椎動物就超過43種，主流的可變剪接專業(yè)數(shù)據(jù)庫數(shù)據(jù)量也在呈幾何級數(shù)增長，例如ASTD(http：∥www.ebi.ac.uk/astd/relnotes.html)擁有人類、大鼠和小鼠的超過4萬條基因的可變剪接數(shù)據(jù)，這都為大規(guī)模從可變剪接角度分析基因功能提供了有用的素材.與此同時，系統(tǒng)性地對腫瘤相關(guān)基因進(jìn)行分析的想法也日益成熟，并已經(jīng)從Ras，Wnt等經(jīng)典基因家族擴(kuò)展到多條腫瘤基因通路(如KEGG腫瘤基因及242個通路)、基于文獻(xiàn)檢索的整合數(shù)據(jù)庫(如Cancer Gene數(shù)據(jù)庫，http：∥cbio.mskcc.org/Cancer Gene)和特定腫瘤相關(guān)基因數(shù)據(jù)庫(如http：∥www.megabionet.org/bio/hlung).2005年12月美國國立衛(wèi)生研究院(NIH)開始啟動腫瘤基因組計(jì)劃(The Cancer Genome Atlas)，相信隨著類似研究計(jì)劃的不斷推出和基因組信息的不斷完善，腫瘤相關(guān)基因的信息也將更加豐富.

選擇性剪接數(shù)據(jù)如ASTD由于主要基于一代測序的EST數(shù)據(jù)，并且已停止更新，現(xiàn)已難以滿足研究需求.最新的基于RNA-seq選擇性剪接的數(shù)據(jù)庫DBATE(http：∥bioinformatica.uniroma2.it/DBATE/)[9]充分利用了高通量測序的數(shù)據(jù)分析并展示了基因和轉(zhuǎn)錄本的表達(dá)量，但并未細(xì)化到外顯子的表達(dá)水平，且涵蓋的物種較少，并存在不能穩(wěn)定訪問的情況.由于很多癌癥相關(guān)基因的不正確地剪接已經(jīng)被證實(shí)與癌癥的發(fā)生和發(fā)展相關(guān)，而目前還未有專門針對癌癥相關(guān)基因的選擇性剪接數(shù)據(jù)庫，因此建立一個專門研究癌癥相關(guān)基因的選擇性剪接數(shù)據(jù)庫對于癌癥基因的研究具有重要作用.

本文試圖從選擇性剪接的角度研究癌癥相關(guān)基因，從基因、轉(zhuǎn)錄本以及外顯子3個表達(dá)水平結(jié)合物種間的進(jìn)化關(guān)系，通過數(shù)據(jù)整理和整合，建立友好的訪問形式，為癌癥相關(guān)基因的研究提供一個良好的平臺和新的研究思路.

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

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

當(dāng)前發(fā)布的數(shù)據(jù)庫(Alternative Splicing Exon Evolution Database of cancer relation genes, ASeeDB)主要數(shù)據(jù)來源基于Ensembl 68版本[10]，并輔以NCBI數(shù)據(jù)庫中的數(shù)據(jù)作為補(bǔ)充.現(xiàn)數(shù)據(jù)庫覆蓋了人類、黑猩猩、小鼠、大鼠、鴨嘴獸和斑馬魚6個物種，并提供了這6個物種相關(guān)基因基本信息，包括染色體位置，核酸序列，蛋白質(zhì)序列和轉(zhuǎn)錄本信息等.數(shù)據(jù)庫主要由7個模塊組成，分別為As classifier、PAML、Exon profiler、ESE Finder、Repeatmasker、InterproScan和 exon expression(圖1).通過這7個模塊的分析與統(tǒng)計(jì)，數(shù)據(jù)庫共覆蓋了來自Memorial Sloan-Kettering Cancer Center 癌癥相關(guān)基因數(shù)據(jù)庫中2989個癌癥相關(guān)基因，提供了這些癌癥相關(guān)基因的進(jìn)化關(guān)系以及表達(dá)量等信息.

1.2 方法

1.2.1 數(shù)據(jù)收集

通過癌癥相關(guān)數(shù)據(jù)庫的信息共獲取到2989個癌癥相關(guān)基因作為研究對象.以歐洲生物信息學(xué)中心數(shù)據(jù)庫Ensembl 68版本作為數(shù)據(jù)庫的主要基礎(chǔ)數(shù)據(jù)來源，下載2989個癌癥相關(guān)基因的染色體位置，核酸序列，氨基酸序列以及人類、黑猩猩、小鼠、大鼠、鴨嘴獸和斑馬魚的直系同源關(guān)系.NMD(Nonsense-mediated mRNA Decay)和不能翻譯成蛋白質(zhì)的轉(zhuǎn)錄本從數(shù)據(jù)中過濾去除.通過整理統(tǒng)計(jì)將這些信息儲存于一個臨時數(shù)據(jù)庫待后續(xù)分析.由于Ensembl 數(shù)據(jù)庫中只包含位于primary assembly基因組上的基因，所以NCBI上的數(shù)據(jù)同時用來補(bǔ)充到數(shù)據(jù)庫中以保證ASeeDB數(shù)據(jù)庫數(shù)據(jù)的完整性.

1.2.2 外顯子區(qū)域分析和注釋

由于同一個基因中的外顯子可能會因?yàn)檫x擇性剪接而導(dǎo)致其在不同的轉(zhuǎn)錄本中具有不同的長度，這給建立一一對應(yīng)的直系同源關(guān)系帶來了巨大的難度.為了解決這個問題，我們定義了“外顯子區(qū)域”，通過將在不同轉(zhuǎn)錄本中有相互重疊的外顯子合并成一個理論上的“外顯子區(qū)域”(圖2，見第652頁).每個外顯子區(qū)域在染色體上的起始和終止位置被記錄下來，同時記錄每個外顯子區(qū)域包含的真實(shí)外顯子.通過這種方式就可以大大降低各個物種建立直系同源外顯子的復(fù)雜性并為之后的進(jìn)化選擇壓分析(Ka/Ks)提供良好的基礎(chǔ).

在建立外顯子區(qū)域之后，通過AS classifier、ESEfinder[11]、RepeatMasker[12]和InterProScan[13]4個模塊對外顯子區(qū)域進(jìn)行注釋.AS classifier的功能是將外顯子分類.雖然Ensembl也提供了外顯子分類信息，但是由于存在許多實(shí)驗(yàn)證據(jù)不充分或可靠度較低的不能翻譯成蛋白質(zhì)的轉(zhuǎn)錄本，因此會對選擇性剪接外顯子產(chǎn)生過高預(yù)估.因此我們利用AS classifier 對外顯子進(jìn)行重新分類.我們將外顯子分成5種類型：Constitutive Exon、Cassette Exon、3′ Splicing Usage、5′ Splicing Usage and Intron Retention[14].ESE Finder 模塊利用ESEfinder 3.0對外顯子區(qū)域進(jìn)行選擇性剪接元件的掃描，參數(shù)使用默認(rèn)值，模塊通過抓取和整理網(wǎng)頁結(jié)果將其存入數(shù)據(jù)庫中.外顯子區(qū)域中的重復(fù)序列對于建立直系同源外顯子也有較大影響，因此我們利用RepeatMasker工具對外顯子區(qū)域的序列進(jìn)行了重復(fù)序列的掃描.為了更好理解選擇性剪接外顯子的存在對于其翻譯的蛋白質(zhì)功能的影響，利用InterproScan模塊通過InterproScan工具對外顯子區(qū)域的序列進(jìn)行了蛋白質(zhì)結(jié)構(gòu)域的預(yù)測，同時對預(yù)測結(jié)果進(jìn)行整理并整合到數(shù)據(jù)庫中.

1.2.3 基因?qū)虻闹毕低赐怙@子區(qū)域的建立

為了建立直系同源外顯子區(qū)域，我們利用Ensembl 直系同源數(shù)據(jù)庫獲取2989個癌癥相關(guān)基因在人、黑猩猩、小鼠、大鼠、鴨嘴獸和斑馬魚中的直系同源基因關(guān)系.考慮到直系同源外顯子的復(fù)雜性，只有在Ensembl 數(shù)據(jù)庫中標(biāo)注為one-to-one 類型的直系同源基因才被用于下一步的分析.對于來自于NCBI數(shù)據(jù)庫的基因序列，其直系同源關(guān)系通過NCBI的homolog數(shù)據(jù)庫進(jìn)行建立.在獲取了一個相對完整的癌癥相關(guān)基因的直系同源基因列表后，通過在直系同源基因間的外顯子區(qū)域進(jìn)行雙向blastn[15]來建立直系同源外顯子區(qū)域.直系同源外顯子區(qū)域之間首先必須滿足互為其比對結(jié)果的最佳匹配，其次比對的E-value值必須小于1×10-5[16].對于沒有比對上的外顯子區(qū)域，我們通過blastn對其在直系同源基因的內(nèi)含子區(qū)域進(jìn)行搜索.通過這種方式有可能找到一些之前未被注釋的外顯子以及一些失去功能的外顯子.同樣，E-value值小于1×10-5作為有意義的結(jié)果的閾值.

1.2.4 進(jìn)化選擇壓

直系同源外顯子區(qū)域的構(gòu)建為外顯子水平的進(jìn)化選擇壓的計(jì)算建立了基礎(chǔ).為了保證外顯子區(qū)域具有正確的開放閱讀框，我們通過選擇該外顯子區(qū)域中的最長的真實(shí)外顯子序列做為該區(qū)域的代表.根據(jù)建立好的直系同源外顯子區(qū)域的對應(yīng)關(guān)系，將這些真實(shí)的外顯子翻譯成的氨基酸序列先通過Muscle[17]進(jìn)行對位排列(alignment)，再通過EMBOSS transalign[18]利用排列好的氨基酸序列指導(dǎo)外顯子核酸序列進(jìn)行對位排列，這樣有效避免了讀碼框的偏移.然后我們利用PAML中的yn00 算法對排列好的外顯子核酸序列進(jìn)行Ka/Ks計(jì)算[19].Ka/Ks值是進(jìn)化選擇壓的一個特征值，在多數(shù)情況下較低的Ka/Ks值(例如<1)代表該區(qū)域處于一個純化選擇的狀態(tài)，而較高的Ka/Ks值(例如>1)則表示該區(qū)域的外顯子序列處于一個較高的進(jìn)化選擇壓力中[20].PAML的yn00共有5種模型來計(jì)算Ka/Ks值，根據(jù)算法時間與名稱分別稱之為NG86、LWL85、LWL85m、LPB93和yn00.5種模型計(jì)算的值通過整理存入數(shù)據(jù)庫中，提供網(wǎng)頁展示.

1.2.5 表達(dá)量水平

為了研究每個癌癥相關(guān)基因在不同物種不同組織間的表達(dá)模式，RNA-Seq高通量數(shù)據(jù)被用于進(jìn)行表達(dá)量的分析.NCBI的SRA數(shù)據(jù)庫中豐富的高通量數(shù)據(jù)資源為表達(dá)量計(jì)算提供了基礎(chǔ)，但SRA中大量的RNA-Seq序列也為篩選數(shù)據(jù)帶來一定困難.為了避免參差不齊的數(shù)據(jù)對數(shù)據(jù)分析的影響，我們設(shè)定了篩選標(biāo)準(zhǔn)：1) 讀長大于50bp;2) pair-end數(shù)據(jù)作為首選;3) 必須是組織樣本.我們從NCBI的SRA數(shù)據(jù)庫中下載了人、黑猩猩、小鼠、大鼠、鴨嘴獸和斑馬魚6個物種的多種正常組織(平均每個物種8個正常組織)的RNA-seq數(shù)據(jù).

同一物種的不同組織的RNA-Seq的原始數(shù)據(jù)選自于SRA數(shù)據(jù)庫中的同一個研究項(xiàng)目，即為分析提供了較為一致的實(shí)驗(yàn)條件，同時也為不同組織之間的表達(dá)水平比較提供了基礎(chǔ).使用TopHat[21]軟件將RNA-Seq的讀長定位到各個物種的基因組上，參數(shù)采用默認(rèn)設(shè)置，軟件使用到的各個物種的基因組數(shù)據(jù)和基因注釋文件來則自于Enseml 68版本.由于RNA-Seq實(shí)驗(yàn)中的PCR擴(kuò)增會對后續(xù)基因表達(dá)量的計(jì)算帶來偏差，因此在計(jì)算表達(dá)量之前首先進(jìn)行PCR擴(kuò)增數(shù)據(jù)的清除步驟.通過Samtools rmdup[22]功能去除pair-end RNA-Seq數(shù)據(jù)的PCR擴(kuò)增產(chǎn)生的讀長，其原理是將定位于基因組上的相同位置的多對讀長只保留其中一對從而達(dá)到去除PCR擴(kuò)增影響的效果.基因和轉(zhuǎn)錄本的表達(dá)量通過Cufflinks[23]軟件進(jìn)行計(jì)算，參數(shù)使用默認(rèn)設(shè)置，Cufflinks通過基因組注釋分別統(tǒng)計(jì)出基因和轉(zhuǎn)錄本的表達(dá)量.由于本次研究更多的關(guān)注外顯子水平，所以我們同時進(jìn)行了外顯子區(qū)域的表達(dá)量計(jì)算.基于Tophat分析結(jié)果中的讀長定位信息以及外顯子區(qū)域的坐標(biāo)，利用Perl自編腳本對外顯子區(qū)域的表達(dá)量進(jìn)行分析.對于定位到多處的讀長，根據(jù)Cufflinks軟件的處理方式將這些讀長平均分配到定位到的所有位置，根據(jù)以上規(guī)則計(jì)算出外顯子區(qū)域的表達(dá)量.最終我們從基因、轉(zhuǎn)錄本、外顯子區(qū)域3個水平統(tǒng)計(jì)了癌癥相關(guān)基因的表達(dá)量.

由于RNA-Seq進(jìn)行的是cDNA的測序，在cDNA片段定位到基因組的過程中，如果該片段處于跨越兩個外顯子的位置，那么就需要將該片段斷開以保證其能跨越內(nèi)含子后定位到基因組上.TopHat在進(jìn)行讀長定位時會同時嘗試去發(fā)現(xiàn)這樣斷開的讀長并確定其剪接位點(diǎn).剪接位點(diǎn)的信息對于我們了解癌癥相關(guān)基因的基因表達(dá)模式具有重要的意義.通過TopHat獲得的剪接位點(diǎn)結(jié)合Ensembl的基因注釋信息，我們將剪接位點(diǎn)分為已知的“known”和未被注釋過的新的“novel”剪接位點(diǎn).跨域剪接位點(diǎn)的讀長作為該位點(diǎn)的支持證據(jù)用于評價(jià)該位點(diǎn)的可信度.這里值得注意的是，由于基因組上存在重復(fù)序列等的影響，可能存在可以定位到多處的讀長，因此對于位于這些位置的新的剪接位點(diǎn)需要更加注意其可靠性，建議新的剪接位點(diǎn)需要10個讀長以上的支持方可用于實(shí)驗(yàn)驗(yàn)證.

SRA數(shù)據(jù)庫中使用Illumina平臺進(jìn)行RNA-Seq測序的人類癌癥組織數(shù)據(jù)用于本研究的癌癥轉(zhuǎn)錄組的分析.目前ASeeDB覆蓋了肺癌和前列腺癌的相關(guān)數(shù)據(jù)(表1).

表1 癌癥樣本信息

我們從SRA數(shù)據(jù)庫中挑選了來自“Lung Cancer Sequencing Project”的高質(zhì)量的肺癌轉(zhuǎn)錄組數(shù)據(jù)(ERP001058)，前列腺癌數(shù)據(jù)來自SRA中“Complete transcriptomic landscape of prostate cancer in Chinese population using RNA-seq”(ERP000550).和計(jì)算正常組織表達(dá)量類似，TopHat軟件用于讀長在基因組上的定位，參數(shù)默認(rèn).實(shí)驗(yàn)中的PCR 擴(kuò)增影響通過Samtools的rmdup 功能進(jìn)行去除.Cuffdif 用于計(jì)算不同基因和轉(zhuǎn)錄的表達(dá)量，參數(shù)默認(rèn)；而外顯子的表達(dá)量通過自編Perl腳本根據(jù)TopHat的讀長定位信息和基因注釋文件進(jìn)行計(jì)算.對于定位到多處的讀長，我們將其平均分配到可能定位的位置.使用Fisher精確檢驗(yàn)進(jìn)行癌癥樣本與正常對照樣本間表達(dá)量差異的統(tǒng)計(jì)檢驗(yàn)，Benjamini-Hochberg修正用于表達(dá)量差異計(jì)算的假陽性控制.

2 結(jié)果與分析

2.1 數(shù)據(jù)庫搜索界面

數(shù)據(jù)庫提供了簡單方便的搜索功能(圖3)，用戶可以選擇基因名字，例如“ABTB1”，Ensembl ID 或者NCBI Unigene ID作為搜索的關(guān)鍵詞，并且數(shù)據(jù)庫支持模糊查詢.同時用戶也可以限制所要檢索的特定物種以獲得更為簡潔的結(jié)果.

2.2 數(shù)據(jù)庫搜索結(jié)果頁面

ASeeDB為用戶提供了詳盡簡潔直觀的搜尋結(jié)果頁面.圖4和圖5展示了一個ASeeDB數(shù)據(jù)庫“ZFYVE9”基因的搜索結(jié)果的頁面.圖4(a)顯示的是搜索基因的基本信息，包括基因的基因名、基因ID、在基因組上的位置和轉(zhuǎn)錄本信息，通過點(diǎn)擊其他物種的名字可以方便的獲得該基因的直系同源基因情況，該基因的轉(zhuǎn)錄本ID和蛋白質(zhì)ID指向相應(yīng)的核酸和蛋白質(zhì)序列；圖4(b)展示的是該基因的外顯子區(qū)域的信息，包括外顯子區(qū)域的對位排列和外顯子注釋信息，其中綠色的方塊代表每個物種直系同源基因的外顯子區(qū)域，灰色的方塊表示該基因的內(nèi)含子與其直系同源基因的外顯子區(qū)域比對有高度相似(E<1×10-3)的序列.外顯子區(qū)域的對位排列使得用戶可以直觀地觀察到外顯子的保守情況，比如圖上列出的人ZFYVE9基因的3～19外顯子區(qū)域在進(jìn)化過程中相對保守，通常情況下相對保守的外顯子更可能具有重要的功能.點(diǎn)擊外顯子區(qū)域的數(shù)字可以展示該基因的轉(zhuǎn)錄本信息以及外顯子區(qū)域的構(gòu)建信息等.如果對該基因的某個外顯子區(qū)域的對位排列的結(jié)果感興趣可以使用該區(qū)域的blast功能進(jìn)行進(jìn)一步證實(shí)，選擇好感興趣的外顯子區(qū)域號和物種，點(diǎn)擊check就可以自動把序列遞交到NCBI進(jìn)行檢索.圖4(c)展示了外顯子區(qū)域的注釋信息，包括外顯子分類、ESE、重復(fù)序列、外顯子區(qū)域的結(jié)構(gòu)域預(yù)測和外顯子區(qū)域的Ka/Ks情況，通過這一部分可以直觀地了解檢索基因中比較重要的外顯子，特別在癌癥研究中，如果這些外顯子發(fā)生突變更可能導(dǎo)致蛋白功能的嚴(yán)重缺失.圖4(c)右側(cè)的柱狀圖是對外顯子區(qū)域的可視化展示，橫坐標(biāo)為外顯子區(qū)域的標(biāo)號，縱坐標(biāo)為每100bp ESE預(yù)測的個數(shù).紅色柱子表示該區(qū)域的外顯子為選擇性剪接外顯子，而綠色的柱子則代表組成型外顯子，同樣地，柱狀圖的高低可以直觀地展示哪些外顯子區(qū)域更可能會被剪接.Ka/Ks部分則提供了Ka/Ks>1 的外顯子區(qū)域信息，可以幫助用戶快速地定位到相應(yīng)的外顯子區(qū)域，這部分外顯子區(qū)域可能處于一個比較強(qiáng)的進(jìn)化選擇壓下，通過點(diǎn)擊該表下方鏈接可以得到由5種模型計(jì)算所得到的全部外顯子區(qū)域的Ka/Ks結(jié)果.

圖5(見第656頁)展示的是檢索到的癌癥相關(guān)基因的正常組織和癌癥組織的表達(dá)信息.表達(dá)量信息的第一部分給出了EBI芯片數(shù)據(jù)鏈接，剪接位點(diǎn)的情況，癌癥表達(dá)量數(shù)據(jù)鏈接和表達(dá)量數(shù)據(jù)來源.通過點(diǎn)擊EBI芯片數(shù)據(jù)鏈接可以方便地鏈接到該基因在EBI Atlas中的芯片表達(dá)信息，以與下面的高通量數(shù)據(jù)做比較.剪接位點(diǎn)的鏈接連接到詳細(xì)的剪接位點(diǎn)信息頁面，而癌癥鏈接提供了癌癥和對照組的詳細(xì)表達(dá)量信息.數(shù)據(jù)來源表格提供了該物種每個組織的RNA-Seq的數(shù)據(jù)，通過SRA accession可以進(jìn)入SRA數(shù)據(jù)庫查看原始RNA-Seq數(shù)據(jù)情況.基因在各個組織中的表達(dá)量通過柱狀圖表示，橫坐標(biāo)為各個組織，縱坐標(biāo)為基因的表達(dá)量，通過FPKM來代表.用戶可以根據(jù)檢索基因在不同組織的中表達(dá)情況，在癌癥研究中選擇合適的組織樣本，并可以通過相應(yīng)的直系同源基因的表達(dá)情況選擇合適的動物模型等.基因轉(zhuǎn)錄本的表達(dá)量通過折線圖來展示；同樣地，橫坐標(biāo)為組織，縱坐標(biāo)為FPKM.由于有些基因的轉(zhuǎn)錄本較多，在圖上很難清楚地觀察某個轉(zhuǎn)錄本在各個組織中的變化趨勢，這時可以通過點(diǎn)擊折線圖的圖例來增減圖中所展示的轉(zhuǎn)錄本的數(shù)據(jù).這部分信息可以幫助用戶在研究基因表達(dá)量的時候選擇合適的轉(zhuǎn)錄本，比如進(jìn)行qPCR實(shí)驗(yàn)時選擇高表達(dá)量的轉(zhuǎn)錄本作為設(shè)計(jì)引物的標(biāo)準(zhǔn).不同組織的外顯子區(qū)域的表達(dá)量同樣采用折線圖進(jìn)行展示，橫坐標(biāo)為外顯子區(qū)域的標(biāo)號，縱坐標(biāo)為FPKM，不同組織通過不同線段的顏色來代表，點(diǎn)擊折線圖圖例可以對展示的組織數(shù)據(jù)進(jìn)行刪減.

癌癥組織表達(dá)數(shù)據(jù)頁面提供了配對的對照和患病3組數(shù)據(jù)，分別對對照和患病組表達(dá)量進(jìn)行了基因、轉(zhuǎn)錄本、外顯子區(qū)域的計(jì)算，通過柱狀圖進(jìn)行展示可以方便地觀察對照和患病在這三個層次上的差異，且提供了詳細(xì)的統(tǒng)計(jì)檢驗(yàn)值，方便用戶篩選到顯著差異的基因.

2.3 數(shù)據(jù)庫總結(jié)

通過ASeeDB數(shù)據(jù)庫，我們利用比較基因組學(xué)的分析方法結(jié)合高通量RNA-seq測序數(shù)據(jù)，系統(tǒng)地總結(jié)了癌癥相關(guān)基因的選擇性剪接外顯子的功能和進(jìn)化關(guān)系.數(shù)據(jù)庫基本數(shù)據(jù)見表2.

表2 數(shù)據(jù)庫基本信息統(tǒng)計(jì)

ASeeDB包含癌癥相關(guān)基因外顯子區(qū)域的進(jìn)化保守性、Ka/Ks值以及基因、轉(zhuǎn)錄本、外顯子區(qū)域的表達(dá)量統(tǒng)計(jì)等信息，結(jié)合這些信息用戶可以方便地找到一些有研究價(jià)值的基因或者外顯子.外顯子區(qū)域蛋白質(zhì)結(jié)構(gòu)域預(yù)測可以指出其可能的功能，而外顯子是否存在選擇性剪接的信息又可以用來推斷包含它的轉(zhuǎn)錄本是否具有相似的功能，以及剪接是否會對蛋白質(zhì)功能產(chǎn)生影響.數(shù)據(jù)庫提供的物種間的序列比對可以幫助用戶發(fā)現(xiàn)沒有注釋的外顯子區(qū)域或者是保留的失去功能的外顯子.相比于基因?qū)用娴腒a/Ks，數(shù)據(jù)庫提供的外顯子層面的Ka/Ks，對于發(fā)現(xiàn)適應(yīng)性進(jìn)化事件具有更高的敏感性，可以更方便地來探索基因中未知功能的區(qū)域[24].

3 討 論

ASeeDB 是一個基于選擇性剪接分析工具通過比較基因組學(xué)方法對癌癥相關(guān)基因研究的數(shù)據(jù)庫.傳統(tǒng)的選擇性剪接數(shù)據(jù)庫通常包含各個物種所有基因的選擇性剪接數(shù)據(jù).但是在癌癥研究中，很多正常基因的不正常的剪接形式或者剪接調(diào)控通常是由于腫瘤發(fā)生過程中的結(jié)果而不是原因，這會對不同癌癥的研究帶來更大的復(fù)雜性.而很多癌癥相關(guān)基因的不正確的剪接已經(jīng)被證實(shí)與癌癥的發(fā)生和發(fā)展相關(guān)，因此ASeeDB主要關(guān)注癌癥相關(guān)基因的選擇性剪接事件.

Ensembl和NCBI作為主要的生物學(xué)數(shù)據(jù)庫具有多樣的基因數(shù)據(jù)資源，但也有各自的不足.例如，Ensembl只包含有位于primary assembly參考基因組上的基因，而NCBI的homology 數(shù)據(jù)庫沒有Ensembl的直系同源數(shù)據(jù)庫清楚明了.因此Ensembl和NCBI通過整合來保證ASeeDB數(shù)據(jù)的準(zhǔn)確和完整.NCBI的SRA數(shù)據(jù)庫提供了豐富的二代測序數(shù)據(jù).但是由于二代測序數(shù)據(jù)的多樣性，包括讀段長度、平臺以及測序深度的不同會對后續(xù)的生物信息分析帶來各種各樣的偏向性.為了盡量避免這種情況的發(fā)生，我們制定了同一的篩選標(biāo)準(zhǔn)：必須是Illumina平臺且讀長必須超過50bp.通過篩選，6個物種共49個正常組織的RNA-seq數(shù)據(jù)(平均每個物種8個組織)，以及2種癌癥(6個SRA RNA-Seq數(shù)據(jù))被用于ASeeDB數(shù)據(jù)庫表達(dá)量部分的分析.ASeeDB可以方便地查詢癌癥相關(guān)基因的進(jìn)化保守性，進(jìn)化選擇壓和表達(dá)量等信息，除了一些基本的統(tǒng)計(jì)信息，數(shù)據(jù)庫其他方面也可以應(yīng)用于特定基因的分析.通過對數(shù)據(jù)庫中的信息進(jìn)行數(shù)據(jù)挖掘，我們成功地發(fā)現(xiàn)了一些基因的外顯子發(fā)生了具有一定生物學(xué)意義的選擇性剪接，例如PPAR-gamma的4個選擇性剪接外顯子承受著不同的進(jìn)化選擇壓力，其中Ka/Ks>1的一個外顯子可能與嚙齒類的脂肪儲存相關(guān).表3列出了我們通過數(shù)據(jù)挖掘發(fā)現(xiàn)的Ka/Ks>1的一些外顯子信息.

表3 ASeeDB 數(shù)據(jù)庫中提示的可能帶有特殊功能的選擇性剪接基因/外顯子

[1] CHOW L T, GELINAS R E, BROKER T R,etal. An amazing sequence arrangement at the 5′ ends of adenovirus 2 messenger RNA [J].Cell,1977,12(1)：1-8.

[2] BERGET S M, MOORE C, SHARP P A. Spliced segments at the 5′ terminus of adenovirus 2 late mrna [J].ProcNatlAcadSciUSA,1977,74(8)：3171-3175.

[3] CARSTENS R P, WAGNER E J, GARCIA-BLANCO M A. An intronic splicing silencer causes skipping of the iiib exon of fibroblast growth factor receptor 2 through involvement of polypyrimidine tract binding protein [J].MolecularandCellularBiology,2000,20(19)：7388-7400.

[4] POZZOLI U, SIRONI M. Silencers regulate both constitutive and alternative splicing events in mammals [J].CellularandMolecularLifeSciences,2005,62(14)：1579-1604.

[5] WOOLARD J, WANG W Y, BEVAN H S,etal. VEGF165b, an inhibitory vascular endothelial growth factor splice variant：mechanism of action,invivoeffect on angiogenesis and endogenous protein expression [J].CancerResearch,2004,64(21)：7822-7835.

[6] SCHAAL T D, MANIATIS T. Multiple distinct splicing enhancers in the protein-coding sequences of a constitutively spliced pre-mrna [J].MolecularandCellularBiology,1999,19(1)：261-273.

[7] LANGE T, GUTTMANN-RAVIV N, BARUCH L,etal. VEGF162, a new heparin-binding vascular endothelial growth factor splice form that is expressed in transformed human cells [J].TheJournalofBiologicalChemistry,2003,278(19)：17164-17169.

[8] TAZI J, BAKKOUR N, STAMM S. Alternative splicing and disease [J].BiochimicaetBiophysicaActa,2009,1792(1)：14-26.

[9] BIANCHI V, COLANTONI A, CALDERONE A,etal. DBATE：Database of alternative transcripts expression [J].Database：theJournalofBiologicalDatabasesandCuration,2013(6)：1843-1852.

[10] FLICEK P, AMODE M R, BARRELL D,etal. Ensembl 2012 [J].NucleicAcidsResearch,2012,40(Database issue)：84-90.

[11] CARTEGNI L, WANG J, ZHU Z,etal. ESEfinder：A web resource to identify exonic splicing enhancers [J].NucleicAcidsResearch,2003,31(13)：3568-3571.

[12] TARAILO-GRAOVAC M, CHEN N S. Using repeatmasker to identify repetitive elements in genomic sequences [D/OL]. Curr Protoc Bioinformatics,2009.doi：10.1002/0471250953.bi0410s25.

[13] ZDOBNOV E M, APWEILER R. Interproscan—an integration platform for the signature-recognition methods in interpro [J].Bioinformatics,2001,17(9)：847-848.

[14] MCMANUS C J, GRAVELEY B R. RNA structure and the mechanisms of alternative splicing [J].CurrentOpinioninGenetics&Development,2011,21(4)：373-379.

[15] ALTSCHUL S F, MADDEN T L, SCHAFFER A A,etal. Gapped BLAST and PSI-BLAST：a new generation of protein database search programs [J].NucleicAcidsResearch,1997,25(17)：3389-3402.

[16] FU C L, LIN W C. Identification of gene-oriented exon orthology between human and mouse [J].BMCGenomics,2012,13(1)：1-10.

[17] EDGAR R C. MUSCLE：Multiple sequence alignment with high accuracy and high throughput [J].NucleicAcidsResearch,2004,32(5)：1792-1797.

[18] RICE P, LONGDEN I, BLEASBY A. EMBOSS：the European Molecular Biology Open Software Suite [J].TrendsinGenetics,2000,16(6)：276-277.

[19] YANG Z. PAML：A program package for phylogenetic analysis by maximum likelihood [J].ComputerApplicationsintheBiosciences,1997,13(5)：555-556.

[20] YANG Z, BIELAWSKI J P. Statistical methods for detecting molecular adaptation [J].TrendsinEcology&Evolution,2000,15(12)：496-503.

[21] TRAPNELL C, PACHTER L, SALZBERG S L. Tophat：Discovering splice junctions with RNA-Seq [J].Bioinformatics,2009,25(9)：1105-1111.

[22] LI H, HANDSAKER B, WYSOKER A,etal. The Sequence Alignment/Map format and samtools [J].Bioinformatics,2009,25(16)：2078-2079.

[23] TRAPNELL C, ROBERTS A, GOFF L,etal. Differential gene and transcript expression analysis of RNA-seq experiments with tophat and Cufflinks [J].NatureProtocols,2012,7(3)：562-578.

[24] CLARK A G, GLANOWSKI S, NIELSEN R,etal. Inferring nonneutral evolution from human-chimp-mouse orthologous gene trios [J].Science,2003,302(5652)：1960-1963.

[25] CHAO C C, CHANG P Y, LU H H. Human Gas7 isoforms homologous to mouse transcripts differentially induce neurite outgrowth [J].JournalofNeuroscienceResearch,2005,81(2)：153-162.

[26] ARLT A, SCHAFER H. Role of the immediate early response 3 (IER3) gene in cellular stress response, inflammation and tumorigenesis [J].EuropeanJournalofCellBiology,2011,90(6/7)：545-552.

[27] ALBERTI L, BACHELOT T, DUC A,etal. A spliced isoform of interleukin 6 mrna produced by renal cell carcinoma encodes for an interleukin 6 inhibitor [J].CancerResearch,2005,65(1)：2-5.

[28] MICHELS J, JOHNSON P W, PACKHAM G. Mcl-1 [J].TheInternationalJournalofBiochemistry&CellBiology,2005,37(2)：267-271.

ASeeDB：A Comparative Genomic Database for Alternative Splicing Exon Evolution

XU Minjie1,2，*, DOU Tonghai1，*, XU Jiaxi1, JIANG Jianping1,2, LIU Cong1, FU Maobin1,GAO Yuan1, CHEN Chengwen2, ZHANG Liang2, ZHOU Yan1,2

(1. State Key Laboratory of Genetic Engineering, Department of Microbiology and Microbial Engineering,SchoolofLifeSciences,FudanUniversity,Shanghai200438,China; 2.Shanghai-MOSTKeyLaboratoryofHealthandDiseaseGenomics,ChineseNationalHumanGenomeCenteratShanghai,Shanghai201203,China)

Alternative splicing(AS) is an essential eukaryotic gene regulatory mechanism associated with various types of physiological and pathological activities. The incorrect or mis-regulated alternative splicing of cancer genes has been associated with the development of multiple sorts of cancers. Moreover, as the basic element of AS evolution, AS exons show various evolutionary features across species and multiple studies have been made to elucidate these characteristics. In this research, we combine Illumina RNA sequencing data with comparative genomic tools to establish a cancer gene database called Alternative Splicing Exon Evolution Database of Cancer Genes(ASeeDB). We have systematically summarized the common features for 2 989 cancer genes AS exon evolution by combining comparative genomic analysis. General features including evolutionary conservation, domain prediction,Ka/Ksvalue and expression level have been outlined, together with some promising genes/exons for future research. Domain prediction of exons can point out putative function and critical region of cancer genes; Cross species sequence alignments are able to find unannotated exons and dysfunctional remains of exons. In addition, Protein domain predictions point out known putative functional regions in genes. Comparing to averaged gene-levelKa/Ksvalues, Exon-levelKa/Ksvalues are more informative in terms of finding adaptive evolution events, which suggest unknown functional regions in genes.

alternative splicing； evolution； database； exon

0427-7104(2016)05-0649-11

2015-12-22

國家科技基礎(chǔ)專項(xiàng)課題(2009FY120100),國家高技術(shù)研究發(fā)展計(jì)劃(2012AA020409)；國家自然科學(xué)基金(31071158)

徐敏杰(1987—)，男，碩士研究生；*并列第一作者；周 雁，男，副教授，通訊聯(lián)系人，E-mail: zhouy@fudan.edu.cn.

TP 311.135.1

A