Question: How to compare DGE between different groups and create a volcano plot in R based on a GEO2R dataset?
1
ishackm • 100 wrote:
Hi all,
I have the following GEO2R dataset which has 4 different groups. I would like to perform DGE on these groups please on R and create volcano plots to compare the DGE. For example, G0 vs G1.
This is my code so far in R:
> # Version info: R 3.2.3, Biobase 2.30.0, GEOquery 2.40.0, limma 3.26.8
> # R scripts generated Tue Mar 26 09:31:37 EDT 2019
>
> ################################################################
> # Differential expression analysis with limma library(Biobase) library(GEOquery) library(limma)
>
> # load series and platform data from GEO
>
> gset <- getGEO("GSE40595", GSEMatrix =TRUE, AnnotGPL=TRUE) if
> (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx
> <- 1 gset <- gset[[idx]]
>
> # make proper column names to match toptable fvarLabels(gset) <- make.names(fvarLabels(gset))
>
> # group names for all samples gsms <- paste0("00000000111111111111111111111111111111122222233333",
> "333333333333333333333333333") sml <- c() for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }
>
> # log2 transform ex <- exprs(gset) qx <- as.numeric(quantile(ex, 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) || (qx[2] > 0 && qx[2] < 1 && qx[4]
> > 1 && qx[4] < 2) if (LogC) { ex[which(ex <= 0)] <- NaN exprs(gset) <- log2(ex) }
>
> # set up the data and proceed with analysis sml <- paste("G", sml, sep="") # set group names fl <- as.factor(sml) gset$description <-
> fl design <- model.matrix(~ description + 0, gset) colnames(design) <-
> levels(fl) fit <- lmFit(gset, design) cont.matrix <-
> makeContrasts(G3-G0, G1-G0, G2-G1, G3-G2, levels=design) fit2 <-
> contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2, 0.01) tT <-
> topTable(fit2, adjust="fdr", sort.by="B", number=250)
>
> tT <- subset(tT,
> select=c("ID","adj.P.Val","P.Value","F","Gene.symbol","Gene.title"))
> write.table(tT, file=stdout(), row.names=F, sep="\t")
If you could guide me on how I could edit this code, that would be really helpful.
Many Thanks,
Ishack
