顧星博 李 昂 溫 琪 焦辛妮 劉 艷△

Rstudio和隨機叢林在高維全基因組學數(shù)據(jù)分析中的應(yīng)用*

顧星博1，2李 昂1溫 琪1焦辛妮1劉 艷1△

目的 結(jié)合Rstudio和Random Jungle兩款軟件的優(yōu)勢，通過遠程登錄簡便、快捷地實現(xiàn)全基因組學數(shù)據(jù)的分析。方法 在服務(wù)器端搭建Rstudio Server，封裝隨機叢林R程序并對英國威康信托病例對照協(xié)會(WTCCC)高血壓真實數(shù)據(jù)進行分析，通過Rstudio Sweave動態(tài)生成分析結(jié)果。結(jié)果 在客戶端即可通過網(wǎng)絡(luò)瀏覽器登陸服務(wù)器端搭建好的Rstudio Server；數(shù)據(jù)量較大時，封裝后的隨機叢林程序在Rstudio Server中運算速度方面相對于R randomForest包優(yōu)勢明顯。2001名高血壓病例及3004名對照的22條染色體共有490032個SNP位點，根據(jù)真實標簽通過10次隨機叢林過程篩選出穩(wěn)定排序在前20位的SNPs位點，其Cochron-Armitage檢驗的P值也有10個排進前11位。結(jié)論 全基因組學數(shù)據(jù)的處理會消耗很多的時間和內(nèi)存，普通計算機根本無法承受，Random Jungle軟件的命令行運行方式又不易于數(shù)據(jù)處理、算法組合或嵌套及結(jié)果的再現(xiàn)、可視化，在服務(wù)器端搭建Rstudio Server并結(jié)合Random Jungle的分析策略可有效地應(yīng)用于全基因組學數(shù)據(jù)分析，簡化分析過程、提高分析速度和效率、方便實現(xiàn)分析結(jié)果的動態(tài)輸出及再現(xiàn)。

全基因組關(guān)聯(lián)分析 SNPs Rstudio 隨機森林 隨機叢林

隨著人類基因組測序技術(shù)的飛速發(fā)展，全基因組關(guān)聯(lián)研究(genome-wide association studies，GWAS)方法已在全球范圍內(nèi)廣泛應(yīng)用于復(fù)雜疾病關(guān)聯(lián)的候選區(qū)域篩選當中[1]。GWAS所涉及到的(single-nucleotide polymorphisms，SNPs)位點數(shù)目可達數(shù)十萬甚至上百萬，其交互作用中位點的組合更為復(fù)雜，對如此海量數(shù)據(jù)進行處理分析時將面臨一個重要的問題即運算負擔問題，這決定著能否方便、快捷、有效地對數(shù)據(jù)進行處理分析。本研究旨在提出在服務(wù)器端搭建Rstudio Server并與隨機叢林相結(jié)合的一種新的高維數(shù)據(jù)分析策略，使運用R和隨機森林進行GWAS數(shù)據(jù)分析成為可能，同時將分析結(jié)果通過Rstudio Sweave進行動態(tài)呈現(xiàn)。

方 法

1.Rstudio

Rstudio是基于C++開發(fā)的一種R語言的集成開發(fā)環(huán)境(IDE)，相比于R的工作界面其功能更加人性化，除可應(yīng)用于各種操作系統(tǒng)的桌面版本之外，其特有的Server版本可在Linux服務(wù)器上進行配置，用戶通過網(wǎng)頁瀏覽器即可使用Rstudio，使得在高性能計算機上遠程多人同時在線對大型數(shù)據(jù)進行運算及程序開發(fā)成為可能，而且Rstudio同時支持Sweave和Tex，可使輸出結(jié)果動態(tài)呈現(xiàn)并以極為靈活的方式創(chuàng)建出符合出版質(zhì)量的分析報告形式[2]。

(1) Rstudio Server的搭建

最新版本的Rstudio Server v0.98需要Debian 6/Ubuntu 10.04或者Red Hat/CentOS 5.4以上操作系統(tǒng)環(huán)境，在目標服務(wù)器的系統(tǒng)中添加CRAN庫后通過命令安裝最新版本R[3]，然后再通過命令安裝Rstudio Server，并對Rstudio Server進行配置，即在/etc/rstudio目錄下創(chuàng)建并編輯rserver.conf.txt和rsession.conf.txt，這兩個文本文件將分別包含與服務(wù)器運行相關(guān)和與R運行相關(guān)的配置信息(主要的配置信息見表1)，配置完成之后管理員即可在服務(wù)器終端使用通過命令添加用戶。

