I am using edgeR to perform differential expression (DE) analysis on a set of RNA-seq data samples (2 controls; 8 treatments). To correct for batch effects, I am using RUVSeq.
I am able to get a list of DE genes without normalization:
x <- as.factor(rep(c("Ctl","Inf"),c(2,8)))
set <- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, "DE_genelist.txt", sep=""))
Then immediately after creating the "top" object, I use RUVg to normalize:
# [...]
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg <- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, "DE_RUVg_genelist.txt", sep=""))
And I get the error:
Error in as.data.frame.default(x[[i]], optional = TRUE) :
cannot coerce class ‘structure("SeqExpressionSet", package = "EDASeq")’ to a data.frame
I am not sure how to print the list of normalized results like I can with the unnormalized data. Ideally, I would get a file with the same format as the edgeR output (as a .csv or .txt file):
"logFC" "logCPM" "LR" "PValue" "FDR"
"COBLL1" -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
"UBE2D1" -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
"NEK7" -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
"SMC6" -2.300 5.674738981329 61.8130019860261 3.7767230643666e-15 3.4974443325016e-12
How can I get a list of genes as an output after normalization with RUVSeq?