Query regarding limma
2
0
Entering edit mode
@roopa-subbaiaih-5490
Last seen 9.7 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.5k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia

By using write.fit().

Best wishes
Gordon

ADD COMMENT
0
Entering edit mode
boczniak767 ▴ 740
@maciej-jonczyk-3945
Last seen 5 weeks ago
Poland
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
ADD COMMENT
0
Entering edit mode

Thank you both. It worked.

Roopa

ADD REPLY

Login before adding your answer.

Traffic: 503 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