基于組織特異性基因集權(quán)重的分析方法研究

田 甜, 王化琨

(黑龍江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 哈爾濱 150080)

0 引 言

基因集檢驗(yàn)或通路分析是一種功能強(qiáng)大且被廣泛采用的基因分析方法，可用于分析和解釋高維基因組數(shù)據(jù)[1-2]?；蚣瘷z驗(yàn)可使研究人員從研究單個(gè)基因變量的水平拓展到基因集合的多變量水平，進(jìn)而探索具有生物學(xué)意義的基因組關(guān)聯(lián)情況，如：參與特定代謝通路的關(guān)聯(lián)基因。 基于特定功能相關(guān)的基因組變量進(jìn)行的統(tǒng)計(jì)學(xué)檢驗(yàn)比基于單個(gè)基因變量的檢驗(yàn)有更多優(yōu)勢(shì)，包括改善的統(tǒng)計(jì)功效、更直觀的生物學(xué)解釋等[3-5]。鑒于這些優(yōu)勢(shì)，研究人員在過(guò)去的15年間投入了大量的精力，開發(fā)出很多有效的基因集檢驗(yàn)方法[6-9]，如ROAST方法，它是自檢驗(yàn)方法，運(yùn)用了一種基于蒙特卡洛(Monte Carlo)多變量技術(shù)原理，對(duì)數(shù)據(jù)進(jìn)行旋轉(zhuǎn)來(lái)替代傳統(tǒng)的交換排列陣的方法。再如CAMERA方法，它是競(jìng)爭(zhēng)性檢驗(yàn)方法，利用線性模型估算出基因集中平均的基因關(guān)聯(lián)。在改善基因集檢驗(yàn)方法的同時(shí)，研究人員也在建立大型基因集公共存儲(chǔ)庫(kù)方面取得了不錯(cuò)的進(jìn)展, 如最常用的GO數(shù)據(jù)庫(kù)[10-12]。但是基因集檢驗(yàn)在實(shí)際應(yīng)用中仍然受到一些限制，包括基因注釋的質(zhì)量、統(tǒng)計(jì)能力和組織特異性等方面。

組織特異性的存在是取決于構(gòu)成該組織主要成分的細(xì)胞性質(zhì)。隨著微陣列技術(shù)的成熟運(yùn)用，通過(guò)基因表達(dá)的全基因組分析，人們發(fā)現(xiàn)基因表達(dá)的功能與組織密切相關(guān)，并且一些普遍存在的生物過(guò)程也發(fā)生在特定的組織當(dāng)中。隨著人類蛋白圖譜(HPA)[13]和基因組織表達(dá)工程(GTEx)[14]等項(xiàng)目的開展，人們對(duì)特定組織內(nèi)的基因活動(dòng)的認(rèn)識(shí)也逐漸加深。現(xiàn)在，通過(guò)質(zhì)譜分析和免疫組織化學(xué)等技術(shù)，可以對(duì)完整的人類基因組序列分析中鑒定出的約20 687個(gè)蛋白質(zhì)編碼基因的組織特異性活動(dòng)進(jìn)行研究。如在HPA中所述，在所有蛋白質(zhì)編碼基因中，有34%的蛋白質(zhì)在至少一種組織中表達(dá)升高，其中有17%在特定組織中的mRNA水平至少是所有組織中平均水平的5倍，還有約44%的蛋白編碼基因在所有組織中表達(dá)，并且在任何組織中都沒有升高。

