基于組織特異性基因集權重的分析方法研究

田 甜, 王化琨

(黑龍江大學 數學科學學院, 哈爾濱 150080)

0 引 言

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

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

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

1 材料與方法

1.1 數據來源

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

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

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

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

(2) 分子簽名數據庫MSigDB[12]

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

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

1.2 分析方法

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

圖1 組織特異性基因集權重計算的方法示意圖

1.2.1 組織特異性基因權重的計算

(1)

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

1.2.2 組織特異性基因集權重的計算

(2)

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

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

(3)

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

2 結果與分析

2.1 組織特異性基因集權重

在假設檢驗之前對分子簽名數據庫中的13個不同類別基因集都進行了過濾，如表1所示。這是為了解決某些基因集的組織特異性權重相對于其他基因集的組織特異性權重較大的問題，在這種情況下，一些不重要的p值會在基因集檢驗下生成顯著的FDR值，從而影響結果的準確性，所以將基因集權重離散化，即對基因集進行過濾。

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

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

2.2 組織特異性基因集檢驗

當僅檢驗一個基因集時，p值是統計顯著性的適當度量，但是當檢驗包含數千個基因集的大型基因集時，在基因集中可能會出現看似很高的假陽性的p值，這就被稱為多重假設檢驗的問題。

2.2.1 白血病數據

GSE131184數據集是來源美國國家生物技術信息中心(NCBI)網站中的基因表達數據庫(GEO),平臺號是GPL570，該數據集是從急性髓系白血病(AML)或T急性淋巴細胞白血病(T-all)患者的骨髓樣本中獲取的總RNA，樣品分為兩個批次運行(標記為L或M)，共有125個樣本。

2.2.2 重度抑郁癥數據

GSE54563數據集是來源美國國家生物技術信息中心(NCBI)網站中的基因表達數據庫(GEO),平臺號是GPL6947，該數據集是關于重度抑郁癥的，表達數據是從人類大腦前扣帶皮層組織中獲取的，包含了25對共50個樣本的對照樣本和重度抑郁癥樣本。

2.2.3 二型糖尿病數據

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

2.2.4 結果

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

p值加權在以下兩種情況下最有效，一是當基因集檢驗的目的是識別目標組織中起重要作用的基因集，二是在基因集檢驗中因變量與被分析組織的正常活動和組織特異性過程顯著關聯。由于權重較大的基因集更能反映正常活動和組織特異性，在這兩種情況下，p值加權能更好地提高統計功效。

表3 組織特異性基因集檢驗結果

表4列出了使用上述方法進行wFDR分析骨骼肌相對于二型糖尿病產生的兩個基因集發現。

表4 在骨骼肌中與二型糖尿病有關的通路

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

基因集權重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 人類組織的功能表征

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

3 結 論

使用人類蛋白質圖譜的組織特異性基因活性信息和過濾后的分子簽名數據庫中的基因集生成了組織特異性的基因集權重。為了避免權重較大的基因集影響結果的準確性，將分子簽名數據中的基因集合進行過濾。利用組織特異性基因集權重對p值進行加權并用BH方法提供FDR控制，進行wFDR分析，在不同組織的三種疾病下加權的FDR一共有64個發現，而未加權的FDR一共有34個發現。在p值小于0.2的假設下，使用組織特異性基因集權重對p值進行加權能產生更多的發現。使用組織特異性基因集權重可以提高基因集檢驗的統計功效，并且也可以提高在高維基因組數據的實驗中識別生物學上有關的基因集關聯的信息。