Hello everyone,
I have RNA-seq count data that includes 8 different sample groups with unequal replication. I want to visualize this count data for all 8 sample groups in a single heatmap. The heatmap should clearly differentiate between up-regulated and down-regulated genes among the different sample groups. However, the differentially expressed genes identified through pairwise comparisons using DESeq2 analysis are not showing clear differentiation in the heatmap plot. I have tried using log2 TPM normalized values for the heatmap. Is my analysis statistically correct, or should I approach this differently? Please help me.
1 answer
You have to scale the counts first to emphasize differences. Clustering based on unscaled values emphasizes the magnitude of counts, not differences. See for a minimal example Scaling RNA-Seq data before clustering?
Yeah... you're right. Now I have done according to the post. The code I have used is
heat_df_centered <- t(apply(heat_df, 1, function(x) {
if(all(is.na(x))) {
return(x)
} else {
return((x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE))
}
}))
# Handle non-finite values
heat_df_centered[!is.finite(heat_df_centered)] <- 0
Heatmap(heat_df_centered,
cluster_columns = FALSE,
col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
clustering_method_rows = "ward.D2",
name = "log2TPM\nmean centered",
top_annotation = ha)
The code for scaling is t(scale(t(x))), don't use apply here, that is inefficent. Before command is vectorized. If you use DEGs then by definition there cannot be rows that are all zero. If there would be, you could still do scaled[complete.cases(scaled),]. But what you do is basically correct in a technical sense. Now it's down to interpret the results.
Log in to answer this question.
-
log2Fold Change as input for k-means analysis
written by concetta •
-
filter out gene with zero counts in differential gene analysis
written by tommy •
-
Correlation analyses on TCGA data
written by mario.keller.1988 •
-
Using RPKM to compare gene expression between samples
written by kz •
-
Reactome website gene expression analysis
written by camillab. •
-
DESeq2 for pairwise comparison of multiple groups
written by thjnant
-
Unequal sample sizes Deseq2
written by Payal
-
Differentially expressed genes visualization with Cytoscape
written by Sib •
-
Quantile normalizing prior to or after TPM scaling?
written by mforde84
-
RNA-seq replicates pooled before sequencing. How to proceed with DGE analysis?
written by Sentinel156
Show example images and the code used to generate them. We can't see what you see, so it's hard to say if you are plotting things incorrectly or if it's inherent to the variability in your dataset.
Generally, fold-changes can be hard to distinguish from TPMs between samples when TPM variability between genes is large. Gene-wise scaling (e.g. Z score where mean expression across samples is 0) can help highlight the differences.
Also, if you have non-differentially expressed genes included in the heatmap, that can make it harder to see differences in TPMs between samples. You can try a heatmap with only DEGs if that's the case and you insist on showing expression values.
As mentioned though, code and an example image would be helpful in answering your question.
Thank you @rfran010 and jared.andrews07 for your comments. I have used the DEGs only for my heatmap. The TPM normalized values were calculated from the original count data by using the 'bioinfokit=2.1.4' python package.(https://github.com/reneshbedre/bioinfokit). And I have used the below code for creating the heatmap.
Among many of the clustering methods I have found that the 'complete' method was separating the clusters very well. And I got the heatmap as below. Please give your comments about the correctness of this plot and it's interpretation.
