保姆式GEO數(shù)據(jù)挖掘演示
配對GEO數(shù)據(jù)分析
GEO數(shù)據(jù)挖掘系列文
TCGAbiolinks做差異分析
我的一個paired DEGs 的例子
rm(list=ls())
suppressMessages(library(limma))
suppressMessages(library("plyr"))
data_dir<-"your_path"
setwd(data_dir)
#####download data
if(T){
suppressMessages(library(GEOquery))
gset <- getGEO('GSE76161', destdir=".",
AnnotGPL = T,
getGPL = T)
save(gset,file=paste0(human_dir,'/GSE76161_eSet.Rdata'))
}
load(paste0(human_dir,'/GSE76161_eSet.Rdata'))
###Get gene expression and phenotype data
b = gset[[1]]
exprSet=exprs(b)
###get GPL information and clean the data
GPL = b@featureData@data
colnames(GPL)
ids = GPL[,c(1,3)]
ids = ids[ids[,2] != '',]
a<-strsplit(as.character(ids[,2]), "http:///")
tmp <- mapply(cbind, ids[,1],a)
df <- ldply (tmp, data.frame)
ID2gene = df
colnames(ID2gene) = c("id", "gene")
exprSet = exprSet[rownames(exprSet) %in% ID2gene[,1],]
ID2gene = ID2gene[match(rownames(exprSet), ID2gene[,1]),]
dim(ID2gene)
dim(exprSet)
tail(sort(table(ID2gene[,2])), n = 12L)
{
MAX = by(exprSet, ID2gene[,2],
function(x) rownames(x)[ which.max(rowMeans(x))])
MAX = as.character(MAX)
exprSet = exprSet[rownames(exprSet) %in% MAX,]
rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]
}
exprSet[1:5,1:5]
pdata = pData(b)
colnames(pdata)
save(exprSet, pdata, file = paste0(human_dir,'/finalSet.Rdata'))
###DEG analysis
#build the matrix
SibShip <- factor(pdata$characteristics_ch1.3,levels=c("patient_id: 3","patient_id: 4","patient_id: 5"))
Treat<-factor(pdata[,47],levels=c("vehicle (0.1% DMSO)","1 uM OCA (INT-747)"))
design <- model.matrix(~0 + (SibShip+Treat))
rownames(design) = colnames(exprSet)
fit1 <- lmFit(exprSet, design)
#fit2 <- contrasts.fit(fit1, design)
fit2 <- eBayes(fit1)
nrDEG = topTable( fit2,
adjust.method="BH",
coef = "Treat1 uM OCA (INT-747)",
n = nrow(fit2)
)
nrDEG<-nrDEG[order(nrDEG$P.Value),]
write.csv(nrDEG,paste0(data_dir,"/GSE76161_human_DEGs_INT-747_vs_DMSO.csv"))
