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,"%",")")) + ylab(paste0("PC2 (",pc.per,"%",")")) + labs(title="PCA plot", caption=paste0("produced on ", Sys.time())) + coord_fixed() + theme_bw() ggplotly(pca.plot)