盡管對(duì)人類基因的組織特異性活動(dòng)的描述花費(fèi)了大量精力，但是在基因集檢驗(yàn)中卻很少利用組織特異性信息，如分子簽名數(shù)據(jù)庫(kù)[12](MSigDB)完全缺乏有關(guān)基因集的組織特異性和基因集的注釋信息。由于組織特定版本的基因集不易創(chuàng)建，通常使用的作法為：在執(zhí)行基因集檢驗(yàn)時(shí)使用通用的基因集進(jìn)行，而不去考慮實(shí)驗(yàn)組織的類型。如果基因在組織中都是無(wú)所不在的表達(dá)，忽視組織特異性進(jìn)行基因集檢驗(yàn)對(duì)結(jié)果影響不大。但除了持家基因是在所有組織中具有類似的表達(dá)水平之外，大多數(shù)基因在不同組織中表達(dá)水平是有差異性的，而在此條件下進(jìn)行基因集檢驗(yàn)，會(huì)在一定程度上提高Ⅰ型和Ⅱ型錯(cuò)誤率。為了克服這種基因集檢驗(yàn)的缺點(diǎn)，本研究使用來(lái)自人類蛋白圖譜的組織特異性基因活性信息和分子簽名數(shù)據(jù)庫(kù)中的所有過(guò)濾后的基因集合生成了組織特異性基因集權(quán)重的集合[15]，并且利用這些組織特性的基因集權(quán)重對(duì)p值進(jìn)行加權(quán)，以這種方式進(jìn)行組織特異性的基因集檢驗(yàn)。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

(1)人類蛋白質(zhì)圖譜HPA[13]

人類蛋白質(zhì)圖譜是一項(xiàng)基于瑞典的科研計(jì)劃，它始于2003年，旨在利用多種組學(xué)技術(shù)的整合來(lái)繪制細(xì)胞、組織和器官中所有人類蛋白質(zhì)的圖譜，包括基于抗體的成像，基于質(zhì)譜的蛋白質(zhì)組學(xué)、轉(zhuǎn)錄組學(xué)和系統(tǒng)生物學(xué)。本文主要是使用組織圖譜，該圖譜是基于對(duì)37種主要的人體正常組織類型的RNA(RNA-seq) 進(jìn)行深度測(cè)序，并在包含44種不同組織類型的組織微陣列上進(jìn)行免疫組織化學(xué)分析，它包含人類基因在mRNA和蛋白質(zhì)水平上的表達(dá)譜信息。本文關(guān)于人類蛋白質(zhì)編碼基因的組織特異性活性信息是從v18.1版本的人類蛋白質(zhì)圖譜中下載的。

① IHC蛋白豐度數(shù)據(jù)：是基于免疫組化學(xué)和組織微陣列的蛋白在正常人體組織的表達(dá)譜，該數(shù)據(jù)包括基因標(biāo)識(shí)符、組織名稱、注釋細(xì)胞類型、表達(dá)值和表達(dá)值的基因可靠性，根據(jù)表中的蛋白表達(dá)值計(jì)算數(shù)值得分記為ProteinScore(其中未檢測(cè)到=0，低=0.5，中等=1.0,高=2.0)。

② RNA-seq數(shù)據(jù)：是基于RNA-seq的37個(gè)組織的RNA水平，該數(shù)據(jù)包括基因標(biāo)識(shí)符、分析樣本(組織)和每百萬(wàn)轉(zhuǎn)錄本。根據(jù)表中的TPM值，計(jì)算折疊后的TPM值，記為FoldAboveMean, 該組織中的平均TPM值表示為mean(TPM),則FoldAboveMean=TPM/mean(TPM)。然后根據(jù)基因名和組織名將兩個(gè)列表結(jié)合成一個(gè)表。

(2) 分子簽名數(shù)據(jù)庫(kù)MSigDB[12]

