使用R語言metaMA程序包實現(xiàn)基因表達(dá)譜數(shù)據(jù)的Meta分析

李賢菁，鄧巧玲，李嘉興，張晨曦，闕雅婷，翁鴻，黃靜宇，李勝,3,4

基因芯片技術(shù)又稱微陣列技術(shù)[1]，因其具有操作簡單，信息量大，響應(yīng)迅速等特點，在醫(yī)學(xué)和生物學(xué)等領(lǐng)域得到廣泛運用。DNA芯片屬于基因芯片的一種類型，它是將已知序列的核苷酸作為探針，高密度排布在固相載體上。熒光標(biāo)記的cDNA與載體上的已知序列雜交，用相應(yīng)計算機軟件檢測和處理熒光信號強度，從而獲得基因特征性序列或者特異性表達(dá)的相關(guān)生物信息。

隨著分子生物學(xué)和基因組學(xué)的深入研究，基因芯片技術(shù)不斷得到發(fā)展和完善。芯片尺寸增大，種類與數(shù)據(jù)類型也愈加豐富。而如何能夠有效的處理和分析來自基因芯片的大量數(shù)據(jù)，也是研究者們關(guān)注的焦點。R是一款功能強大的統(tǒng)計分析軟件，擁有完整的數(shù)據(jù)處理、統(tǒng)計和繪圖系統(tǒng)以及操作環(huán)境。在生物信息方面，研究者通過開發(fā)不同R包，以滿足分析不同類型基因芯片數(shù)據(jù)的需要。

基于R語言的metaMA包由Guillemette Marot開發(fā)[2]。不同研究中通過比較正常組織和癌組織不同基因表達(dá)量的差異而篩選出的差異表達(dá)基因不完全相同，使用metaMA包進(jìn)行meta分析篩選出差異表達(dá)基因的方法基于更大樣本量，因此可增加結(jié)論的統(tǒng)計功效。此文中涉及metaMA分析流程的R語言程序包主要有GEOquery[3]、org.Hs.eg.db[4]、metaMA[2]、VennDiagram[5]、db1[6]、db2[7]、db3[8]、annaffy[9]。其中，GEOquery用于基因組數(shù)據(jù)的獲取，org.Hs.eg.db用于實現(xiàn)數(shù)據(jù)轉(zhuǎn)換，metaMA用于找出差異表達(dá)的基因，VennDiagram用于繪制多組研究中差異基因的維恩圖，db1、db2和db3用于對3例基因組數(shù)據(jù)中的差異基因進(jìn)行注釋，annaffy用于生成已識別探針的生物信息。本文將以3例口腔鱗狀細(xì)胞癌數(shù)據(jù)集為例，說明如何對公開的不同研究數(shù)據(jù)進(jìn)行統(tǒng)一轉(zhuǎn)化，繼而運用metaMA包找出差異表達(dá)的基因，繪制維恩圖，并通過對差異基因注釋得出結(jié)論。

1 數(shù)據(jù)準(zhǔn)備

1.1 軟件R的安裝及基因組數(shù)據(jù)讀入 從官方網(wǎng)站（https://www.r-project.org）下載最新版R軟件R-3.4.2并進(jìn)行安裝。獲取GEO數(shù)據(jù)，可以使用R語言的GEOquery包。安裝、加載GEOquery包，以及下載相應(yīng)的口腔鱗狀細(xì)胞癌數(shù)據(jù)集。具體代碼如下：

source("https://bioconductor.org/biocLite.R")

biocLite("GEOquery")

library(GEOquery)

data1<- getGEO('GSE9844')

data2 <- getGEO('GSE3524')

data3 <- getGEO('GSE13601')

數(shù)據(jù)集可以直接從GEO（https://www.ncbi.nlm.nih.gov/geo/）下載txt文件，用read.table（）讀入R，或者通過Excel整理出如表二所示的基因表達(dá)矩陣，用read.csv（）讀入R。

1.2 示例數(shù)據(jù)集 以3例口腔鱗狀細(xì)胞癌數(shù)據(jù)集為例，說明metaMA的分析流程。3例數(shù)據(jù)集如表1所示。

表1 示例3例口腔鱗狀細(xì)胞癌數(shù)據(jù)集

1.3 芯片預(yù)處理 基因芯片預(yù)處理過程，即是根據(jù)

