text: user9299: 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.
  
  
  user2388: 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.
  
original_text: 9299: <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.stack.imgur.com/QxbWh.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QxbWh.png" alt="PCA plot" /></a></p>
   2388: <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>
  