<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>

CC BY-SA 4.0
2020-08-21T16:22:31.293Z
2020-08-21T16:22:31.293Z
swbarnes2
<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>

CC BY-SA 4.0
2020-08-20T20:22:44.993Z
Can you show the pca plot you're describing?
winni2k
CC BY-SA 4.0
2020-08-21T09:54:40.947Z
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.
PPK
CC BY-SA 4.0
2020-08-21T10:03:17.350Z
@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?
Reza Rezaei
CC BY-SA 4.0
2020-08-24T10:29:44.127Z
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
StupidWolf
CC BY-SA 4.0
2020-08-24T10:31:31.827Z
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.
StupidWolf
<rna-seq><deseq2><edger><pca><limma>
the variation between treatments is less than the variation between replicates in RNA-seq data
# 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.

**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**: @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?  
> **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.


---