熒光信號生成的數(shù)據(jù)將其轉(zhuǎn)變成基因組數(shù)據(jù)的過程。質(zhì)量控制和數(shù)據(jù)處理對基因芯片分析的有效性至關(guān)重要。通過質(zhì)量控制，可以剔除不合格芯片，保留合格芯片做后續(xù)的數(shù)據(jù)處理分析。數(shù)據(jù)處理包含以下步驟，分別是背景校正，標(biāo)準(zhǔn)化，匯總。背景校正用于去除背景噪聲的影響。標(biāo)準(zhǔn)化用于消除樣品制備等非實驗因素對結(jié)果的影響。最后經(jīng)過匯總得到探針組水平的數(shù)據(jù)。目前較為流行的預(yù)處理方法有dChip,MAS5,RMA。前兩種的背景校正方法均為mas，后一種的背景校正方法為rma。

1.4 根據(jù)基因組數(shù)據(jù)繪制箱線圖 在使用metaMA程序包進(jìn)行Meta分析之前，需要查看以上三個數(shù)據(jù)集基因表達(dá)水平的分布情況，這里繪制各組基因表達(dá)數(shù)據(jù)的箱線圖，能夠直觀比較預(yù)處理效果具體代碼如下：

par(mfrow=c(2,2))

boxplot(data.frame(exprs(data1[[1]])),outline=FALSE)

boxplot(data.frame(exprs(data2[[1]])),outline=FALSE)

boxplot(data.frame(exprs(data3[[1]])),outline=FALSE)

圖1 箱線圖

箱線圖中，橫坐標(biāo)為樣本，縱坐標(biāo)為表達(dá)值，比較表達(dá)值從而獲得結(jié)果變異程度的信息。從獲得箱線圖（圖1A、1B和1C）中可看出，數(shù)據(jù)3（data3）的值遠(yuǎn)大于數(shù)據(jù)1（data1）和數(shù)據(jù)2（data2），變異相對較大。數(shù)據(jù)1[10]、數(shù)據(jù)2[11]與數(shù)據(jù)3[12]的數(shù)據(jù)處理方法不完全相同。數(shù)據(jù)1使用RMA方法對芯片數(shù)據(jù)做預(yù)處理。數(shù)據(jù)2使用Mas5.0軟件對Affymetrix數(shù)據(jù)進(jìn)行初步分析，進(jìn)一步標(biāo)準(zhǔn)化數(shù)據(jù)后，輸出到應(yīng)用程序Genespring 6.0（Silicon Genetics，Redwood City，CA）。數(shù)據(jù)3使用Mas5.0方法的預(yù)處理，未進(jìn)行對數(shù)轉(zhuǎn)換進(jìn)行標(biāo)準(zhǔn)化，這里我們對data3數(shù)據(jù)集進(jìn)行2為底的對數(shù)轉(zhuǎn)換，同時我們再次繪制標(biāo)準(zhǔn)化后的data3基因表達(dá)水平分布圖（圖1D），具體代碼如下：

exprs(data3[[1]])<-log2(exprs(data3[[1]]))

boxplot(data.frame(exprs(data3[[1]])),outline=FALSE)1.5 明確樣本分組 如前所述，本研究使用的3例口腔癌數(shù)據(jù)集包含正常組織和癌組織兩類樣本，在使用metaMA進(jìn)行meta分析前，需要對以上兩類樣本進(jìn)行分組、編碼。我們將正常口腔組織編碼為1，口腔癌組織編碼為0。具體代碼如下：

c1<-as.numeric(pData(data1[[1]])["source_name_ch1"]=="Oral Tongue Squamous Cell

Carcinoma")

c2<-as.numeric(apply(pData(data2[[1]])["descript ion"], 1,toupper)=="SERIES OF 16 TUMORS")

c3<-as.numeric(pData(data3[[1]])["source_name_ch1"]=="Tumor")

classes<-list(c1,c2,c3)

1.6 將基因表達(dá)矩陣轉(zhuǎn)化為統(tǒng)一的Unigene形式如表1所示，這三個口腔癌基因表達(dá)譜分別采用不同平臺的基因芯片進(jìn)行檢測，且三個數(shù)據(jù)集的探針集是不同的，因此需要將其轉(zhuǎn)換成形式統(tǒng)一的基因編號。函數(shù)probe2unigene返回這3組獨立研究中的探針?biāo)鶎?yīng)的所有UNIGENE。 函數(shù)unigene2probe將UNIGENE與Entrez Gene ID相對應(yīng)。 要實現(xiàn)這一轉(zhuǎn)換，需要加載Bioconductor中的org.Hs.eg.db包。具體代碼如下：

source("http://bioconductor.org/biocLite.R")

biocLite("org.Hs.eg.db")

require("org.Hs.eg.db")

x <- org.Hs.egUNIGENE

mapped_genes <- mappedkeys(x)

link<- as.list(x[mapped_genes])

probe2unigene<-function(expset) {

#construction of the map probe->unigene

probes=rownames(exprs(expset))

gene_id=fData(expset)[probes,"ENTREZ_GENE_ID"]

unigene=link[gene_id]

names(unigene)<-probes

probe_unigene=unigene

}