(2)Rstudio Sweave

Sweave是R與LaTeX的結(jié)合，旨在打造一個可重復(fù)的統(tǒng)計分析流程[4]。LaTeX主要用于文章、書籍和研究報告的排版，是Leslie Lamport基于TeX排版語言編寫的一組宏代碼，擁有更為規(guī)范的命令。Sweave最早只是一個R程序包，需通過包中函數(shù)Sweave對已通過文本編輯器編譯過的Rnw格式的noweb文件進行處理，進而生成tex格式的LaTeX文件，最后再通過LaTeX文本編譯器進行編譯生成PDF格式的報告。目前Rstudio的開發(fā)團隊已將Sweave整個“移植”到了Rstudio當中，只需事先在操作系統(tǒng)中安裝好LaTeX組件即可對源代碼直接進行編譯。

2.分析方法

(1)Cochron-Armitage檢驗

Cochron-Armitage檢驗是由William Cochran和Peter Armitage提出的[5]，用于分析一個二分類變量和另外一個二分或多分類變量之間的關(guān)聯(lián)性，是對Pearson卡方檢驗的修改，Cochron-Armitage檢驗是最為常用的GWAS單變量分析方法之一，本研究使用運算效率高的PLINK軟件進行Cochron-Armitage檢驗的計算。

(2)隨機叢林

隨機森林(random forset，RF)[6]作為一門經(jīng)典的數(shù)據(jù)挖掘方法，已在多項研究中應(yīng)用于基因組學數(shù)據(jù)分析[7]，在篩選交互作用時RF比傳統(tǒng)的單位點分析方法更具有優(yōu)勢。但以往的RF分析軟件(如R中randomForest包)源代碼主要由Fortran編寫，會受到運算時間和內(nèi)存的限制。隨機叢林(random jungle，RJ)是近年來發(fā)展起來的RF多核并行運算軟件包，可在Linux和Windows環(huán)境下使用。RJ在繼承了最為常用的R中randomForset包功能的同時對原有隨機森林方法進行了擴展[8]，如除可計算常用Gini重要性評分、Permutation重要性評分之外還可計算Meng式重要性評分及條件重要性評分；在變量篩選方面補充了可靈活設(shè)置的后退法并且對構(gòu)建森林時CART的類型進行了補充。

RJ源代碼由C++編寫，且其程序框架也不同于其他RF分析軟件，在處理高維數(shù)據(jù)時運算時間及內(nèi)存使用上的優(yōu)勢極為突出。但其操作需通過命令行的方式，在數(shù)據(jù)處理時輸入和輸出都是獨立的文本格式文件并不方便后續(xù)的整理、分析及可視化，因此本研究將RJ在Rstudio Server中以函數(shù)的形式進行了封裝，從而提供了一個中間層使RJ和工作站中已搭建的Rstudio Server之間無縫接入。封裝完成之后在Rstudio Server中調(diào)用RJ的數(shù)據(jù)分析結(jié)果可與直接使用RJ軟件結(jié)果一致。

實際數(shù)據(jù)

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

本研究數(shù)據(jù)來源于課題組申請下載的英國威康信托病例對照協(xié)會(WTCCC)于2007年公布的英國白種人全基因組高血壓SNPs數(shù)據(jù)[9]，該數(shù)據(jù)包括2001名高血壓病例及3004名對照(其中1958年英國出生隊列1504人，血站獻血者1500人)，22條染色體共490032個SNP位點。

(2)數(shù)據(jù)預(yù)處理與質(zhì)量控制

首先，利用Python語言編寫的腳本程序?qū)⒃璫hiamo格式數(shù)據(jù)轉(zhuǎn)換成tped/tfam格式數(shù)據(jù)，并參考數(shù)據(jù)下載時WTCCC所提供的剔除的調(diào)查者編號對不滿足GWAS人群分層現(xiàn)象的樣本給予剔除(WTCCC采用MDS plot、主成分分析等方法剔除樣本)；其次，使用PLINKv1.07[10]對數(shù)據(jù)進行質(zhì)量控制，主要剔除標準除參照WTCCC質(zhì)量控制標準之外還將最小等位基因頻率(minor allele frequency，MAF)小于0.05的位點進行了剔除，得到1952個病例和2938個對照個體共360302個SNP位點的ped/map格式數(shù)據(jù)；最后，利用PLINK軟件將ped/map格式數(shù)據(jù)轉(zhuǎn)換成可在R中進行分析的raw格式數(shù)據(jù)。

