diff --git a/DESCRIPTION b/DESCRIPTION
index f88d12584323029412d0741698c6b2d0d469ea8d..69b1543f7dcf1fb35813db854047da418c83a950 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -41,7 +41,8 @@ Imports:
     stats,
     SummarizedExperiment,
     utils,
-    crayon
+    crayon,
+    S4Vectors
 Suggests:
     BiocStyle,
     covr,
@@ -53,7 +54,6 @@ Suggests:
     pscl,
     testthat,
     rmarkdown,
-    S4Vectors,
     scDD,
     scran,
     mfa,
diff --git a/NAMESPACE b/NAMESPACE
index d9df705e0105c334ecd5f24e209a0e75f8a2086e..c9b65bdf8f29dc8cb1fdd86c22292094c7709fd9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -78,8 +78,11 @@ exportClasses(SplatParams)
 exportClasses(ZINBParams)
 importFrom(BiocParallel,SerialParam)
 importFrom(BiocParallel,bplapply)
+importFrom(S4Vectors,"metadata<-")
+importFrom(S4Vectors,metadata)
 importFrom(SingleCellExperiment,"cpm<-")
 importFrom(SingleCellExperiment,SingleCellExperiment)
+importFrom(SingleCellExperiment,cpm)
 importFrom(SummarizedExperiment,"assays<-")
 importFrom(SummarizedExperiment,"colData<-")
 importFrom(SummarizedExperiment,"rowData<-")
@@ -96,15 +99,19 @@ importFrom(checkmate,checkIntegerish)
 importFrom(checkmate,checkNumber)
 importFrom(checkmate,checkNumeric)
 importFrom(ggplot2,aes_string)
+importFrom(ggplot2,coord_fixed)
 importFrom(ggplot2,element_blank)
+importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_abline)
 importFrom(ggplot2,geom_boxplot)
 importFrom(ggplot2,geom_hline)
 importFrom(ggplot2,geom_point)
 importFrom(ggplot2,geom_smooth)
+importFrom(ggplot2,geom_tile)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggtitle)
 importFrom(ggplot2,scale_colour_manual)
+importFrom(ggplot2,scale_fill_distiller)
 importFrom(ggplot2,scale_fill_manual)
 importFrom(ggplot2,scale_x_log10)
 importFrom(ggplot2,scale_y_continuous)
@@ -122,6 +129,7 @@ importFrom(methods,slot)
 importFrom(methods,slotNames)
 importFrom(methods,validObject)
 importFrom(stats,aggregate)
+importFrom(stats,cor)
 importFrom(stats,dnbinom)
 importFrom(stats,ks.test)
 importFrom(stats,median)
diff --git a/R/compare.R b/R/compare.R
index a00984682790070d5ddf39b63a0d4b55d5c4815a..dd556d85548aa7eeff2e5c5063d422d2bf70bf46 100644
--- a/R/compare.R
+++ b/R/compare.R
@@ -50,9 +50,12 @@
 #' names(comparison)
 #' names(comparison$Plots)
 #' @importFrom ggplot2 ggplot aes_string geom_point geom_smooth geom_boxplot
