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.steve@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 <catagui@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]]