(3)分析流程與程序運行背景

第一步，在工作站Rstudio Server中通過RJ篩選重要性評分位于前100位的位點，重復(fù)10次；第二步，利用permutation過程將病例對照標簽打亂，再使用RJ計算打亂后的前100位位點的重要性評分，并與打亂前結(jié)果(即第一步的分析結(jié)果)進行比較，同樣重復(fù)10次；第三步，將打亂前結(jié)果中重要性評分均穩(wěn)定靠前的SNP位點與單變量Cochron-Armitage檢驗的結(jié)果進行比較。

本研究的工作站環(huán)境為Linux操作系統(tǒng)，Intel至強5600，CPU主頻2.93GHz，12核64G內(nèi)存。通過預(yù)實驗，考慮隨機叢林OOB分類錯誤率以及運算負擔，設(shè)定RJ運行的主要參數(shù)mtry=36000(0.1P，P為全部SNPs位點數(shù))，ntree為5000，重要性評分計算方式選擇與R中randomForest包Permutation重要性評分相同的Liaw Score，標簽置換前后20次運算均采用5 CPUs并行。輸出結(jié)果利用自行編寫的Rnw格式的Sweave腳本程序在Rstudio Server中進行編譯直接生成PDF文檔報告，并方便后期的查看和驗證。

結(jié) 果

1.Rstudio Server的啟動與登陸

管理員通過運行rstudio-server start命令啟動Rstudio Server，若在運行過程中對相應(yīng)的配置文件進行了更改，則需運行rstudio-server restart命令重啟服務(wù)，以使服務(wù)器重新讀取更新的配置信息。Rstudio Server服務(wù)運行后，用戶在客戶端的網(wǎng)絡(luò)瀏覽器中通過IP即可登陸，在登陸頁面中輸入帳號和密碼后出現(xiàn)使用界面，使用界面中包含有R控制臺、文件管理器、軟件包管理器、圖形輸出窗口等，即客戶端能調(diào)用遠程目標服務(wù)器所提供的R的各項功能。Rstudio Server可實現(xiàn)多人同時登陸，每個用戶會有自己的獨立目錄，他人不能訪問，可保證隱私性；另外，也可設(shè)立每個用戶都能訪問的公共目錄，方便對他人代碼或程序進行調(diào)用，從而對代碼和程序進行更好的管理與共享。

2.RJ和RF在Rstudio Server中運行速度比較

圖1為在Rstudio Server中利用模擬的SNPs數(shù)據(jù)集比較RJ和RF(R中randomForest包)在串行(非并行)和多核并行時的運算時間(RF采用randomForest包和foreach包進行并行)。模擬數(shù)據(jù)集源自于R中MDR包train數(shù)據(jù)集，該數(shù)據(jù)集包含446名病例和554名對照共5000個SNPs位點和一個二分類表型變量，接近常見疾病的真實數(shù)據(jù)且方便在R中直接利用命令加載。實驗過程中分別保留該數(shù)據(jù)集100、500、1000、5000個位點，建立5000棵樹在單核和10核環(huán)境下運行。

可以看出，在單核環(huán)境下RJ運算速度明顯優(yōu)于RF且隨著SNPs位點數(shù)量的增多，提速效果越明顯(5000位點時RJ耗時172秒，比randomForest快26.34倍)。在10核并行環(huán)境下，當SNPs位點數(shù)為100時RJ耗時48秒比foreach包并行后的randomForest慢4.38倍；SNPs位點數(shù)量為500和1000時，兩種方法運算時間相近(RJ略慢于RF)；而當SNPs位點數(shù)量達到5000時，RJ僅耗時74秒，比RF快42.13倍。

3.WTCCC高血壓GWAS實際數(shù)據(jù)分析

(1) 隨機叢林篩選致病SNPs位點

