Answers.Body: <p>You do the normalization before running your edgeR. The purpose of <code>RUVg</code> 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.</p>
  <p>Using the example dataset in <a href="https://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.R" rel="nofollow noreferrer">vignette</a>:</p>
  <pre><code>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)
  </code></pre>
  <p>You can look at it, it's an expression set with counts etc, not results:</p>
  <pre><code>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:  
  </code></pre>
  <p>You run edgeR now on the results of RUVg:</p>
  <pre><code>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)
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-01T23:35:55.507Z
Answers.Text: This is great! Definitely answers the question. [Although the very last line `topTags(lrt)` will only provide a list of 10 genes, so I'll use `topTags(lrt, n=nrow(set))$table` for a full list].
Answers.UserDisplayName: Gawain
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-01T23:36:39.107Z
Answers.Text: Follow-up question for you then: What is the utility of performing the "upper-quartile" (UQ) normalization before using RUVg and then running edgeR? My thought was, "if I am using RUVg to correct for batch effects, why should I also use UQ?" Perhaps it is important?
Answers.UserDisplayName: Gawain
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-02T14:15:49.983Z
Answers.Text: This is what I know about UQ normalization: scaling with or without UQ before RUVg will usually lead to similar results. It is only important to include the edgeR offset -- `calcNormFactors(y, method="upperquartile")` -- if you decide to scale prior to RUV use. If no UQ scaling is applied, the edgeR offset should be set to 0 (zero).
Answers.UserDisplayName: Gawain
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-02T21:57:49.073Z
Answers.Text: The betweenLaneNormalization() in limma adjust for signal intensity difference (in arrays) or in edgeR, you try to normalize between samples such that most genes are unchanged. RUVg adjust for unseen variation on top of this, for example batch effects. So you don't want to introduce this as a kind of variation for RUVg model to pick up. Hence you do the normalization before. And as you noted, if you use edgeR,  RUV also ask you to provide the offset , which is the normalization factors.
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-02T21:58:24.970Z
Answers.Text: Hope this answers your question. You can try and see for yourself what is the offset, and what RUVg does on top of it
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-03T13:34:11.243Z
Answers.Text: Wonderful explanation! Thanks for your help @StupidWolf!!
Answers.UserDisplayName: Gawain
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-01T22:36:05.630Z
Answers.LastActivityDate: 2020-07-01T22:36:05.630Z
Answers.UserDisplayName: StupidWolf
Answers.Body: <p>I have not used this package but from your code it seems that <code>ruvg</code> is not a table. Instead, it is an R object, which means that you cannot use <code>write.table</code>. I think the results you want is stored in the object. All R objects contain &quot;slots&quot; of data, which can be accessed by <code>@</code>. If I were you, I would type <code>ruvg@</code> and should be able to see which data slots are contained in the object.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-01T21:57:53.413Z
Answers.Text: This certainly gives me more direction than I had before, so I appreciate that. Just looking at all the data stored by the ruvg object (Large SeqExpressionSet (1.1 Mb)) tells me that it may take some time to find what I am looking for in there. I'll report back if I find the results I want. Thanks for your suggestion!
Answers.UserDisplayName: Gawain
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-01T14:45:43.413Z
Answers.LastActivityDate: 2020-07-01T14:45:43.413Z
Answers.UserDisplayName: Phoenix Mu
Body: <p>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.</p>
  <p>I am able to get a list of DE genes without normalization:</p>
  <pre><code>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;))
  </code></pre>
  <p>Then immediately after creating the &quot;top&quot; object, I use RUVg to normalize:</p>
  <pre><code># [...]
  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;))
  </code></pre>
  <p>And I get the error:</p>
  <pre><code>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
  </code></pre>
  <p>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):</p>
  <pre><code>&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
  </code></pre>
  <p><strong>How can I get a list of genes as an output <em>after</em> normalization with RUVSeq?</strong></p>
  
Tags: <r><rna-seq><differential-expression><normalization><edger>
Title: How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?
ThreadText: # How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?
  
  (Asked by: username_0 on 2020-06-30)
  
  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?**
  
  **Tags:** r, rna-seq, differential-expression, normalization, edger
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2020-07-01. Score: 1 | Accepted answer
  
  You do the normalization before running your edgeR. The purpose of `RUVg` is to remove "Remove Unwanted Variation Using Control Genes". 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](https://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.R):
  
  ```
  library(RUVSeq)
  library(zebrafishRNASeq)
  data(zfGenes)
  filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
  filtered <- zfGenes[filter,]
  genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
  spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]
  
  x <- as.factor(rep(c("Ctl", "Trt"), each=3))
  set <- newSeqExpressionSet(as.matrix(filtered),
                             phenoData = data.frame(x, row.names=colnames(filtered)))
  set <- betweenLaneNormalization(set, which="upper")
  
  set1 <- 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 <- model.matrix(~x + W_1, data=pData(set1))
  y <- DGEList(counts=counts(set1), group=x)
  y <- calcNormFactors(y, method="upperquartile")
  y <- estimateGLMCommonDisp(y, design)
  y <- estimateGLMTagwiseDisp(y, design)
  fit <- glmFit(y, design)
  lrt <- glmLRT(fit, coef=2)
  topTags(lrt)
  
  ```
  
  > **Comments:**
  > **username_0**: This is great! Definitely answers the question. [Although the very last line `topTags(lrt)` will only provide a list of 10 genes, so I'll use `topTags(lrt, n=nrow(set))$table` for a full list].  
  > **username_0**: Follow-up question for you then: What is the utility of performing the "upper-quartile" (UQ) normalization before using RUVg and then running edgeR? My thought was, "if I am using RUVg to correct for batch effects, why should I also use UQ?" Perhaps it is important?  
  > **username_0**: This is what I know about UQ normalization: scaling with or without UQ before RUVg will usually lead to similar results. It is only important to include the edgeR offset -- `calcNormFactors(y, method="upperquartile")` -- if you decide to scale prior to RUV use. If no UQ scaling is applied, the edgeR offset should be set to 0 (zero).  
  > **username_1**: The betweenLaneNormalization() in limma adjust for signal intensity difference (in arrays) or in edgeR, you try to normalize between samples such that most genes are unchanged. RUVg adjust for unseen variation on top of this, for example batch effects. So you don't want to introduce this as a kind of variation for RUVg model to pick up. Hence you do the normalization before. And as you noted, if you use edgeR,  RUV also ask you to provide the offset , which is the normalization factors.  
  > **username_1**: Hope this answers your question. You can try and see for yourself what is the offset, and what RUVg does on top of it  
  > **username_0**: Wonderful explanation! Thanks for your help @username_1!!  
  
  ---
  
  ### Answer by username_2 on 2020-07-01. Score: 0
  
  I have not used this package but from your code it seems that `ruvg` is not a table. Instead, it is an R object, which means that you cannot use `write.table`. I think the results you want is stored in the object. All R objects contain "slots" of data, which can be accessed by `@`. If I were you, I would type `ruvg@` and should be able to see which data slots are contained in the object.
  
  > **Comments:**
  > **username_0**: This certainly gives me more direction than I had before, so I appreciate that. Just looking at all the data stored by the ruvg object (Large SeqExpressionSet (1.1 Mb)) tells me that it may take some time to find what I am looking for in there. I'll report back if I find the results I want. Thanks for your suggestion!  
  
  ---
  
  