####======================getElist()========================####
## getElist()通过给定GSE号,自动下载标准格式的soft文件,提取其中的矩阵、注释和target信息,最后形成eset类列表。此函数依赖于GEOquery,可以获得详尽的表型、注释、处理流程内容,十分强大。不过,如果是外地下载soft文件,应该改为GSEXXXXX.soft.gz的名称才能正确识别。如果想保存eset文件为rda,可以选择savefile = T。如果有个性化GSE.soft文件保存空间,可以自定义gse.space属性。
## 升级情况
# 2018-09-26:修复了有多种格式的pdata(比如肿瘤和正常标本列数、列名不致)无法合并的bug(from GSE17154)
# 2018-09-26:修复了有GSE中原作者可能存在的GSM名字错误导致eset无法正常生成的bug(from GSE17154)
#' @export
getElist <- function(gse.names,
gse.space=NULL,
savefile = F){
##加载必要的包
library(pacman)
p_load("stringr","GEOquery","DT","Biobase")
####生成私人空间
if(is.null(gse.space)){
public.workspace = "E:/RCloud/database/DataDownload/GEOqueryDownload" # 这个因人而异
gse.space = str_c(public.workspace,"/",gse.names)
dir.create(gse.space,recursive = T)
} else {
dir.create(gse.space,recursive = T)
}
##下载GSE数据
gse <- getGEO(gse.names,GSEMatrix=F,destdir = gse.space)
##判断平台信息,选择其中一个平台
gsmplatforms <- lapply(GSMList(gse),function(x){Meta(x)$platform_id})
platforms <- NULL;
for (i in 1:length(gsmplatforms)) {
g.i <- gsmplatforms[[i]]
platforms <- c(platforms,g.i)
}
platforms.type = unique(platforms) #全部相同
print(str_c("platform information:",platforms.type))
##构建Elist
eset2 <- list()
#i=1
for (i in 1:length(platforms.type)) {
#选择其中一个平台
gsmlist <- Filter(function(gsm){Meta(gsm)$platform_id==platforms.type[i]},GSMList(gse));length(gsmlist)
GSM.names <- names(gsmlist) #记录名
## featureData(GPL平台)
print(str_c("开始第",i,"个featureData构建..."))
gpl <- GPLList(gse)[[i]]
gpl1 <- Table(gpl)
l1 <- gpl1$ID %in% NA #有时候这一列有NA值,影响rownames注释。删除。
gpl1 <- gpl1[!l1,]
rownames(gpl1) <- gpl1$ID
gpl1$ID <- rownames(gpl1) #有时ID是数字,这给后面的注释带来不良影响。因此统一ID类为字符型变量。
feature <- as(gpl1,"AnnotatedDataFrame")
print(str_c("完成第",i,"个featureData构建!"))
## exprs
print(str_c("开始第",i,"个exprs构建..."))
probesets <- gpl1$ID
data.matrix <- do.call('cbind',lapply(gsmlist,function(x)
{tab <- Table(x)
mymatch <- match(probesets,tab$ID_REF)
return(tab$VALUE[mymatch])
}))
data.matrix <- apply(data.matrix,2,function(x){as.numeric(as.character(x))})
rownames(data.matrix) <- probesets
colnames(data.matrix) <- names(gsmlist)
#View(data.matrix[1:10,1:10])
data.matrix <- data.matrix[,GSM.names]
print(str_c("完成第",i,"个exprs构建!"))
## phenoData
print(str_c("开始第",i,"个phenoData构建..."))
pdata <- GSMList(gse)
## 制定统一列名:有时这里有列名经常写错(比如多打一个字母)。这里统一用第一个array的信息作用列名。
cols <- list();len <- NULL;
for(a in 1:length(pdata)){
pdata.i <- pdata[a]
pdata.i2 <- pdata.i[[1]]@header[["characteristics_ch1"]]
#vt <- pdata.i2;z="node: 0"
extra.pdata <- function(vt,n,split="[:] "){
vt <- as.character(vt)
a <- NULL
for(z in vt){
z.i <- unlist(strsplit(z,split))
a <- c(a,z.i[n])
}
return(a)
}
col.a <- extra.pdata(pdata.i2,1)
len <- c(len,length(col.a))
cols <- c(cols,list(col.a))
}
len.types <- unique(len)
len.ps <- match(len.types,len)
list.colnames <- cols[len.ps];names(list.colnames) <- len.types
## 提取多个芯片的信息。
df.pdata <- data.frame()#j=2
for (j in 1:length(pdata)) {
pdata.i <- pdata[j]
pdata.i2 <- pdata.i[[1]]@header[["characteristics_ch1"]]
df.j <- data.frame(pattern = extra.pdata(pdata.i2,1),
value = extra.pdata(pdata.i2,2))
colnames(df.j)[2] <- names(pdata.i)
df.j <- as.data.frame(t(df.j))
colnames(df.j) <- list.colnames[[match(ncol(df.j),len.types)]]
if(ncol(df.j) == 1){
#即只有一个变量,此时数据框变化要多费周折
coln <- colnames(df.j)
df.j <- df.j[2,]
df.j <- as.character(df.j)
df.j <- as.data.frame(df.j)
colnames(df.j) <- coln;rownames(df.j) <- names(pdata.i)
} else {
df.j <- df.j[-1,]
}
library(plyr)
df.pdata <- rbind.fill(df.pdata,df.j)
}
rownames(df.pdata) <- GSM.names
pheno <- as(df.pdata,"AnnotatedDataFrame")
print(str_c("完成第",i,"个phenoData构建!"))
## protocalData
print(str_c("开始第",i,"个protocalData构建..."))
library(plyr)
proData <- NULL
for(z in 1:length(gsmlist)){
gsm1 <- gsmlist[[z]]
gsm.df <- data.frame(v1=gsm1@dataTable@columns[["Description"]]);rownames(gsm.df) <- as.character(gsm1@dataTable@columns[["Column"]])
#提取header的信息。如果某一项有2个及以上值,合并
header <- gsm1@header;length(header)
gsm.df2 <- NULL
for(a in 1:length(header)){
header.i <- header[[a]];len=length(header.i)
header.i <- paste0(header.i,collapse = "_")
header.i2 <- data.frame(v1=header.i,v2=len);rownames(header.i2)[1] <- names(header[a])
gsm.df2 <- rbind(gsm.df2,header.i2)
#print(paste0("a",a))
}
gsm.df2 <- gsm.df2[order(gsm.df2$v2,decreasing = T),]
gsm.df2 <- subset(gsm.df2,select = "v1")
#进一步处理header信息
gsm.df4 <- rbind(gsm.df,gsm.df2)
colnames(gsm.df4)[1] <- names(gsmlist[z])
gsm.df5 <- as.data.frame(t(gsm.df4))
#合并数据
proData <- rbind.fill(proData,gsm.df5)
#print(paste0("z",z))
}
rownames(proData) <- GSM.names
proData <- as(proData,"AnnotatedDataFrame")
print(str_c("完成第",i,"个protocalData构建!"))
## 生成ExpressionSet object
eset <- new('ExpressionSet',
exprs=data.matrix,
phenoData=pheno,
featureData = feature,
protocolData = proData)
print(str_c("完成第",i,"个eset构建!"))
eset2 <- c(eset2,list(eset))
names(eset2)[i] <- names(GPLList(gse))[i]
}
colnames(data.matrix)
rownames(df.pdata)
##输出结果
if(length(eset2) == 1){
eset3 <- eset2[[1]]
print(str_c(gse.names,"为单注释平台数据。"))
} else {
eset3 <- eset2
print(str_c(gse.names,"为多注释平台数据。"))
}
## End
if(savefile == F){
tuneR::play(music)
return(eset3)
} else {
eset <- eset3
save(eset,file = str_c(gse.space,"/",gse.names,"_eset.rda"))
tuneR::play(music)
return(eset3)
}
}
## 例子
# gse.names = "GSE11121"
# View(exprs(eset)[1:10,]) #矩阵信息
# View(fData(eset)[1:10,]) #注释信息
# View(pData(eset)) #表型信息
# View(protocolData(eset)@data) #流程信息
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.