標簽打亂前后20次5 CPUs并行運算，共耗時207.39h，平均每次耗時10.37h。RJ篩選的重要性評分位于前100位位點的分布圖(見圖2)中，橫坐標為重要性評分值，縱坐標為概率密度值。由圖中可以看出打亂病例對照標簽前位點的重要性評分相對較高，而打亂病例對照標簽后位點的重要性評分都處于很低的水平，提示利用RJ篩選出的SNPs位點對病例與對照有很好的分類能力。

(2)隨機叢林結(jié)果與Cochron-Armitage檢驗篩選結(jié)果的比較

通過隨機叢林篩選出10次平均重要性評分穩(wěn)定排序在前20位的SNPs位點(見表2)共匹配到12個基因(位點與基因匹配采用了課題組利用dbSNP數(shù)據(jù)庫自行更新的基因映射信息，若1個SNP位點在某個基因內(nèi)部或距離較近則認為此SNP位點與該基因相關(guān))，TUSC3基因和LCORL基因各有2個位點，LOC100422352基因和MIR378D1基因各有4個位點，其他8個基因各有1個位點；對比Cochron-Armitage檢驗的P值排序結(jié)果可見，有9個位點的Cochron-Armitage檢驗的P值排序也在前10位，且TUSC3基因、LCORL基因、LOC100422352基因和MIR378D1基因各有1個位點(rs7837736、rs1553460、rs10843660、rs16837871和其他8個基因的6個位點(rs1528356、rs4131463、rs825148、rs17116117、rs6840033、rs10499044)Cochron-Armitage檢驗P值有統(tǒng)計學意義(P<5×10-7)。

(3)Rstudio Sweave結(jié)果輸出

在Rstudio中通過運行已編寫完成的Sweave腳本程序，將數(shù)據(jù)分析代碼嵌入到PDF文檔報告即在報告中動態(tài)插入R輸出結(jié)果的同時，也就將結(jié)果的計算過程明確地記錄下來，方便過后進行查看和驗證；另外，還可以免去繁瑣的復(fù)制、粘貼和重新格式化結(jié)果的過程。

討 論

1.高維組學數(shù)據(jù)的數(shù)據(jù)量龐大，其數(shù)據(jù)分析過程也更為復(fù)雜，通常需要團隊多人多階段協(xié)作完成。R在處理大數(shù)據(jù)時表現(xiàn)的并不盡如人意[11]，而Rstudio Server的搭建可為團隊協(xié)作提供一個統(tǒng)一的程序運行及文件管理平臺，在高性能計算機上可有效地避免因本地設(shè)備條件限制、機器性能不足等可能導(dǎo)致的大型統(tǒng)計運算工作難以實現(xiàn)的問題。在工作站搭建Rstudio Server之后，客戶端用戶不再需要一系列R的安裝與配置，只需在有互聯(lián)網(wǎng)的前提下使用網(wǎng)頁瀏覽器即可在熟悉的軟件環(huán)境中進行數(shù)據(jù)的處理分析，使用戶在擺脫了硬件架構(gòu)、系統(tǒng)平臺束縛的同時，也可實現(xiàn)數(shù)據(jù)分析與計算的高效化、快捷化。

2.在Rstudio Server中封裝后的隨機叢林R程序可有效地結(jié)合Rstudio與Random Jungle兩個軟件的優(yōu)勢，在數(shù)據(jù)量較大時運算速度方面的優(yōu)勢極為明顯，且方便數(shù)據(jù)處理、程序調(diào)用、結(jié)果輸出及可視化，亦可令整個流程便于驗證和重現(xiàn)；即便是隨機叢林與其他算法進行更為復(fù)雜的組合或嵌套時，也可在R中方便地得以實現(xiàn)。而Sweave則可令整個分析流程動態(tài)實現(xiàn)，即只需在Rstudio中對事先編寫好的相應(yīng)程序進行改動，即對相應(yīng)部分的代碼稍作調(diào)整之后點擊編譯即可自動生成分析報告，從而使整個分析流程更加高效、便捷、可重復(fù)。

3.本文在預(yù)分析時在不同γ2值條件下，提取的標簽SNP對RJ分析結(jié)果的影響較大，且為盡量保證不遺漏一些效應(yīng)相對較弱的位點，因此事先并未進行prune SNP，直接利用封裝后的RJ程序?qū)?jīng)質(zhì)量控制之后的360302個SNPs位點進行分析。5 CPUs并行運算平均僅耗時為10.37h/每次，且隨機叢林所篩選的SNPs位點對疾病具有很好的分類能力，在病例對照標簽打亂前后位點的重要性評分差距明顯，重要性評分排序前20的位點中Cochron-Armitage檢驗的P值也有10個排進前11位，尤其是重要性評分排序前8的位點兩種方法結(jié)果一致性很好，說明隨機叢林能夠很好地篩選出主效應(yīng)相對較強的位點；另一方面， Cochron-Armitage檢驗的P值還有10個位點排名較為靠后， 提示隨機叢林能發(fā)現(xiàn)一些 Cochron-Armitage檢驗顯示的效應(yīng)相對較弱的位點，分析其原因可能在于位點間存在交互作用。盡管如此，基于隨機森林的特征篩選方法在結(jié)果解釋方面仍面臨著挑戰(zhàn)，其篩選交互作用的能力易受位點主效應(yīng)的影響，在高維數(shù)據(jù)分析過程中更多的時候可能僅是扮演了一個初篩的角色[12]。

