Comparison of gene counts and DESeq results: Magenetique vs. nf-core/rnaseq

library(ggplot2)
library(ggcorrplot)
library(BioVenn)
library(dplyr)

# Set correct dir
setwd("/s/public_webshare/public/ghga/magnetique-harmonized-rnaseq/nf-core_vs_magnetique//")

Correlation of counts

Example for one sample match

# Load original data
load("../orignal_data/data/MAGNetApp/data/MAGNet.RData")
data_mag <- countData
# Subset for C00039
countC00039_mag <- data_mag[, c("C00039"), drop=F]
colnames(countC00039_mag) <- 'magnetique'

# Repeat for nf-core/rnaseq 
load("../DESeq2/data_prep/MAGNet.RData")
data_nf <- countData
countC00039_nf <- data_nf[, c("C00039"), drop=F]
colnames(countC00039_nf) <- 'nf-core'

# Merge and plot
merge <- merge(countC00039_mag, countC00039_nf, by = "row.names")
ggplot(merge, aes(x = magnetique, y = `nf-core`)) + geom_point(color="#2171b5") + theme_bw() +
  geom_abline(slope=1, intercept=0, linetype="dotted") +
  scale_x_log10() + scale_y_log10()

# correlation
cor(merge$magnetique, merge$`nf-core`)

[1] 0.7970865

All vs all

# save common gene/sample list to subset full data
common_genes <- merge$Row.names
# only choose a subset of common samples
common_samples_30 <- intersect(colnames(data_mag), colnames(data_nf))[1:30]
common_samples_end30 <- tail(intersect(colnames(data_mag), colnames(data_nf)), 30)
data_mag_sub <- data_mag[common_genes, common_samples_30]
data_nf_sub <- data_nf[common_genes, rev(colnames(data_mag_sub))]

# correlation matrix
cor_mat <- cor(data_mag_sub, data_nf_sub)

# plot first 30 samples
ggcorrplot(cor_mat) + scale_fill_gradient2(limit = c(0.5,1), low = "blue", high =  "red", mid = "white", midpoint = 0.75)

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Using normalised count data

vsd_nf <- readRDS("../DESeq2/harmonized-deseq/harmonized_DESeqTransform_vsd.rds")
vsd_nf_mat <- assay(vsd_nf)
normC00039_nf <- vsd_nf_mat[, c("C00039"), drop=F]
colnames(normC00039_nf) <- 'nf-core'

dds_mag <- readRDS("../DESeq2/magnetiqueCode2022/MAGNet_DESeqDataSet.rds")
vsd_mag <- vst(dds_mag, blind=FALSE)
vsd_mag_mat <- assay(vsd_mag)
normC00039_mag <- vsd_mag_mat[, c("C00039"), drop=F]
colnames(normC00039_mag) <- 'magnetique'

# Merge and plot
merge_norm <- merge(normC00039_mag, normC00039_nf, by = "row.names")
ggplot(merge_norm, aes(x = magnetique, y = `nf-core`)) + geom_point(color="#2171b5") + theme_bw() +
  geom_abline(slope=1, intercept=0, linetype="dotted")

cor(merge_norm$magnetique, merge_norm$`nf-core`)

[1] 0.9833694

All vs all - norm

# save common gene/sample list to subset full data
common_genes <- merge_norm$Row.names
# only choose a subset of common samples
common_samples_30 <- intersect(colnames(vsd_mag_mat), colnames(vsd_nf_mat))[1:30]
common_samples_end30 <- tail(intersect(colnames(vsd_mag_mat), colnames(vsd_nf_mat)), 30)
vsd_mag_sub <- vsd_mag_mat[common_genes, common_samples_30]
vsd_nf_sub <- vsd_nf_mat[common_genes, rev(colnames(vsd_mag_sub))]

# correlation matrix
cor_mat_norm <- cor(vsd_mag_sub, vsd_nf_sub)

# plot first 30 samples
ggcorrplot(cor_mat_norm) + scale_fill_gradient2(limit = c(0.8,1), low = "blue", high =  "red", mid = "white", midpoint = 0.9)

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

PCA plots

