Hi, This is a bioconductor-related question so please post any follow up questions on that mailing list. You can sign up to that list here:
https://stat.ethz.ch/mailman/listinfo/bioconductor Comments in line: On Sun, Apr 6, 2014 at 10:41 PM, Catalina Aguilar Hurtado <cata...@gmail.com> wrote: > Hi all, > > I have a RNAseq data to analyse were I have a control and a one treatment > for different individuals. I need to block the effects of the individual, > but I am having several troubles to get the data that I need. I am using > voom because my data is very heterogeneous and voom seams to do a good job > normalising my reads. > > I am having the following issues: > > 1. > > I want to get the differentially expressed genes (DEGs) of my treatment > not of my control. I don't understand after the eBayes analysis why I get > the coefficients for both. I have tried a > makeContrasts (TreatvsCont= > c2-co, levels = design) to subtract the control effect but then I get 0 > DEGs. > 2. > > I am not sure when to include the 0 (null model) in the model formula, I > have read examples for both types of models. > > This are my targets, with my column names of my counts, individual and > condition > >>targets > > Individual condition > > A1 1 co > > A2 3 co > > A4 4 co > > A5 5 co > > E1 1 c2 > > E2 2 c2 > > E3 3 c2 > > E4 4 c2 > > E5 5 c2 > > This is the code I have been trying: > >>co2=as.matrix(read.table("2014_04_02_1h_PB.csv",header=T, sep=",", > row.names=1)) > >>nf = calcNormFactors (co2) > >>targets= read.table ("targets.csv", header = T, sep=",",row.names=1) > >>treat <- factor (targets$condition, levels= c("co", "c2")) > >>design <- model.matrix(~0+treat) > >>colnames (design) <- levels (treat) > >>y <- voom(co2,design,lib.size=colSums(co2)*nf) > >>corfit <- duplicateCorrelation(y,design,block=targets$Individual) > >>fit <- > lmFit(y,design,block=targets$Individual,correlation=corfit$consensus) > >>fit2<- eBayes (fit) > >>results_trt <- topTable (fit2, coef="c2", n=nrow (y), sort.by="none") > > >From which gives me 18,000 genes with adj.P.Val < 0.01 out of 22,000 genes > that I have in total. Which makes no sense.. 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 (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 -- Steve Lianoglou Computational Biologist Genentech ______________________________________________ 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.