[1] Moore JH,Asselbergs FW,Williams SM.Bioinformatics challenges for genome-wide association studies.Bioinformatics,2010,26(4):445-455.

[2]Rstudio.http://www.rstudio.com/products/rstudio/.

[3]Rstudio Server.http://www.rstudio.com/products/rstudio/download-server/.

[4]謝益輝.Sweave：打造一個可重復(fù)的統(tǒng)計研究流程.http:// cos.name/2010/11/reproducible-research-in-statistics/

[5]Armitage P.Tests for linear trends in proportions and frenquencies.Biometrics,1955,11(3):375-386.

[6]Breiman L.Random forests.Machine learning,2001,45(1):5-32.

[7]武曉巖，李康.基因表達數(shù)據(jù)判別分析的隨機森林方法.中國衛(wèi)生統(tǒng)計，2006，23(6)：491-494.

[8] Schwarz DF,Konig IR,Ziegler A.On safari to Random Jungle:a fast implementation of Random Forests for high-dimensional data.Bioinformatics,2010,29(14):1752-1758.

[9]WellcomeTrust Case Control,Consortium.Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls.Nature,2007,447(7145):661-678.

[10]Purcell S,B Neale,K Todd-Brown,et al.PLINK:a tool set for whole-genomeassociation and population- based linkage analyses.The American Journal of Human Genetics,2007,81(3):559-575.

[11]VincentZoonekynd.Statistics with R(2007).

[12] Lin HY,Chen YA,Tsai,et al.TRM:a powerful two-stage machine learning approach for identifying SNP-SNP interactions.Annals of human genetics,2012,76(1):53-62.

(責任編輯：鄧 妍)

The Application of Rstudio and Random Jungle in High-dimensional Genome-wide Data

Gu Xingbo,Li Ang,Wen Qi,et al.

(School of Public Health，Harbin Medical University(150081),Heilongjiang)

Objective Combining the advantages of Rstudio and Random Jungle software to realize the efficient analysis of genome-wide data via remote login.Methods To set up the Rstudio Server,use the packaged Random Jungle R program to analyze the Wellcome Trust Case Control Consortium hypertension GWAS data and generate the dynamic analysis results.Results Client side people can login the Rstudio Server through a web browser; When the amount of data is large,the packaged Random Jungle grogram speed obviously faster than the R randomForest package.2001 cases of HT affected individuals and 3004 controls of 22 chromosomes with total 490032 SNPs.The top 20 SNPs identified by RJ which use actual label on 10 times that also have 10 in top 11 of Cochron-Armitage testPvalues.Conclusion Genome-wide data processing will cost a lot of time and memory and it incompetent to afford by ordinary computers completely.Random Jungle Software’s command-line operation mode may not be easy to data processing,algorithm combinations,reproduction of results and visualizations.The provided analysis of the strategy which combines the Rstudio and Random Jungle software can be well applied in GWAS data analysis.This strategy can simplify the analysis process,improving the speed and efficiency of the data analysis.Furthermore,it is convenient to achieve dynamic output and reproduction.

GWAS;SNPs;Rstudio;Random Forest;Random Jungle

*國家自然科學基金(30972537,81172741)

1.哈爾濱醫(yī)科大學衛(wèi)生統(tǒng)計學教研室(150081)

2.哈爾濱醫(yī)科大學附屬第一醫(yī)院心內(nèi)科

△通信作者：劉艷，E-mail:liuyan@ems.hrbmu.edu.cn；