unigene2probe<-function(map) {

suppressWarnings(x <- cbind(unlist(map),names(map)))

unigene_probe=split(x[,2], x[,1])

}

需要注意的是，無法與UNIGENE相匹配的探針將不被用于Meta分析。因為一個unigene可能會對應(yīng)多個探針，選擇一個恰當(dāng)?shù)姆椒▉砗喜ⅹ毩⒀芯康奶结樦陵P(guān)重要。下面代碼中，根據(jù)已被證實有效的參數(shù)mergemeth的默認(rèn)值，convert2metaMA選擇通過計算它們的平均值來總結(jié)對應(yīng)于相同UNIGENE的探針的值。convert2metaMA函數(shù)為每個研究構(gòu)建矩陣，矩陣中只包含所有研究里共有的UNIGENE對應(yīng)的行（基于UNIGENE轉(zhuǎn)換）。這些矩陣存儲在函數(shù)返回的值的子集列表中。convert2metaMA還會返回conv.unigene，作用是在進(jìn)行meta分析后能夠返回原始探針的名稱。具體代碼如下：

convert2metaMA<-function (listStudies,mergeme th=mean) {

if (!(class(listStudies) %in% c("list"))) {

stop("listStudies must be a list") }

conv_unigene=lapply(listStudies,

FUN=function(x) unigene2probe(probe2unigene(x)))

id=lapply(conv_unigene,names)

inter=Reduce(intersect,id)

if(length(inter)<=0){stop("no common genes")}

print(paste(length(inter),"genes in common"))

esets=lapply(1:length(listStudies),FUN=function(i){

l=lapply(conv_unigene[[i]][inter],

FUN=function(x) exprs(listStudies[[i]])[x,,drop=T RUE])

esetsgr=t(sapply(l,FUN=function(ll) if(is.null(dim(ll))){ll}

else{apply(ll,2,mergemeth)}))

esetsgr

})

return(list(esets=esets,conv.unigene=conv_unigene))

}

conv<-convert2metaMA(list(data1[[1]],data2[[1]],data3[[1]]))

輸出結(jié)果為

[1]"7745 genes in common"#說明此三個基因組中存在相交情況的基因數(shù)目位為7745。

esets<-conv$esets

conv_unigene<-conv$conv.unigene

上圖中，表2為data1的基因表達(dá)矩陣部分，其中，第一行是樣本編號，第一列為探針編號。表3為合并后的基因表達(dá)矩陣部分，第一列為探針對應(yīng)的UNIGENE，數(shù)目即為三組研究中相交的基因數(shù)目。

2 應(yīng)用metaMA包進(jìn)行分析

加載metaMA包時需要同時用到包limma和包SMVar。包limma屬于R的基礎(chǔ)包之一。而包SMVar屬于Bioconductor，如未安裝，則需要安裝之后才能順利加載metaMA包。

source("http://bioconductor.org/biocLite.R")

biocLite("SMVar")

install.packages(“metaMA”)

library(metaMA)

res<-pvalcombination(esets=esets,classes=classes)

輸出結(jié)果如表4

通過length函數(shù)，我們可以查看本次Meta分析所得出的差異基因的數(shù)目，具體代碼如下：

length(res$Meta)

輸出結(jié)果為

[1]2481#說明3組數(shù)據(jù)中，差異表達(dá)基因數(shù)目為2481。

表2 數(shù)據(jù)1的表達(dá)譜

表3 轉(zhuǎn)化為UNIGENE后的表達(dá)譜

表4 Meta分析后相關(guān)與差異基因相關(guān)的統(tǒng)計結(jié)果

此外我們還可以查看經(jīng)Meta分析后的差異基因結(jié)果，具體代碼如下：

Hs.Meta<-rownames(esets[[1]])[res$Meta]

3 繪制維恩圖

通過繪制維恩圖，可以直觀比較獨立研究或者M(jìn)eta分析中差異表達(dá)基因的數(shù)目。安裝VennDiagram包可以實現(xiàn)該功能。代碼如下：

install.packages("VennDiagram")

library(VennDiagram)

venn.plot<-venn.diagram(x=list(study1=res$stud y1,study2=res$study2,

study3=res$study3, meta=res$Meta), filename =NULL, col = "black",

fill = c("blue", "red", "purple","green"),margin=0.05, alpha = 0.6)

jpeg("venn_jpeg.jpg")

grid.draw(venn.plot)

dev.off()

維恩圖中，藍(lán)色，紫色，紅色區(qū)域分別表示在三個獨立研究中篩選出的差異基因，綠色表示經(jīng)Meta分析篩選出的2481個差異基因。由維恩圖知，經(jīng)由metaMA包篩選出而不能在獨立研究中被識別的差異基因數(shù)目為219（圖2）。