# nf-core
pca_plt_nf <- plotPCA(vsd_nf, intgroup=c("Etiology"), returnData = T)
percentVar <- round(100 * attr(pca_plt_nf, "percentVar"))
ggplot(pca_plt_nf, aes(PC1, PC2, color=Etiology)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() +
  theme_bw()

# Mag
pca_plt_mag <- plotPCA(vsd_mag, intgroup=c("Etiology"), returnData = T)
percentVar <- round(100 * attr(pca_plt_mag, "percentVar"))
ggplot(pca_plt_mag, aes(PC1, PC2, color=Etiology)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() +
  theme_bw()

Comaprison of DESeq2 results

DCM vs HCM

# load results tables
deseq_DCMvsHCM_nf <- readRDS("../DESeq2/harmonized-deseq/harmonized_DESeqResults.rds")$DCMvsHCM
deseq_DCMvsHCM_mag <- readRDS("../orignal_data/data/MAGNetApp/data/DGE/MAGNet_DCMvsHCM_GeneTonic.rds")$res_de

# filter
DEgenes_DCMvsHCM_mag <- rownames(subset(deseq_DCMvsHCM_mag, padj <= 0.05))
DEgenes_DCMvsHCM_nf <- rownames(subset(deseq_DCMvsHCM_nf, padj <= 0.05))

# compare gene lists
overlap_DCMvsHCM <- data.frame(genes = intersect(DEgenes_DCMvsHCM_mag, DEgenes_DCMvsHCM_nf))
overlap_DCMvsHCM$type <- "overlapping"
nf_DCMvsHCM <- data.frame(genes = setdiff(DEgenes_DCMvsHCM_nf,rownames(overlap_DCMvsHCM)))
nf_DCMvsHCM$type <- "nf-core"
mag_DCMvsHCM <- data.frame(genes = setdiff(DEgenes_DCMvsHCM_mag,rownames(overlap_DCMvsHCM)))
mag_DCMvsHCM$type <- "Magnetique"
DCMvsHCM <- rbind(overlap_DCMvsHCM, nf_DCMvsHCM, mag_DCMvsHCM)
paste0("Proportion overlapping DE genes DCMvsHCM: ", nrow(overlap_DCMvsHCM)/nrow(DEgenes_DCMvsHCM_mag))

[1] "Proportion overlapping DE genes DCMvsHCM: "

# venn
biovenn <- draw.venn(DEgenes_DCMvsHCM_mag, DEgenes_DCMvsHCM_nf, NULL, title="DCM vs HCM", nrtype="abs", subtitle = NULL, xtitle = "Magnetique", ytitle = "nf-core")

[1] "x total: 201"
[1] "y total: 149"
[1] "z total: 0"
[1] "x only: 80"
[1] "y only: 28"
[1] "z only: 0"
[1] "x-y total overlap: 121"
[1] "x-z total overlap: 0"
[1] "y-z total overlap: 0"
[1] "x-y only overlap: 121"
[1] "x-z only overlap: 0"
[1] "y-z only overlap: 0"
[1] "x-y-z overlap: 0"

# 
ggplot(deseq_DCMvsHCM_merge_df, aes(x=log2FoldChange, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
    guides(color = "none") +
  xlim(c(-3,3)) +
  facet_grid(~method)

ggplot(deseq_DCMvsHCM_merge_df, aes(x=baseMean, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
    guides(color = "none") +
  xlim(c(10, 1000)) +
  facet_grid(~method, scales = "free_y")

DCM vs NFD

# load results tables
deseq_DCMvsNFD_nf <- readRDS("../DESeq2/harmonized-deseq/harmonized_DESeqResults.rds")$DCMvsNFD
deseq_DCMvsNFD_mag <- readRDS("../orignal_data/data/MAGNetApp/data/DGE/MAGNet_DCMvsNFD_GeneTonic.rds")$res_de

# filter
DEgenes_DCMvsNFD_mag <- rownames(subset(deseq_DCMvsNFD_mag, padj <= 0.05))
DEgenes_DCMvsNFD_nf <- rownames(subset(deseq_DCMvsNFD_nf, padj <= 0.05))

# compare gene lists
overlap_DCMvsNFD <- data.frame(genes = intersect(DEgenes_DCMvsNFD_mag, DEgenes_DCMvsNFD_nf))
overlap_DCMvsNFD$type <- "overlapping"
nf_DCMvsNFD <- data.frame(genes = setdiff(DEgenes_DCMvsNFD_nf,rownames(overlap_DCMvsNFD)))
nf_DCMvsNFD$type <- "nf-core"
mag_DCMvsNFD <- data.frame(genes = setdiff(DEgenes_DCMvsNFD_mag,rownames(overlap_DCMvsNFD)))
mag_DCMvsNFD$type <- "Magnetique"
DCMvsNFD <- rbind(overlap_DCMvsNFD, nf_DCMvsNFD, mag_DCMvsNFD)
paste0("Proportion overlapping DE genes DCMvsNFD: ", nrow(overlap_DCMvsNFD)/nrow(DEgenes_DCMvsNFD_mag))

[1] "Proportion overlapping DE genes DCMvsNFD: "

# venn
biovenn <- draw.venn(DEgenes_DCMvsNFD_mag, DEgenes_DCMvsNFD_nf, NULL, title="DCM vs NFD", nrtype="abs", subtitle = NULL, xtitle = "Magnetique", ytitle = "nf-core")

[1] "x total: 16562"
[1] "y total: 16722"
[1] "z total: 0"
[1] "x only: 2335"
[1] "y only: 2495"
[1] "z only: 0"
[1] "x-y total overlap: 14227"
[1] "x-z total overlap: 0"
[1] "y-z total overlap: 0"
[1] "x-y only overlap: 14227"
[1] "x-z only overlap: 0"
[1] "y-z only overlap: 0"
[1] "x-y-z overlap: 0"

ggplot(deseq_DCMvsNFD_merge_df, aes(x=log2FoldChange, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
    guides(color = "none") +
  xlim(c(-3,3)) +
  facet_wrap(~method)

ggplot(deseq_DCMvsNFD_merge_df, aes(x=baseMean, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
    guides(color = "none") +
  xlim(c(1,100)) +
  facet_wrap(~method)

HCM vs NFD

# load results tables
deseq_HCMvsNFD_nf <- readRDS("../DESeq2/harmonized-deseq/harmonized_DESeqResults.rds")$HCMvsNFD
deseq_HCMvsNFD_mag <- readRDS("../orignal_data/data/MAGNetApp/data/DGE/MAGNet_HCMvsNFD_GeneTonic.rds")$res_de

# filter
DEgenes_HCMvsNFD_mag <- rownames(subset(deseq_HCMvsNFD_mag, padj <= 0.05))
DEgenes_HCMvsNFD_nf <- rownames(subset(deseq_HCMvsNFD_nf, padj <= 0.05))

# compare gene lists
overlap_HCMvsNFD <- intersect(DEgenes_HCMvsNFD_mag, DEgenes_HCMvsNFD_nf)
paste0("Proportion overlapping DE genes HCMvsNFD: ", length(overlap_HCMvsNFD)/length(DEgenes_HCMvsNFD_nf))

[1] "Proportion overlapping DE genes HCMvsNFD: 0.828318108543794"

# venn
biovenn <- draw.venn(DEgenes_HCMvsNFD_mag, DEgenes_HCMvsNFD_nf, NULL, title="HCM vs NFD", nrtype="abs", subtitle = NULL, xtitle = "Magnetique", ytitle = "nf-core")

[1] "x total: 10607"
[1] "y total: 11166"
[1] "z total: 0"
[1] "x only: 1358"
[1] "y only: 1917"
[1] "z only: 0"
[1] "x-y total overlap: 9249"
[1] "x-z total overlap: 0"
[1] "y-z total overlap: 0"
[1] "x-y only overlap: 9249"
[1] "x-z only overlap: 0"
[1] "y-z only overlap: 0"
[1] "x-y-z overlap: 0"

length(intersect(DEgenes_HCMvsNFD_mag, DEgenes_HCMvsNFD_nf))

[1] 9249