Hi Steve, Thanks for you reply. I have tried what you suggested and checked the limma guide sessions. Now I am getting the results from only one treatment which is great but the problem is all my adjusted Pvalues (adj.P.Val) are >0.05.
Don't really know how to solve this, I know my data works because I used DESeq and got around 200 differentially expressed genes FDR < 0.05. >From: Steve Lianoglou <lianoglou.st...@gene.com> >Date: Mon, Apr 7, 2014 at 5:09 PM >Subject: Re: [R] How to make a proper use of blocking in limma using >voom >To: Catalina Aguilar Hurtado <cata...@gmail.com> >Cc: r help <r-help@r-project.org> >Hi, >This is because you defined your model matrix to have no intercept >term (that's what the 0 in `~ 0 + treat` is doing). >You will have to test for a particular contrast (not just a coeficient >in your design) in order to get what you are after. Sections 9.7 ( >Multi-level Experiments) and 16.3 ( >Comparing Mammary Progenitor Cell Populations with Illumina BeadChips) >in the lima user's guide may be the most useful for you to follow >along at this point: >http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf<http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf> >(look for the `makeContrasts` call) >After you call `fit <- lmFit( ... )` in your code above you might do something like: >R> cm <- makeContrasts(co2Vsco=treatco2 - treatco0, levels=design) >R> fit2 <- eBayes(contrasts.fit(fit, cm)) >R> res <- topTable(fit2, coef='co2Vsco') >Note that the `treatco2` and `treatco` are only correct if these are >the column names of your design matrix -- substitute with the >appropriate names for your example, if necessary. >> Thanks in advance for the help. >Noodle on that a bit and if you still have questions, please do post a >follow up question on the bioconductor list. >Btw, to help make your question more interpretable, since we don't >have your targets file, I think it would be easier for us if you >copy/paste the output of `dput(targets)` and `dput(design)` after you >create those object in a follow up email if it's necessary to write >one. >HTH, >-steve Here is the desing, targets and the script: R> design co c2 1 1 0 2 1 0 3 1 0 4 1 0 5 0 1 6 0 1 7 0 1 8 0 1 9 0 1 R> targets Individual condition B1 1 co B3 3 co B4 4 co B5 5 co D1 1 c2 D2 2 c2 D3 3 c2 D4 4 c2 D5 5 c2 R>co=as.matrix(read.table("2014_04_02_6h_PB.csv",header=T, sep=",", row.names=1)) R>keep <- rowSums(co)>=9 R>nkeep <- sum(keep) R>co2 <- co[keep,] R>nf = calcNormFactors (co2) R>targets= read.table ("targets.csv", header = T, sep=",",row.names=1) R>targets R>ind <- factor (targets$Individual) R>treat <- factor (targets$condition, levels= c("co", "c2")) R>design <- model.matrix(~0 +treat) R>colnames (design) <- levels (treat) R>y <- voom(co2,design,lib.size=colSums(co2)*nf) R>corfit <- duplicateCorrelation(y,design,block=ind) R>fit <- lmFit(y,design,block=ind,correlation=corfit$consensus) R>cm <- makeContrasts (co2Vsco=c2-co, levels=design) R>fit2 <- eBayes(contrasts.fit(fit,cm)) R>res <- topTable (fit2,coef="co2Vsco", n=nrow (y), sort.by=ânoneâ) R>res <- data.frame (res) R>res2 <- subset (res, adj.P.Val < 0.1) R> res2 [1] ID logFC AveExpr t P.Value adj.P.Val B <0 rows> (or 0-length row.names) Thanks again, Catalina [[alternative HTML version deleted]]
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