注釋齊全的基因集代表了生物過(guò)程的整體，對(duì)于解釋大規(guī)模基因組數(shù)據(jù)至關(guān)重要，分子簽名數(shù)據(jù)庫(kù)是此類集合中使用最廣泛的存儲(chǔ)庫(kù)之一。分子簽名數(shù)據(jù)庫(kù)是將人類基因從位置、功能、代謝途徑和靶標(biāo)結(jié)合等多種角度出發(fā)，構(gòu)建了許多基因集合，其中的一個(gè)基因集合中包含了具有相近位置或類似功能的許多基因。分子簽名數(shù)據(jù)庫(kù)涵蓋了很多種類的基因集合和更廣泛的基因集來(lái)源和類型，包括從原始研究出版物中提取的簽名以及從GO[10]和KEGG[11]等專業(yè)資源中提取的整套集合。該數(shù)據(jù)庫(kù)中的基因集是通過(guò)人工篩選和自動(dòng)計(jì)算兩種方法獲得的。最初的分子簽名數(shù)據(jù)庫(kù)是在2005年發(fā)布在GSEA軟件上，共有1 325套，并且該數(shù)據(jù)庫(kù)是不斷在更新的。

本文使用的基因集是從分子簽名數(shù)據(jù)庫(kù)(MSigDB)的v7.1版本下載的，包含了13個(gè)不同類別的集合，共有25 824個(gè)基因集。

1.2 分析方法

本文使用了來(lái)自人類蛋白質(zhì)圖譜的組織特異性基因活性信息來(lái)計(jì)算分子簽名數(shù)據(jù)庫(kù)中集合的組織特異性基因集權(quán)重，該方法如圖1所示。由圖可知，來(lái)自分子簽名數(shù)據(jù)庫(kù)的集合被表示為一個(gè)矩陣，行代表基因集，列代表基因，若基因和基因集之間存在注釋，則元素記為1，然后使用人類蛋白質(zhì)圖譜中在不同組織的基因活性信息，根據(jù)公式(1)計(jì)算出組織特異性基因權(quán)重，最后將分子簽名數(shù)據(jù)庫(kù)中的基因集合作為輸入，利用下述步驟，得到組織特異性基因集權(quán)重。

圖1 組織特異性基因集權(quán)重計(jì)算的方法示意圖

1.2.1 組織特異性基因權(quán)重的計(jì)算

(1)

式中：ei,t表示基因i在組織t中的RNA-seq數(shù)據(jù)折疊后的TPM值，單位為每千堿基轉(zhuǎn)錄本片段數(shù)/百萬(wàn)片段映射，若基因i在組織t中的RNA-sep數(shù)據(jù)缺失，則ei,t記為1；ai,t表示基因i在組織t中的IHC蛋白豐度數(shù)據(jù)，若基因i在組織t中的蛋白質(zhì)豐度數(shù)據(jù)缺失，則ai,t記為1。式(1)生成了組織特異性基因權(quán)重，需要在蛋白質(zhì)和RNA水平的證據(jù)才能產(chǎn)生非零值。

1.2.2 組織特異性基因集權(quán)重的計(jì)算

(2)

式中m表示注釋到基因集j的基因，m={i=1,2,…,p}且Gj,i=1，其中Gj,i=1表示基因i注釋到基因集j中，|m|為此集合的尺寸。mc為補(bǔ)集，mc={1,2,…,p}且Gj,i=0，其中Gj,i=0表示基因i未注釋到基因集j中，|mc|為補(bǔ)集的尺寸。

采用單邊雙樣本的t檢驗(yàn)，得到了p值，則有：

(3)

式中pvalj,t是t檢驗(yàn)中得到的p值，是以多大的誤差拒絕原假設(shè)，原假設(shè)為組織t中注釋到基因集j的基因的平均權(quán)重等于不在基因集j中的平均權(quán)重。取p值的負(fù)對(duì)數(shù)就得到了組織特異性的基因集權(quán)重,并且本文還在假設(shè)檢驗(yàn)前對(duì)MSigDB數(shù)據(jù)庫(kù)中的13個(gè)不同類別的基因集都進(jìn)行了過(guò)濾，保留了基因集中基因個(gè)數(shù)在10～200之間的基因集合。

2 結(jié)果與分析

2.1 組織特異性基因集權(quán)重

