Answers.Body: <p>If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, <strong>I'd believe that</strong>, even if looking at the first few PCs doesn't look promising.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-21T16:22:31.293Z
Answers.LastActivityDate: 2020-08-21T16:22:31.293Z
Answers.UserDisplayName: swbarnes2
Body: <p>I have a set of <code>RNA-seq</code> samples from targeting different proteins in a complex with <code>siRNAs</code>. However, the <code>PCA</code> plot shows that <strong>the distance between the two replicates of each treatment is more than the distance between the treatments</strong>. So, there is <strong>no clear border between treatments</strong> to cluster them separately.</p>
  <p>If I go through my usual RNA-seq data analysis pipeline using <code>EdgeR</code> normalization and <code>Limma/Voom</code> DEG analysis, I get a lot of DEGs. However, when I try to visualize them using a heatmap, the replicates do not cluster together, as expected.</p>
  <p>What should I do in this situation? I mean should I trust the DEGs, or try to still extract co-expressed gene modules from this data?</p>
  <p>my code for normalization:</p>
  <pre><code>myDGEList &lt;- DGEList(Txi_gene$counts)
  log2.cpm &lt;- cpm(myDGEList, log=TRUE)
  
  cpm &lt;- cpm(myDGEList)
  
  # prepare your grouping matrix
  group &lt;- factor(targets$group)
  design &lt;- model.matrix(~group)
  
  # filtering
  keepers &lt;- filterByExpr(cpm, design)
  myDGEList.filtered &lt;- myDGEList[keepers,]
  log2.cpm.filtered &lt;- cpm(myDGEList.filtered, log=TRUE)
  
  myDGEList.filtered.norm &lt;- calcNormFactors(myDGEList.filtered, method = &quot;TMM&quot;)
  log2.cpm.filtered.norm &lt;- cpm(myDGEList.filtered.norm, log=TRUE)
  </code></pre>
  <p>my code for PCA plot:</p>
  <pre><code>pca.res &lt;- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
  pc.var &lt;- pca.res<span class="math-container">$sdev^2
  pc.per &lt;- round(pc.var/sum(pc.var)*100, 1) 
  pca.res.df &lt;- as_tibble(pca.res$</span>x)
  
  # PC1 vs PC2
  pca.plot &lt;- ggplot(pca.res.df) +
    aes(x=PC1, y=PC2, label=sampleLabels, color = group, shape = batch) +
    geom_point(size=4) +
    xlab(paste0(&quot;PC1 (&quot;,pc.per[1],&quot;%&quot;,&quot;)&quot;)) + 
    ylab(paste0(&quot;PC2 (&quot;,pc.per[2],&quot;%&quot;,&quot;)&quot;)) +
    labs(title=&quot;PCA plot&quot;,
         caption=paste0(&quot;produced on &quot;, Sys.time())) +
    coord_fixed() +
    theme_bw()
  
  ggplotly(pca.plot)
  </code></pre>
  <p>the PCA plot figure. The legend is cropped out for privacy reasons. Colors represent replicates. Shapes represent batch effect which was has been fixed using <code>combat_seq</code> on raw-data.
  <a href="https://i.sstatic.net/QxbWh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/QxbWh.png" alt="PCA plot" /></a></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-20T20:22:44.993Z