圖2 經(jīng)分析后的差異表達(dá)基因的維恩圖

4 進(jìn)一步篩選差異基因

經(jīng)metaMA分析獲得的2481個差異基因，需進(jìn)一步進(jìn)行篩選得出關(guān)鍵核心基因?？蓮腡CGA數(shù)據(jù)庫中下載頭頸部腺癌的基因表達(dá)矩陣，根據(jù)臨床信息從中篩選出口腔鱗狀細(xì)胞癌數(shù)據(jù)納入研究。本次研究共篩選出211個癌組織樣本和27個正常組織對照。對其進(jìn)行差異基因表達(dá)分析，獲得P值＜0.01的差異基因后，取這些基因與上述經(jīng)metaMA分析獲得的2481個基因的交集。將交集中的基因結(jié)合TCGA數(shù)據(jù)庫提供的臨床信息進(jìn)行生存分析，以確定關(guān)鍵核心基因。

5 注釋

5.1 獲取芯片類型 在找到差異表達(dá)基因后，用getanndb函數(shù)獲得相應(yīng)的芯片類型。

getanndb<-function(expset) {

gpl_name=annotation(expset)

gpl=getGEO(gpl_name)

title=Meta(gpl)$title

db=paste(gsub("[-|_]","",

tolower(strsplit(title,"[[]|[]]")

[[1]][[2]])),"db",sep=".")

return(db)

}

db1<-getanndb(data1[[1]])

db1

#[1]"hgu133plus2.db"

db2<-getanndb(data2[[1]])

db2

#[1]"hgu133a.db"

db3<-getanndb(data3[[1]])

db3

#[1]"hgu95av2.db"

5.2 使用注釋包對差異表達(dá)基因進(jìn)行注釋

#source("http://bioconductor.org/biocLite.R")

#biocLite(db1)

#biocLite(db2)

#biocLite(db3)

library(db1,character.only=TRUE)

library(db2,character.only=TRUE)

library(db3,character.only=TRUE)

5.3 獲取注釋文件 在下面代碼中，origId.Meta可返回各獨立研究中基因的原始IDs，原始IDs與進(jìn)行Meta分析的共同差異表達(dá)基因的Hs.Meta IDs相對應(yīng)。由于幾個原始的探針匹配相同的UNIGENE，origId.Meta可能會比Hs.Meta包含更多的ID。annaffy包中的函數(shù)aafTableAnn和saveHTML可用于保存當(dāng)前工作目錄中的一個命名為“annotation.html”的 HTML文件，該文件包含本次分析的差異基因及其相應(yīng)的生物信息。

origId.Meta<-lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta])))

source("http://bioconductor.org/biocLite.R")

biocLite("annaffy")

library(annaffy)

annlist<-lapply(1:length(origId.Meta),FUN=function(i)

aafTableAnn(origId.Meta[[i]],chip=get(paste0("db", i))))

annot<-do.call(rbind,annlist)

saveHTML(annot,file="annotation.html",title="Responder genes")

6 討論

基因芯片技術(shù)發(fā)展迅速，被用于生物學(xué)、醫(yī)學(xué)等多個領(lǐng)域研究。Meta分析是基于多項研究成果的薈萃分析。目前用于獲取Meta分析所需的基因表達(dá)數(shù)據(jù)的主要平臺包括GEO（https://www.ncbi.nlm.nih.gov/geo/）及ArrayExpress（http://www.ebi.ac.uk/arrayexpress/）等。按照常規(guī)研究思路，通過meta分析，篩選出與疾病相關(guān)的差異表達(dá)基因后，需要對差異表達(dá)的基因進(jìn)行預(yù)測和功能分析。探究這些基因的潛在的生物學(xué)功能，有助于分析解釋某些疾病的發(fā)生機理，為治療提供新思路[13]。

嚴(yán)格的文章篩選和有效的質(zhì)量評價是獲得可靠結(jié)論的必要條件，此外還需要識別和控制發(fā)表偏倚。在輸入Meta分析的研究數(shù)據(jù)時，需要將不同平臺、不同類型芯片之間的差異納入考慮范圍。有時還需選擇合適方法處理尚未進(jìn)行預(yù)處理的數(shù)據(jù)。本文的3組基因組數(shù)據(jù)定量分析了口腔鱗狀細(xì)胞癌相關(guān)基因的基因表達(dá)情況，并沒有重點討論標(biāo)準(zhǔn)化方法，僅是對數(shù)據(jù)3進(jìn)行以2為底的對數(shù)轉(zhuǎn)換。在進(jìn)行嚴(yán)謹(jǐn)?shù)幕蛐酒治鲞^程中，需要對數(shù)據(jù)處理方法進(jìn)行評估和優(yōu)化。