But how to ignore this thing and run the whole code? that is my question
I am trying to execute this GEO2R script in RStudio version 4.0.
# Differential expression analysis with limma
library(Biobase)
library(GEOquery)
library(limma)
# load series and platform data from GEO
gset <- getGEO("GSE89333", 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 <- "0011"
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
# Create the study design
design <- model.matrix(~ description + 0, gset)
colnames(design) <- levels(fl)
# Fit the linear model on the study's data
fit <- lmFit(gset, design)
# Make individual contrasts
cont.matrix <- makeContrasts(G1-G0, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
# Applying the empirical Bayes method to the fitted values
# Acts as an extra normalisation step and aims to bring the different probe-wise variances to common values
fit2 <- eBayes(fit2, 0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="B", number=250)
tT <- subset(tT, select=c("ID","Gene.symbol","logFC"))
write.table(tT, file=stdout(), row.names=F, sep="\t")
But it is giving me error in gsms (group samples) when I am running the script with different microarray platforms. How can i exclude this step? Please help me
1 answer
This part of the code:
# group names for all samples
gsms <- "0011"
sml <- c()
for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }
And then later:
# set up the data and proceed with analysis
sml <- paste("G", sml, sep="") # set group names
fl <- as.factor(sml)
gset$description <- fl
the above code results in:
gset$description
[1] G0 G0 G1 G1 Levels: G0 G1
This is a round-about way of creating the factor describing the treatments - check the Samples section at the GEO Accession Viewer:
GSM2366230 scrambled_day3_rep1
GSM2366231 scrambled_day3_rep2
GSM2366232 siCHRNA5_day3_rep1
GSM2366233 siCHRNA5_day3_rep2
Which means, these two snippets are specific to this one experiment (GSE89333), and you will need to modify the code to properly describe other experiments.
And my answer is you can't ignore, and you will have to examine experiment by experiment, because you have to understand the experimental design in order to code proper contrasts and perform a meaningful analysis.
Log in to answer this question.
Please someone help me regarding this. It really urgent
I answered your question bellow. But the contributors to the forum are volunteers and answer here on their free time, it really isn't nice to try to hurry people into answering a question, specially on a Sunday.
Please take a minute to read this more detailed comment on this issue: C: cuffdiff error. IT hangs