Comments.Text: Can you show the pca plot you're describing?
Comments.UserDisplayName: winni2k
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-21T09:54:40.947Z
Comments.Text: Do you have a plot before Combat batch correction? It may be easier to add the batch information to the edgeR model formula and not correct in in the counts.
Comments.UserDisplayName: PPK
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-21T10:03:17.350Z
Comments.Text: @PPK I do have a plot but my replicates are in different batches, so the PC1 positive side is replicate 1 (batch1) of each treatment and PC1 negative side is replicate 2 (batch2) of each treatment in not batch corrected PCA plot. At first, I used EdgeR model for batch effect and just used combat for visualization but I comapred them and the reults were almost identical, so I used combat-seq on my raw data before any normalization. Do you still think I should add a PCA without batch correction to the question?
Comments.UserDisplayName: Reza Rezaei
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-24T10:29:44.127Z
Comments.Text: Can you show the PCAplot without correction, I think your plot above looks weird. Also when you mentioned you found DEGs, which comparison did you find it for? I see many groups here. For the comparion you are interested in, make a histogram of the p-values and see whether there is an enrichment in low p-values
Comments.UserDisplayName: StupidWolf
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-24T10:31:31.827Z
Comments.Text: The answer to your question is, they should be ok, unless you are underestimating the variance by not including relevant covariates. Still I think you have to interpret the data to see whether it makes sense. It would help a lot if you make some of these diagnostic plots, otherwise it is pretty much guessing.
Comments.UserDisplayName: StupidWolf
Tags: <rna-seq><deseq2><edger><pca><limma>
Title: the variation between treatments is less than the variation between replicates in RNA-seq data
ThreadText: # the variation between treatments is less than the variation between replicates in RNA-seq data
  
  (Asked by: username_0 on 2020-08-20)
  
  I have a set of `RNA-seq` samples from targeting different proteins in a complex with `siRNAs`. However, the `PCA` plot shows that **the distance between the two replicates of each treatment is more than the distance between the treatments**. So, there is **no clear border between treatments** to cluster them separately.
  
  If I go through my usual RNA-seq data analysis pipeline using `EdgeR` normalization and `Limma/Voom` DEG analysis, I get a lot of DEGs. However, when I try to visualize them using a heatmap, the replicates do not cluster together, as expected.
  
  What should I do in this situation? I mean should I trust the DEGs, or try to still extract co-expressed gene modules from this data?
  
  my code for normalization:
  
  ```
  myDGEList <- DGEList(Txi_gene$counts)
  log2.cpm <- cpm(myDGEList, log=TRUE)
  
  cpm <- cpm(myDGEList)
  
  # prepare your grouping matrix
  group <- factor(targets$group)
  design <- model.matrix(~group)
  
  # filtering
  keepers <- filterByExpr(cpm, design)
  myDGEList.filtered <- myDGEList[keepers,]
  log2.cpm.filtered <- cpm(myDGEList.filtered, log=TRUE)
  
  myDGEList.filtered.norm <- calcNormFactors(myDGEList.filtered, method = "TMM")
  log2.cpm.filtered.norm <- cpm(myDGEList.filtered.norm, log=TRUE)
  
  ```
  
  my code for PCA plot:
  
  ```
  pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
  pc.var <- pca.res$sdev^2
  pc.per <- round(pc.var/sum(pc.var)*100, 1) 
  pca.res.df <- as_tibble(pca.res$x)
  
  # PC1 vs PC2
  pca.plot <- ggplot(pca.res.df) +
    aes(x=PC1, y=PC2, label=sampleLabels, color = group, shape = batch) +
    geom_point(size=4) +
    xlab(paste0("PC1 (",pc.per[1],"%",")")) + 
    ylab(paste0("PC2 (",pc.per[2],"%",")")) +
    labs(title="PCA plot",
         caption=paste0("produced on ", Sys.time())) +
    coord_fixed() +
    theme_bw()
  
  ggplotly(pca.plot)
  
  ```
  
  the PCA plot figure. The legend is cropped out for privacy reasons. Colors represent replicates. Shapes represent batch effect which was has been fixed using `combat_seq` on raw-data.
  [![PCA plot](https://i.sstatic.net/QxbWh.png)](https://i.sstatic.net/QxbWh.png)
  
  **Tags:** rna-seq, deseq2, edger, pca, limma
  
  > **Comments:**
  > **username_1**: Can you show the pca plot you're describing?  
  > **username_2**: Do you have a plot before Combat batch correction? It may be easier to add the batch information to the edgeR model formula and not correct in in the counts.  
  > **username_0**: @username_2 I do have a plot but my replicates are in different batches, so the PC1 positive side is replicate 1 (batch1) of each treatment and PC1 negative side is replicate 2 (batch2) of each treatment in not batch corrected PCA plot. At first, I used EdgeR model for batch effect and just used combat for visualization but I comapred them and the reults were almost identical, so I used combat-seq on my raw data before any normalization. Do you still think I should add a PCA without batch correction to the question?  
  > **username_3**: Can you show the PCAplot without correction, I think your plot above looks weird. Also when you mentioned you found DEGs, which comparison did you find it for? I see many groups here. For the comparion you are interested in, make a histogram of the p-values and see whether there is an enrichment in low p-values  
  > **username_3**: The answer to your question is, they should be ok, unless you are underestimating the variance by not including relevant covariates. Still I think you have to interpret the data to see whether it makes sense. It would help a lot if you make some of these diagnostic plots, otherwise it is pretty much guessing.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_4 on 2020-08-21. Score: 1
  
  If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, **I'd believe that**, even if looking at the first few PCs doesn't look promising.
  
  
  ---
  
  