Question

Query regarding limma

0

Entering edit mode

Roopa Subbaiaih ▴ 130

@roopa-subbaiaih-5490

Last seen 11.0 years ago

United States

Hi All,

I was doing microarray analysis where I compare healthy with diseased samples. The script which I use is

getwd()
setwd(dir="/CRSP 406-11/DEMOS/GSE14905-a")
ls()
#-----------------------------------------------#
library(affy)
eset = justRMA()
f <- factor(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,

2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2),
labels=c("Healthy", "unaffected"))
design <- model.matrix(~ 0 + f)
design
colnames(design) <-c("Healthy","unaffected")
design
library(limma)
fit <- lmFit(eset, design)
library(hgu133plus2.db)
fit$genes$Symbol <- getSYMBOL(fit$genes$ID,"hgu133plus2.db")
contrast.matrix <-makeContrasts(affected-Healthy,levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
topTable(fit2,coef=1,p=0.05, adjust="fdr")
results <- decideTests(fit2, adjust="fdr", p=0.05)
summary(results)
write.table(results,file="myresults.txt")

The results table shows ~10,000 genes to be upregulated and ~12,000 genes to be down regulated.

My question is how can I get fold change values associated with these genes in the results file?

Thanks in advance, Roopa

Microarray limma • 1.8k views

ADD COMMENT • link updated 11.0 years ago by Gordon Smyth 53k • written 12.8 years ago by Roopa Subbaiaih ▴ 130

score 0 · Answer 1 · 2013-01-31

0

Entering edit mode

Gordon Smyth 53k

@gordon-smyth

Last seen 17 minutes ago

WEHI, Melbourne, Australia

By using write.fit().

Best wishes
Gordon

ADD COMMENT • link 12.8 years ago • updated 11.0 years ago Gordon Smyth 53k

Gordon Smyth · Answer 2 · 2013-01-31

Hi Roopa, > results <- decideTests(fit2, adjust="fdr", p=0.05) > summary(results) > write.table(results,file="myresults.txt") > > The results table shows ~10,000 genes to be upregulated and ~12,000 > genes > to be down regulated. > > My question is how can I get fold change values associated with these > genes > in the results file? You can create object with needed columns binded: e.g. x=cbind(fit$coefficients,fit$p.value,p.adjust(fit$p.value,"BH"),separa te) Where "fit" is the result of lmFit and then eBayes commands - it contains "coefficients" column with (mean) log2 fold change for each gene. *In fact it haven't been clear for me at first - but I compared it to the output of topTable (logFC column) and it is equal.* Back to the table - proposed example will give you: log2 fold change, raw p-value, corrected p-value (here Benjamini-Hochberg, equall to "fdr"), and finally change direction for significant genes (where separate is result of decideTests) Before "write.table" you can change column names with colnames(x)=c(...) HTH -- Maciej Jonczyk, Department of Plant Molecular Ecophysiology Faculty of Biology, University of Warsaw 02-096 Warsaw, Miecznikowa 1 Poland -- This email was Anti Virus checked by Astaro Security Gateway. http://www.astaro.com