ArrayExpress

常用函数

1 ae2bioc

功能：将MAGE-TAB文件从raw data转换为bioconductor对象。

返回值： AffyBatch, ExpressionSet 或 NChannelSet。assayData储存表达值，phenoData储存sdrf，featureData储存arf，experimentData储存idf。

ae2bioc(mageFiles, dataCols = NULL, drop = TRUE)

2 ArrayExpress

功能：从ArrayExpress数据库下载raw data。

返回值：AffyBatch、 ExpressionSet 或 NChannelSet。

ArrayExpress(accession, path = tempdir(), save = FALSE, dataCols = NULL, drop = TRUE)

3 getAE

功能：从ArrayExpress数据库下载并提取MAGE-TAB文件。

返回值：一个含有下载并提取的文件名称的list。

getAE(accession, path = getwd(), type = "full", extract = TRUE, local = FALSE, sourcedir = path) 

4 getcolproc

功能：提取processed MAGE-TAB文件的列名。

返回值：一个含有processed MAGE-TAB文件的列名的list。

getcolproc(files, path)

5 getcolraw

功能：提取raw MAGE-TAB文件的列名。

返回值：一个含有raw MAGE-TAB文件的列名的list。

getcolraw(path, rawfiles)

6 procset

功能：将本地的processed MAGE-TAB文件转化成ExpressionSet。

返回值：ExpressionSet。

procset(files, procol)

7 queryAE

功能：搜索ArrayExpress数据库并下载搜索结果。

返回值：包含搜索结果信息的dataframe。

queryAE(keywords = NULL, species = NULL) 

该函数在我的电脑上运行出错。在github上找到源代码，将其中的query = try(download.file(qr, queryfilename, mode="wb"))改成query = try(download.file(qr, queryfilename, mode="wb", method="auto"))，手动输入代码，创建函数，再次执行命令。成功。

修改后的代码：

getelt = function(x, node, element){
elt = sapply(1:length(xmlRoot(x)), function(i){
  if(length(grep(node,names(xmlRoot(x)[[i]])))  != 0)
    unlist(xmlElementsByTagName(xmlRoot(x)[[i]], node))[names(unlist(xmlElementsByTagName(xmlRoot(x)[[i]], node)))== element]  else "NA" })

sizeelt = sapply(1:length(xmlRoot(x)), function(i){
  if(length(grep(node,names(xmlRoot(x)[[i]])))  != 0)
    length(unlist(xmlElementsByTagName(xmlRoot(x)[[i]], node))[names(unlist(xmlElementsByTagName(xmlRoot(x)[[i]], node)))== element])  else "NA" })

if(inherits(elt, "list") || (inherits(elt, "character") && max(sizeelt[sizeelt!="NA"]) == 1))
{
    elt2 = lapply(elt, function(i) if(length(i) == 0) "NA" else i)

    elt3 = unlist(lapply(elt2, function(i) do.call("paste",c(as.list(i),sep=" | "))))
    } 
if(inherits(elt, "matrix") && max(sizeelt) > 1)
   elt3 = do.call("paste",c(as.list(elt),sep=" | "))

names(elt3) = NULL
return(elt3)
}

geteltmulti = function(x, node, element1, element2){
elt = sapply(1:length(xmlRoot(x)), function(i){
 if(length(grep(node,names(xmlRoot(x)[[i]])))  != 0){
      lapply(1:length(xmlElementsByTagName(xmlRoot(x)[[i]], node)), function(j) {
    extr = unlist(xmlElementsByTagName(xmlRoot(x)[[i]], node)[[j]])

    e1 = extr[names(extr)== element1]
    e2 = extr[names(extr)== element2]

    paste(e1, e2, sep="=")
   })} else "NA" })

if(inherits(elt, "list"))
{
    elt2 = lapply(elt, function(i) if(length(i) == 0) "NA" else i)
    elt3 = lapply(elt2, function(i) unlist(i))
    elt4 = unlist(lapply(elt3, function(i) do.call("paste",c(as.list(i),sep=" | "))))
} else elt4 = do.call("paste",c(as.list(do.call("paste",c(as.list(elt),sep=" | ")), sep="||")))

names(elt4) = NULL
return(elt4)
}