-#' scale_y_continuous scale_y_log10 scale_x_log10 xlab ylab ggtitle
-#' theme_minimal scale_colour_manual scale_fill_manual
-#' @importFrom SingleCellExperiment cpm<-
+#' geom_tile scale_y_continuous scale_y_log10 scale_x_log10 scale_colour_manual
+#' scale_fill_manual scale_fill_distiller coord_fixed facet_wrap xlab ylab
+#' ggtitle theme_minimal
+#' @importFrom S4Vectors metadata<- metadata
+#' @importFrom SingleCellExperiment cpm<- cpm
+#' @importFrom stats cor
 #' @export
 compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
                         fits = TRUE, colours = NULL) {
@@ -81,22 +84,37 @@ compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
         sce <- addFeatureStats(sce, "cpm", log = TRUE)
         n.features <- colData(sce)$total_features_by_counts
         colData(sce)$PctZero <- 100 * (1 - n.features / nrow(sce))
+
+        var.genes <- rev(order(rowData(sce)$VarLogCPM))[1:100]
+        var.cpm <- log2(cpm(sce)[var.genes, ] + 1)
+        var.cors <- as.data.frame.table(cor(t(var.cpm), method = "spearman"))
+        colnames(var.cors) <- c("GeneA", "GeneB", "Correlation")
+        var.cors$VarGeneA <- rep(paste0("VarGene", 1:100), 100)
+        var.cors$VarGeneB <- rep(paste0("VarGene", 1:100), each = 100)
+        var.cors$Dataset <- name
+        var.cors <- var.cors[, c("Dataset", "GeneA", "GeneB", "VarGeneA",
+                                 "VarGeneB", "Correlation")]
+        metadata(sce)$VarGeneCorrelation <- var.cors
+
         sces[[name]] <- sce
     }
 
     features <- rowData(sces[[1]])
     cells <- colData(sces[[1]])
+    var.cors <- metadata(sces[[1]])$VarGeneCorrelation
 
     if (length(sces) > 1) {
         for (name in names(sces)[-1]) {
             sce <- sces[[name]]
             features <- rbindMatched(features, rowData(sce))
             cells <- rbindMatched(cells, colData(sce))
+            var.cors <- rbindMatched(var.cors, metadata(sce)$VarGeneCorrelation)
         }
     }
 
     features$Dataset <- factor(features$Dataset, levels = names(sces))
     cells$Dataset <- factor(cells$Dataset, levels = names(sces))
+    var.cors$Dataset <- factor(var.cors$Dataset, levels = names(sces))
     features <- data.frame(features)
     cells <- data.frame(cells)
 
@@ -172,6 +190,18 @@ compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
         ggtitle("Mean-zeros relationship") +
         theme_minimal()
 
+    var.correlation <- ggplot(var.cors,
+                              aes_string(x = "VarGeneA", y = "VarGeneB",
+                                         fill = "Correlation")) +
+        geom_tile() +
+        scale_fill_distiller(palette = "RdBu", limits = c(-1, 1)) +
+        coord_fixed() +
+        facet_wrap(~ Dataset) +
+        ggtitle("Correlation - 100 variable genes") +
+        theme_minimal() +
+        theme(axis.title = element_blank(),
+              axis.text = element_blank())
+
     if (fits) {
         mean.var <- mean.var + geom_smooth(method = "gam",
                                            formula = y ~ s(x, bs = "cs"))
@@ -187,7 +217,8 @@ compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
                                     LibrarySizes = libs,
                                     ZerosGene = z.gene,
                                     ZerosCell = z.cell,
-                                    MeanZeros = mean.zeros))
+                                    MeanZeros = mean.zeros,
+                                    VarGeneCor = var.correlation))
 
     return(comparison)
 }
diff --git a/tests/testthat/test-compare.R b/tests/testthat/test-compare.R
index 29540dfc3be64a9620582d1f3d7dfffbc3da2026..960cbf44b4c8abf9ab607297b505e8d25f302484 100644
--- a/tests/testthat/test-compare.R
+++ b/tests/testthat/test-compare.R
@@ -12,9 +12,9 @@ test_that("compareSCEs works", {
                         names(comparison)))
     checkmate::expect_class(comparison$ColData, "data.frame")
     checkmate::expect_class(comparison$RowData, "data.frame")
-    expect_length(comparison$Plots, 7)
+    expect_length(comparison$Plots, 8)
     expect_true(all(c("Means", "Variances", "MeanVar", "LibrarySizes",
-                      "ZerosGene", "ZerosCell", "MeanZeros") %in%
+                      "ZerosGene", "ZerosCell", "MeanZeros", "VarGeneCor") %in%
                         names(comparison$Plots)))
     for (plot in names(comparison$Plots)) {
         checkmate::expect_class(comparison$Plots[[plot]], "ggplot")