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.
Here is the desing, targets and the script:
> 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
> 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
>co=as.matrix(read.table("2014_04_02_6h_PB.csv",header=T, sep=",",
row.names=1))
>keep <- rowSums(co)>=9
>nkeep <- sum(keep)
>co2 <- co[keep,]
>nf = calcNormFactors (co2)
>targets= read.table ("targets.csv", header = T, sep=",",row.names=1)
>targets
>ind <- factor (targets$Individual)
>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=ind)
>fit <- lmFit(y,design,block=ind,correlation=corfit$consensus)
>cm <- makeContrasts (co2Vsco=c2-co, levels=design)
>fit2 <- eBayes(contrasts.fit(fit,cm))
>res <- topTable (fit2,coef="co2Vsco", n=nrow (y), sort.by=ânoneâ)
>res <- data.frame (res)
>res2 <- subset (res, adj.P.Val < 0.1)
> 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]]
______________________________________________
[email protected] 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.
