Limma: how to make a table of logFC of certain genes (I need to define logFC of 62 genes specific for the cells of interest)
1
0
Entering edit mode
Элина • 0
@0a782e5a
Last seen 22 months ago
Russia

I'm using limma package to count differential gene expression. I need to make a table of logFC of genes specific for my cells of interest (there are 62 of them). What should I write to make such a table (where should I place the names of the genes)? My code at the moment: Code should be placed in three backticks as shown below

gset <- getGEO("GSE24742", GSEMatrix =TRUE, AnnotGPL=TRUE)
if (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

fvarLabels(gset) <- make.names(fvarLabels(gset))
make.names(fvarLabels(gset))

gsms <- "010101010101010101010101"
sml <- strsplit(gsms, split="")[[1]]
sml

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)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprs(gset) <- log2(ex) }

gs <- factor(sml)
gs
groups <- make.names(c("control","test"))
levels(gs) <- groups
gset$group <- gs
design <- model.matrix(~group + 0, gset)
colnames(design) <- levels(gs)
View(design)
fit <- lmFit(gset, design)  

cts <- paste(groups[1], groups[2], sep="-")
cont.matrix <- makeContrasts(test-control, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
cont.matrix

fit2 <- eBayes(fit2, 0.01)
tT1 <- topTable(fit2, adjust="fdr", sort.by="B", number=50000)
tT1 <- subset(tT1, select=c("ID","adj.P.Val","P.Value","t","B","logFC","Gene.symbol","Gene.title"))
write.table(tT1, file=stdout(), row.names=F, sep="\t")
View(tT1)

I guess the function "topTable" must be modified somehow or replaced with something else. I will be very grateful for any help.

DifferentialExpression limma R TissueMicroarrayData • 822 views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia

Making a table for specified genes

It's easy to subset the toptable for any desired set of genes:

i <- tT1$Gene.symbol %in% MyGeneList
tT1.MyGenes <- tT1[i, ]

The samples are paired

However your current analysis is incorrect because it ignores the pairing of the samples. There is also a strong mean-variance trend. A correct analysis would be:

gset <- getGEO("GSE24742", GSEMatrix =TRUE, AnnotGPL=TRUE)[[1]]
exprs(gset) <- logs(exprs(gset))
Patient <- gl(12,2,24)
Treatment <- gl(2,1,24)
plotMDS(gset,label=paste(Patient,Treatment,sep="."))
design <- model.matrix(~Treatment+Patient)
fit <- lmFit(gset,design)
fit <- eBayes(fit,trend=TRUE)
tab <- topTable(fit,coef="Treatment2",n=Inf)

Beware that this dataset shows no significant results

You should be aware that this dataset shows no statistically significant responses to treatment. The original publication used ordinary t-tests and failed to adjust the p-values for multiple testing. In my opinion, claiming any results from this dataset would be scientifically unsound.

ADD COMMENT

Login before adding your answer.

Traffic: 396 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6