Rstudio和隨機叢林在高維全基因組學數據分析中的應用*

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

Rstudio和隨機叢林在高維全基因組學數據分析中的應用*

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

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

全基因組關聯分析 SNPs Rstudio 隨機森林 隨機叢林

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

方 法

1.Rstudio

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

(1) Rstudio Server的搭建

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

(2)Rstudio Sweave

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

2.分析方法

(1)Cochron-Armitage檢驗

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

(2)隨機叢林

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

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

實際數據

(1)數據來源

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

(2)數據預處理與質量控制

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

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

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

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

結 果

1.Rstudio Server的啟動與登陸

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

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

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

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

3.WTCCC高血壓GWAS實際數據分析

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

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

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

通過隨機叢林篩選出10次平均重要性評分穩定排序在前20位的SNPs位點(見表2)共匹配到12個基因(位點與基因匹配采用了課題組利用dbSNP數據庫自行更新的基因映射信息，若1個SNP位點在某個基因內部或距離較近則認為此SNP位點與該基因相關)，TUSC3基因和LCORL基因各有2個位點，LOC100422352基因和MIR378D1基因各有4個位點，其他8個基因各有1個位點；對比Cochron-Armitage檢驗的P值排序結果可見，有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值有統計學意義(P<5×10-7)。

(3)Rstudio Sweave結果輸出

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

討 論

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

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

3.本文在預分析時在不同γ2值條件下，提取的標簽SNP對RJ分析結果的影響較大，且為盡量保證不遺漏一些效應相對較弱的位點，因此事先并未進行prune SNP，直接利用封裝后的RJ程序對經質量控制之后的360302個SNPs位點進行分析。5 CPUs并行運算平均僅耗時為10.37h/每次，且隨機叢林所篩選的SNPs位點對疾病具有很好的分類能力，在病例對照標簽打亂前后位點的重要性評分差距明顯，重要性評分排序前20的位點中Cochron-Armitage檢驗的P值也有10個排進前11位，尤其是重要性評分排序前8的位點兩種方法結果一致性很好，說明隨機叢林能夠很好地篩選出主效應相對較強的位點；另一方面， Cochron-Armitage檢驗的P值還有10個位點排名較為靠后， 提示隨機叢林能發現一些 Cochron-Armitage檢驗顯示的效應相對較弱的位點，分析其原因可能在于位點間存在交互作用。盡管如此，基于隨機森林的特征篩選方法在結果解釋方面仍面臨著挑戰，其篩選交互作用的能力易受位點主效應的影響，在高維數據分析過程中更多的時候可能僅是扮演了一個初篩的角色[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：打造一個可重復的統計研究流程.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]武曉巖，李康.基因表達數據判別分析的隨機森林方法.中國衛生統計，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.哈爾濱醫科大學衛生統計學教研室(150081)

2.哈爾濱醫科大學附屬第一醫院心內科

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