在假設(shè)檢驗(yàn)之前對(duì)分子簽名數(shù)據(jù)庫(kù)中的13個(gè)不同類別基因集都進(jìn)行了過(guò)濾，如表1所示。這是為了解決某些基因集的組織特異性權(quán)重相對(duì)于其他基因集的組織特異性權(quán)重較大的問(wèn)題，在這種情況下，一些不重要的p值會(huì)在基因集檢驗(yàn)下生成顯著的FDR值，從而影響結(jié)果的準(zhǔn)確性，所以將基因集權(quán)重離散化，即對(duì)基因集進(jìn)行過(guò)濾。

表1 MSigDB v7.1版本的所有13個(gè)不同集合

使用上述分析方法，為分子簽名數(shù)據(jù)庫(kù)中的13個(gè)不同類別基因集合和人類蛋白質(zhì)圖譜所支持的37種人類組織類型，生成了組織特異性基因集權(quán)重。人類蛋白質(zhì)圖譜支持的組織類型有脂肪組織、大腦皮層、骨髓、肝臟、腎臟等37種組織類型。這些組織特異性的基因集權(quán)重，在下文中還會(huì)使用到。

2.2 組織特異性基因集檢驗(yàn)

當(dāng)僅檢驗(yàn)一個(gè)基因集時(shí)，p值是統(tǒng)計(jì)顯著性的適當(dāng)度量，但是當(dāng)檢驗(yàn)包含數(shù)千個(gè)基因集的大型基因集時(shí)，在基因集中可能會(huì)出現(xiàn)看似很高的假陽(yáng)性的p值，這就被稱為多重假設(shè)檢驗(yàn)的問(wèn)題。

2.2.1 白血病數(shù)據(jù)

GSE131184數(shù)據(jù)集是來(lái)源美國(guó)國(guó)家生物技術(shù)信息中心(NCBI)網(wǎng)站中的基因表達(dá)數(shù)據(jù)庫(kù)(GEO),平臺(tái)號(hào)是GPL570，該數(shù)據(jù)集是從急性髓系白血病(AML)或T急性淋巴細(xì)胞白血病(T-all)患者的骨髓樣本中獲取的總RNA，樣品分為兩個(gè)批次運(yùn)行(標(biāo)記為L(zhǎng)或M)，共有125個(gè)樣本。

2.2.2 重度抑郁癥數(shù)據(jù)

GSE54563數(shù)據(jù)集是來(lái)源美國(guó)國(guó)家生物技術(shù)信息中心(NCBI)網(wǎng)站中的基因表達(dá)數(shù)據(jù)庫(kù)(GEO),平臺(tái)號(hào)是GPL6947，該數(shù)據(jù)集是關(guān)于重度抑郁癥的，表達(dá)數(shù)據(jù)是從人類大腦前扣帶皮層組織中獲取的，包含了25對(duì)共50個(gè)樣本的對(duì)照樣本和重度抑郁癥樣本。

2.2.3 二型糖尿病數(shù)據(jù)

GSE73034數(shù)據(jù)集是來(lái)源美國(guó)國(guó)家生物技術(shù)信息中心(NCBI)網(wǎng)站中的基因表達(dá)數(shù)據(jù)庫(kù)(GEO),平臺(tái)號(hào)是GPL6480，該數(shù)據(jù)集是關(guān)于二型糖尿病的，該數(shù)據(jù)比較了瘦、肥胖胰島素敏感(OIS)、肥胖胰島素抵抗(OIR)和肥胖T2D患者肌肉活檢的基因表達(dá)差異。表達(dá)數(shù)據(jù)從人類骨骼肌中獲取，每組有7個(gè)樣本，共28個(gè)樣本。

2.2.4 結(jié)果

