0
$\begingroup$

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

Improve this question
$\endgroup$
6
  • $\begingroup$ Can you show the pca plot you're describing? $\endgroup$
    – winni2k
    Aug 20, 2020 at 20:22
  • $\begingroup$ Please include plots and all relevant code. It is difficult to argue only on words in this context. $\endgroup$
    – user3051
    Aug 21, 2020 at 7:25
  • $\begingroup$ 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. $\endgroup$
    – PPK
    Aug 21, 2020 at 9:54
  • $\begingroup$ @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? $\endgroup$
    – Reza Rezaei
    Aug 21, 2020 at 10:03
  • $\begingroup$ 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 $\endgroup$
    – StupidWolf
    Aug 24, 2020 at 10:29

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.