开始你的创作:GEO2R使用
library(GEOquery)
library(limma)
library(edgeR)
library(tidyverse)
(function(){
GSE ="GSE126094"
message("开始分析: ",GSE," Hsa_circ_101555 functions as a competing endogenous RNA of miR-597-5p to promote colorectal cancer progression [Circular RNAs]")
if(!file.exists(paste0("GEO-COAD/",GSE,".rds"))){
gset <- getGEO(GSE,AnnotGPL=T,destdir = "GEO-COAD")
saveRDS(gset,paste0("GEO-COAD/",GSE,".rds"))
}
gset <- readRDS(paste0("GEO-COAD/",GSE,".rds"))
gset <- gset[[1]]
gset_expr <- exprs(gset)%>%
as.data.frame()
message(GSE," 使用GPL平台 ",gset@annotation)
#boxplot(log2(gset_expr),outline=FALSE)
#title(main=GSE)
gset_pd <- pData(gset)%>%
dplyr::select(sample_id=geo_accession,sample_name=title)%>%
mutate(group=case_when(grepl("normal",sample_name)~"Normal",
T ~"Tumor"))
gset_annot <- gset@featureData@data
gset_expr_annot <- gset_expr %>%
mutate(symbol=gset_annot$circRNA)%>%
rownames_to_column("id")%>%
dplyr::select(-id) %>%
filter(symbol!="") %>%
stats::aggregate(. ~ symbol, data =., FUN = mean)%>%
tibble::column_to_rownames("symbol")
saveRDS(gset_expr_annot,paste0("GEO-COAD/gset_expr_annot_circRNA_",GSE,".rda"))
saveRDS(gset_pd,paste0("GEO-COAD/metadata_circRNA_",GSE,".rda"))
})()
(function(){
GSE ="GSE138589"
message("开始分析: ",GSE,"circRNA Microarray for human colorectal adenocarcinoma and normal mucosa")
if(!file.exists(paste0("GEO-COAD/",GSE,".rds"))){
gset <- getGEO(GSE,AnnotGPL=T,destdir = "GEO-COAD")
saveRDS(gset,paste0("GEO-COAD/",GSE,".rds"))
}
gset <- readRDS(paste0("GEO-COAD/",GSE,".rds"))
gset <- gset[[1]]
gset_expr <- exprs(gset)%>%
as.data.frame()
message(GSE," 使用GPL平台 ",gset@annotation)
#boxplot(gset_expr,outline=FALSE)
#title(main=GSE)
gset_pd <- pData(gset)%>%
dplyr::select(sample_id=geo_accession,sample_name=title)%>%
mutate(group=case_when(grepl("normal",sample_name)~"Normal",
T ~"Tumor"))
gset_annot <- gset@featureData@data
gset_expr_annot <- gset_expr %>%
mutate(symbol=gset_annot$circRNA)%>%
rownames_to_column("id")%>%
dplyr::select(-id) %>%
filter(symbol!="") %>%
stats::aggregate(. ~ symbol, data =., FUN = mean)%>%
tibble::column_to_rownames("symbol")
saveRDS(gset_expr_annot,paste0("GEO-COAD/gset_expr_annot_circRNA_",GSE,".rda"))
saveRDS(gset_pd,paste0("GEO-COAD/metadata_circRNA_",GSE,".rda"))
})
(function(){
GSE ="GSE142837"
message("开始分析: ",GSE," circRNA expression data from human colorectal cancer tissues and the correponding non-cancerous tissues")
if(!file.exists(paste0("GEO-COAD/",GSE,".rds"))){
gset <- getGEO(GSE,AnnotGPL=T,destdir = "GEO-COAD")
saveRDS(gset,paste0("GEO-COAD/",GSE,".rds"))
}
gset <- readRDS(paste0("GEO-COAD/",GSE,".rds"))
gset <- gset[[1]]
gset_expr <- exprs(gset)%>%
as.data.frame()
message(GSE," 使用GPL平台 ",gset@annotation)
#boxplot(log2(gset_expr),outline=FALSE)
#title(main=GSE)
gset_pd <- pData(gset)%>%
dplyr::select(sample_id=geo_accession,sample_name=title)%>%
mutate(group=case_when(grepl("non-cancerous",sample_name)~"Normal",
T ~"Tumor"))
gset_annot <- gset@featureData@data
gset_expr_annot <- gset_expr %>%
mutate(symbol=gset_annot$circRNA)%>%
rownames_to_column("id")%>%
dplyr::select(-id) %>%
filter(symbol!="") %>%
stats::aggregate(. ~ symbol, data =., FUN = mean)%>%
tibble::column_to_rownames("symbol")
saveRDS(gset_expr_annot,paste0("GEO-COAD/gset_expr_annot_circRNA_",GSE,".rda"))
saveRDS(gset_pd,paste0("GEO-COAD/metadata_circRNA_",GSE,".rda"))
})()
library("maptools") # point labels without overlaps
library(umap)
(function(){
GEOS <- c("GSE126094","GSE142837")
Normalized <- c(F,F)
for(i in 1:length(GEOS)){
GEO <- GEOS[i]
message("Normalized:",Normalized[i])
message("开始分析:",GEO)
expr <- readRDS(paste0("GEO-COAD/gset_expr_annot_circRNA_",GEO,".rda"))
metadata <- readRDS(paste0("GEO-COAD/metadata_circRNA_",GEO,".rda"))
message("分组信息列名与表达矩阵的行名相同",identical(colnames(expr),rownames(metadata)))
group <- factor(metadata$group, levels = c("Normal", "Tumor"))
## 检测是否需要log2转换
qx <- as.numeric(quantile(expr, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) {
message("##############",GEO,"log2 转换")
expr[which(expr <= 0)] <- NaN
expr <- log2(expr)
}
## 质控
if(T){
par(mfrow=c(1,3),mar=c(7,2,1,6))
ord <- order(group)
boxplot(expr[,ord], boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
plotDensities(expr, group=group, main=GEO, legend ="topright")
ump <- umap(t(expr), n_neighbors = 9, random_state = 123)
plot(ump$layout, main="UMAP plot, nbrs=9", xlab="", ylab="", col=group, pch=20, cex=1.5)
legend("topright", legend=levels(group), pch=20,
col=1:nlevels(group), title="Group", pt.cex=1.5,inset=c(-0.6,0))
pointLabel(ump$layout, labels = rownames(ump$layout), method="SANN", cex=0.6)
}
design <- model.matrix( ~ 0 + group)
colnames(design) <- levels(group)
message("Tumor-Normal")
contrast.matrix <- makeContrasts("Tumor-Normal", levels = design)
## 是否强制标准化
if(!Normalized[i]){
message(GEO," 进行Normalized")
if(F){
dge <- DGEList(counts=expr)
dge <- calcNormFactors(dge)
v <- voom(dge,design,plot=TRUE, normalize="quantile")
}
if(F){
expr <- expr[complete.cases(expr), ]
v <- vooma(log2(expr), design, plot=T)
expr <- normalizeBetweenArrays(expr) # normalize data
par(mfrow=c(1,1),mar=c(7,3,2,2))
boxplot(log2(v[['E']]), boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
}
expr<- normalizeBetweenArrays(expr)
par(mfrow=c(1,3),mar=c(7,2,1,6))
ord <- order(group)
boxplot(expr[,ord], boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
plotDensities(expr, group=group, main=GEO, legend ="topright")
ump <- umap(t(expr), n_neighbors = 9, random_state = 123)
plot(ump$layout, main="UMAP plot, nbrs=9", xlab="", ylab="", col=group, pch=20, cex=1.5)
legend("topright", legend=levels(group), pch=20,
col=1:nlevels(group), title="Group", pt.cex=1.5, inset=c(-0.6,0))
pointLabel(ump$layout, labels = rownames(ump$layout), method="SANN", cex=0.6)
}
fit <- lmFit(expr, design)%>%
contrasts.fit(contrast.matrix)%>%
eBayes()
#contrast <- paste0(rev(levels(group)), collapse = "-")
res <- topTable(fit, coef = 1, n = Inf)
saveRDS(res,file = paste0("GEO-COAD/limma_",GEO,".rda"))
}
if(F){
GEO <- GEOS[1]
expr <- readRDS(paste0("GEO-COAD/gset_expr_annot_circRNA_",GEO,".rda"))
metadata <- readRDS(paste0("GEO-COAD/metadata_circRNA_",GEO,".rda"))
message("分组信息列名与表达矩阵的行名相同",identical(colnames(expr),rownames(metadata)))
group <- factor(metadata$group, levels = c("Normal", "Tumor"))
par(mfrow=c(1,3),mar=c(7,2,1,5))
ord <- order(group)
boxplot(log2(expr[,ord]), boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
plotDensities(log2(expr), group=group, main=GEO, legend ="topright")
#par(mfrow=c(1,1),mar=c(2,2,2,2))
ump <- umap(t(log2(expr)), n_neighbors = 9, random_state = 123)
#par(mar=c(3,3,2,6), xpd=TRUE)
plot(ump$layout, main="UMAP plot, nbrs=9", xlab="", ylab="", col=group, pch=20, cex=1.5)
legend("topright", legend=levels(group), pch=20,
col=1:nlevels(group), title="Group", pt.cex=1.5, inset=c(-0.3,0))
pointLabel(ump$layout, labels = rownames(ump$layout), method="SANN", cex=0.6)
design <- model.matrix( ~ 0 + group)
colnames(design) <- levels(group)
contrast <- paste0(rev(levels(group)), collapse = "-")
message(contrast)
contrast.matrix <- makeContrasts(contrast, levels = design)
if(F){
dge <- DGEList(counts=log2(expr))
dge <- calcNormFactors(dge)
v <- voom(dge,design,plot=TRUE, normalize="quantile")
v <- vooma(expr, design, plot=T)
par(mfrow=c(1,1),mar=c(7,3,2,2))
boxplot(log2(v[['E']]), boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
}
expr<- normalizeBetweenArrays(log2(expr))
boxplot(expr, boxwex=0.6, notch=T, main=GEO, outline=FALSE, las=2, col=group[ord])
fit <- lmFit(expr, design)%>%
contrasts.fit(contrast.matrix)%>%
eBayes()
res <- topTable(fit, coef = 1, n = Inf)
saveRDS(res,file = paste0("GEO-COAD/limma_",GEO,".rda"))
}
})()