以上三個(gè)數(shù)據(jù)集的基因表達(dá)譜不能直接使用，所以利用R語(yǔ)言中的Bioconductor進(jìn)行數(shù)據(jù)預(yù)處理，將表達(dá)譜中的探針集注釋到對(duì)應(yīng)的基因上，得到處理好后的歸一化基因表達(dá)數(shù)據(jù)，然后使用兩階段的競(jìng)爭(zhēng)性基因集檢驗(yàn)方法CAMERA進(jìn)行基因集檢驗(yàn)[19]，該方法可以隨基因集成員之間的相關(guān)性進(jìn)行調(diào)整，這種競(jìng)爭(zhēng)性的基因集檢驗(yàn)形式假定了基因水平權(quán)重的獨(dú)立性，并且適用于檢驗(yàn)包含許多基因集的大型數(shù)據(jù)庫(kù)中哪些基因集相對(duì)于其它基因集整體變化更為顯著，也適用于找出具有意義的基因集，此方法由Limma包中的Camera功能實(shí)現(xiàn)[20]。本次基因集檢驗(yàn)使用到的基因集合是分子簽名數(shù)據(jù)庫(kù)中過(guò)濾后的C2.CP集合，并使用了Camera中的默認(rèn)設(shè)置，然后將使用Camera方法后得到的FDR值與使用BH方法控制的加權(quán)FDR值進(jìn)行比較，該結(jié)果見表3。這三種類型的疾病在q值小于0.2的情況下的加權(quán)的FDR分析中產(chǎn)生了更多的發(fā)現(xiàn)，如表中所示白血病基因表達(dá)數(shù)據(jù)在未加權(quán)的情況下有32個(gè)發(fā)現(xiàn)，加權(quán)后有44個(gè)發(fā)現(xiàn)；重度抑郁癥基因表達(dá)數(shù)據(jù)在未加權(quán)的情況下有1個(gè)發(fā)現(xiàn)，而加權(quán)后有18個(gè)發(fā)現(xiàn)；二型糖尿病基因表達(dá)數(shù)據(jù)集在未加權(quán)的情況下產(chǎn)生了1個(gè)發(fā)現(xiàn)，而加權(quán)后有2個(gè)發(fā)現(xiàn)。在加權(quán)后wFDR分析產(chǎn)生了額外的基因集發(fā)現(xiàn)，一般來(lái)說(shuō)，這個(gè)結(jié)果在生物學(xué)上是合理的。

p值加權(quán)在以下兩種情況下最有效，一是當(dāng)基因集檢驗(yàn)的目的是識(shí)別目標(biāo)組織中起重要作用的基因集，二是在基因集檢驗(yàn)中因變量與被分析組織的正?；顒?dòng)和組織特異性過(guò)程顯著關(guān)聯(lián)。由于權(quán)重較大的基因集更能反映正常活動(dòng)和組織特異性，在這兩種情況下，p值加權(quán)能更好地提高統(tǒng)計(jì)功效。

表3 組織特異性基因集檢驗(yàn)結(jié)果

表4列出了使用上述方法進(jìn)行wFDR分析骨骼肌相對(duì)于二型糖尿病產(chǎn)生的兩個(gè)基因集發(fā)現(xiàn)。

表4 在骨骼肌中與二型糖尿病有關(guān)的通路

表5 肺(左)、腎臟(右)在C2.CP中權(quán)重前十的基因集

基因集權(quán)重REACTOME_ORGANIC_CATION_ANION_ZWITTERION_TRANSPORT128REACTOME_TRANSPORT_OF_BILE_SALTS_AND_ORGANIC_ACIDS_METAL_IONS_AND_AMINE_COMPOUNDS89KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM77REACTOME_GLYOXYLATE_METABOLISM_AND_GLYCINE_DEGRADATION71REACTOME_TRANSPORT_OF_INORGANIC_CATIONS_ANIONS_AND_AMINO_ACIDS_OLIGOPEPTIDES58REACTOME_SLC_TRANSPORTER_DISORDERS49KEGG_ARGININE_AND_PROLINE_METABOLISM39REACTOME_NA_CL_DEPENDENT_NEUROTRANSMITTER_TRANSPORTERS38REACTOME_PYRIMIDINE_CATABOLISM33KEGG_PROXIMAL_TUBULE_BICARBONATE_RECLAMATION29

