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

```{r}
library(ggplot2)
library(ggcorrplot)
library(BioVenn)
library(dplyr)
```

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

## Correlation of counts
### Example for one sample match
```{r}
# 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'
```

```{r}
# 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'
```


```{r}
# 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()
```

```{r}
# correlation
cor(merge$magnetique, merge$`nf-core`)
```
### All vs all
```{r}
# 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))]

```

```{r}
# correlation matrix
cor_mat <- cor(data_mag_sub, data_nf_sub)
```
```{r}
# plot first 30 samples
ggcorrplot(cor_mat) + scale_fill_gradient2(limit = c(0.5,1), low = "blue", high =  "red", mid = "white", midpoint = 0.75)
```
### Using normalised count data
```{r}
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'
```

```{r}
# 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")
```
```{r}
cor(merge_norm$magnetique, merge_norm$`nf-core`)
```
### All vs all - norm
```{r}
# 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)
```
```{r}
# 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)
```

### PCA plots
```{r}
# 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()
```
```{r}
# 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
```{r}
# 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

```

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

```{r}
# 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))
```
```{r}
# venn
biovenn <- draw.venn(DEgenes_DCMvsHCM_mag, DEgenes_DCMvsHCM_nf, NULL, title="DCM vs HCM", nrtype="abs", subtitle = NULL, xtitle = "Magnetique", ytitle = "nf-core")
```
```{r}
# distribution of p-values, l2FCs and counts by overlapping or unique
# plots of values are performed for each gene set against both DESeq2 results from Magnetique and nf-core pipelines
deseq_DCMvsHCM_mag_df <- as.data.frame(deseq_DCMvsHCM_mag[,c("baseMean","log2FoldChange","padj")]) %>%
  tibble::rownames_to_column(., var = "genes")
deseq_DCMvsHCM_nf_df <- as.data.frame(deseq_DCMvsHCM_nf[,c("baseMean","log2FoldChange","padj")]) %>%
  tibble::rownames_to_column(., var = "genes")
deseq_DCMvsHCM_mag_df <- left_join(deseq_DCMvsHCM_mag_df, DCMvsHCM, by="genes") %>%
  replace(is.na(.), "none")
deseq_DCMvsHCM_nf_df <- left_join(deseq_DCMvsHCM_nf_df, DCMvsHCM, by="genes") %>%
  replace(is.na(.), "none")
deseq_DCMvsHCM_mag_df$method <- "Magnetique"
deseq_DCMvsHCM_nf_df$method <- "nf-core"
deseq_DCMvsHCM_merge_df <- rbind(deseq_DCMvsHCM_mag_df, deseq_DCMvsHCM_nf_df)

ggplot(deseq_DCMvsHCM_merge_df, aes(x=padj, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
  #scale_fill_discrete(labels = c("Distinct", "Overlapping")) +
    guides(color = "none") +
  facet_wrap(~method) +
  geom_vline(xintercept = 0.05, linetype="dotted")
```
```{r}
# same as above but exclude overlapping
ggplot(subset(deseq_DCMvsHCM_merge_df, type != "overlapping"), aes(x=padj, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
  #scale_fill_discrete(labels = c("Distinct", "Overlapping")) +
    guides(color = "none") +
  facet_wrap(~method) +
  xlim(c(0,0.2)) +
    geom_vline(xintercept = 0.05, linetype="dotted")
```
```{r}
# 
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)
```

```{r}
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
```{r}
# 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

```

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

```{r}
# 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))
```

```{r}
# venn
biovenn <- draw.venn(DEgenes_DCMvsNFD_mag, DEgenes_DCMvsNFD_nf, NULL, title="DCM vs NFD", nrtype="abs", subtitle = NULL, xtitle = "Magnetique", ytitle = "nf-core")
```
```{r}
# distribution of p-values and l2FCs by overlapping or not
deseq_DCMvsNFD_mag_df <- as.data.frame(deseq_DCMvsNFD_mag[,c("baseMean","log2FoldChange","padj")]) %>%
  tibble::rownames_to_column(., var = "genes")
deseq_DCMvsNFD_nf_df <- as.data.frame(deseq_DCMvsNFD_nf[,c("baseMean","log2FoldChange","padj")]) %>%
  tibble::rownames_to_column(., var = "genes")
deseq_DCMvsNFD_mag_df <-  left_join(deseq_DCMvsNFD_mag_df, DCMvsNFD, by="genes") %>%
  replace(is.na(.), "none")
deseq_DCMvsNFD_nf_df <- left_join(deseq_DCMvsNFD_nf_df, DCMvsNFD, by="genes") %>%
  replace(is.na(.), "none")
deseq_DCMvsNFD_mag_df$method <- "Magnetique"
deseq_DCMvsNFD_nf_df$method <- "nf-core"
deseq_DCMvsNFD_merge_df <- rbind(deseq_DCMvsNFD_mag_df, deseq_DCMvsNFD_nf_df)

ggplot(deseq_DCMvsNFD_merge_df, aes(x=padj, group=type, color=type, fill=type)) + 
  geom_density(alpha=0.6) + 
  theme_bw() +
    guides(color = "none") +
  facet_wrap(~method)
  
```
```{r}
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)
```

```{r}
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
```{r}
# 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

```

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

```{r}
# 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))
```

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

```{r}
length(intersect(DEgenes_HCMvsNFD_mag, DEgenes_HCMvsNFD_nf))
```