queryAE = function(keywords = NULL, species = NULL){

if(is.null(keywords) && is.null(species))
    stop("No keywords or species specified")

baseURL = "https://www.ebi.ac.uk/arrayexpress/xml/v2/experiments";

qr = paste(baseURL,"?keywords=",keywords,"&species=",species,sep="")
qr=URLencode(qr)
queryfilename = paste("query",keywords,species,".xml",sep="")
query = try(download.file(qr, queryfilename, mode="wb", method="auto"))

x = xmlTreeParse(queryfilename)

ID = sapply(1:length(xmlRoot(x)), function(i) unlist(xmlElementsByTagName(xmlRoot(x)[[i]], "accession"))["accession.children.text.value"])
names(ID) = NULL
x2 = xmlTreeParse(queryfilename, useInternalNodes = TRUE)
ra = getNodeSet(x2,"/experiments//raw[@count]")
##Raw = sapply(ra, function(r) xmlGetAttr(r, "count"))
Rawids = sapply(ra, function(r) xmlGetAttr(r, "name"))
Rawids = gsub(".raw.*.zip","",Rawids)
##names(Raw) = Rawids
Raw = rep(NA,length(ID))
names(Raw) = ID
Raw[which(ID %in% Rawids)] = 'yes'
Raw[which(!ID %in% Rawids)] = 'no'


pr = getNodeSet(x2,"/experiments//fgem[@count]")
##Processed = sapply(pr, function(p) xmlGetAttr(p, "count"))
Procids = sapply(pr, function(r) xmlGetAttr(r, "name"))
Procids = gsub(".processed.*.zip","",Procids)
##names(Processed) = Procids
Processed = rep(NA,length(ID))
names(Processed) = ID
Processed[which(ID %in% Procids)] = 'yes'
Processed[which(!ID %in% Procids)] = 'no'

date = getelt(x, node = "releasedate",  element = "releasedate.children.text.value")

pmid = getelt(x, node = "bibliography",
  element = "bibliography.children.accession.children.text.value")

spec = getelt(x, node = "species",
  element = "species.children.text.value")

experimentdesign = getelt(x, node = "experimentdesign",
  element = "experimentdesign.children.text.value")   

experimentalfactor = geteltmulti(x, node = "experimentalfactor",
  element1 = "children.name.children.text.value",
  element2 = "children.value.children.text.value")

xmlparsed = data.frame(ID = ID, Raw = Raw[ID], Processed = Processed[ID], ReleaseDate = date, PubmedID = pmid, Species = spec, ExperimentDesign = experimentdesign, ExperimentFactors = experimentalfactor)
return(xmlparsed)
}

example：

1 搜索ArrayExpress数据库

library("ArrayExpress")
sets <- queryAE(keywords = "pneumonia", species = "homo+sapiens")

2 下载数据

rawset <- ArrayExpress("E-MEXP-1422")

3 下载需要提供colnames的数据

eset <- try(ArrayExpress("E-MEXP-1870"))

报错，存在重复的行名。但范例说要手动输入列名。

eset = ArrayExpress("E-MEXP-1870",
 dataCols=list(R="ScanArray Express:F650 Mean",
 G="ScanArray Express:F550 Mean",
 Rb="ScanArray Express:B650 Mean",
 Gb="ScanArray Express:B550 Mean"))

仍旧报错。同样的错误。

4 下载ArrayExpress Datasets

mexp1422 = getAE("E-MEXP-1422", type = "full")

5 转换本地 raw MAGE-TAB 文件为R对象

使用上一步getAE的返回值。

rawset= ae2bioc(mageFiles = mexp1422)

6 转换本地 processed MAGE-TAB 文件为R对象

6.1 确定要提取的列名

cn = getcolproc(mexp1422)
show(cn)

6.2 转换本地 processed MAGE-TAB 文件为R对象

proset = procset(mexp1422, cn[2])