2.3 人類組織的功能表征

組織特異性基因集權(quán)重可以用來(lái)分析相關(guān)人類組織的功能表征。具體說(shuō)，對(duì)分子簽名數(shù)據(jù)庫(kù)中的某類集合如C2.CP集合，按照1.2節(jié)所述步驟，對(duì)分配給給定組織的每個(gè)基因集合的權(quán)重進(jìn)行排序，權(quán)重大的基因集更能反映組織中活躍的主要生物過(guò)程。如對(duì)分子簽名數(shù)據(jù)庫(kù)中的典型通路(C2.CP)中的肺組織和腎臟組織進(jìn)行分析，表5顯示了給定肺組織、腎臟組織在分子簽名數(shù)據(jù)庫(kù)中的C2.CP集合中排名前十的基因集合。對(duì)給定的分子簽名數(shù)據(jù)庫(kù)中的集合名可以通過(guò)GSEA軟件查找，分子簽名數(shù)據(jù)庫(kù)中對(duì)該集合的描述，或者通過(guò)外部鏈接到其他數(shù)據(jù)庫(kù)(如Reactome信號(hào)通路數(shù)據(jù)庫(kù))的描述信息。Reactome是一個(gè)免費(fèi)開源的通路數(shù)據(jù)庫(kù)，提供直觀的生物信息學(xué)工具，用于可視化、解釋和分析通路相關(guān)知識(shí)，以支持基礎(chǔ)研究、基因組分析、建模和系統(tǒng)生物學(xué)研究等。根據(jù)這些數(shù)據(jù)庫(kù)中的信息可知，權(quán)重高的基因集合捕捉到了與該組織有關(guān)的已知生物學(xué)特性和過(guò)程，識(shí)別出了肺組織中的代謝通路，如REACTOME_SURFACTANT_METABOLISM是與肺組織表面活性物質(zhì)有關(guān)的新陳代謝通路，也識(shí)別出與腎臟有關(guān)的運(yùn)輸通路，如REACTOME_ORGANIC_CATION_ANION_ZWITTERION_TRANSPORT是有機(jī)陽(yáng)離子/陰離子/兩性離子轉(zhuǎn)運(yùn)通路，這些轉(zhuǎn)運(yùn)蛋白在腎臟組織中表達(dá)。

3 結(jié) 論

使用人類蛋白質(zhì)圖譜的組織特異性基因活性信息和過(guò)濾后的分子簽名數(shù)據(jù)庫(kù)中的基因集生成了組織特異性的基因集權(quán)重。為了避免權(quán)重較大的基因集影響結(jié)果的準(zhǔn)確性，將分子簽名數(shù)據(jù)中的基因集合進(jìn)行過(guò)濾。利用組織特異性基因集權(quán)重對(duì)p值進(jìn)行加權(quán)并用BH方法提供FDR控制，進(jìn)行wFDR分析，在不同組織的三種疾病下加權(quán)的FDR一共有64個(gè)發(fā)現(xiàn)，而未加權(quán)的FDR一共有34個(gè)發(fā)現(xiàn)。在p值小于0.2的假設(shè)下，使用組織特異性基因集權(quán)重對(duì)p值進(jìn)行加權(quán)能產(chǎn)生更多的發(fā)現(xiàn)。使用組織特異性基因集權(quán)重可以提高基因集檢驗(yàn)的統(tǒng)計(jì)功效，并且也可以提高在高維基因組數(shù)據(jù)的實(shí)驗(yàn)中識(shí)別生物學(xué)上有關(guān)的基因集關(guān)聯(lián)的信息。