How to create a list of differentially expressed (DE) genes after normalization with RUVSeq? 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 &lt;- as.factor(rep(c(&quot;Ctl&quot;,&quot;Inf&quot;),c(2,8)))
set &lt;- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design &lt;- model.matrix(~x, data=pData(set))
y &lt;- DGEList(counts=counts(set), group=x)
y &lt;- calcNormFactors(y, method=&quot;upperquartile&quot;)
y &lt;- estimateGLMCommonDisp(y, design)
y &lt;- estimateGLMTagwiseDisp(y, design)
fit &lt;- glmFit(y, design)
lrt &lt;- glmLRT(fit, coef=2)
top &lt;- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, &quot;DE_genelist.txt&quot;, sep=&quot;&quot;))

Then immediately after creating the &quot;top&quot; object, I use RUVg to normalize:
# [...]
top &lt;- topTags(lrt, n=nrow(set))$table
empirical &lt;- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg &lt;- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, &quot;DE_RUVg_genelist.txt&quot;, sep=&quot;&quot;))

And I get the error:
Error in as.data.frame.default(x[[i]], optional = TRUE) : 
  cannot coerce class ‘structure(&quot;SeqExpressionSet&quot;, package = &quot;EDASeq&quot;)’ 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):
&quot;logFC&quot; &quot;logCPM&quot; &quot;LR&quot; &quot;PValue&quot; &quot;FDR&quot;
&quot;COBLL1&quot; -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
&quot;UBE2D1&quot; -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
&quot;NEK7&quot; -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
&quot;SMC6&quot; -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?
You do the normalization before running your edgeR. The purpose of RUVg is to remove &quot;Remove Unwanted Variation Using Control Genes&quot;. In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.
Using the example dataset in vignette:
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter &lt;- apply(zfGenes, 1, function(x) length(x[x&gt;5])&gt;=2)
filtered &lt;- zfGenes[filter,]
genes &lt;- rownames(filtered)[grep(&quot;^ENS&quot;, rownames(filtered))]
spikes &lt;- rownames(filtered)[grep(&quot;^ERCC&quot;, rownames(filtered))]

x &lt;- as.factor(rep(c(&quot;Ctl&quot;, &quot;Trt&quot;), each=3))
set &lt;- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set &lt;- betweenLaneNormalization(set, which=&quot;upper&quot;)

set1 &lt;- RUVg(set, spikes, k=1)

You can look at it, it's an expression set with counts etc, not results:
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples 
  element names: counts, normalizedCounts, offset 
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  

You run edgeR now on the results of RUVg:
design &lt;- model.matrix(~x + W_1, data=pData(set1))
y &lt;- DGEList(counts=counts(set1), group=x)
y &lt;- calcNormFactors(y, method=&quot;upperquartile&quot;)
y &lt;- estimateGLMCommonDisp(y, design)
y &lt;- estimateGLMTagwiseDisp(y, design)
fit &lt;- glmFit(y, design)
lrt &lt;- glmLRT(fit, coef=2)
topTags(lrt)