figure1.Rmd

---
title: "The main figure1"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE)
knitr::opts_chunk$set(fig.width=12, fig.height=7)
```

```{r,message=FALSE,warning=FALSE}
#library(splatter)
library(Seurat)
library(speckle)
library(DCATS)
library(ggplot2)
library(tidyverse)
library(MCMCpack)
library(pROC)
library(patchwork)
library(ggpubr)
library(lattice)
```

```{r}
source("functionsV2.r")
options(future.globals.maxSize = 20000 * 1024^2) # 20G memory
```

```{r}
theme_set(theme_classic()+
            theme(panel.border = element_blank(),
          legend.key = element_blank(),
          axis.ticks = element_blank(),
          panel.grid = element_blank(),
          panel.grid.minor = element_blank(), 
          panel.grid.major = element_blank(),
          panel.background = element_blank(),
          legend.background = element_blank(),
          plot.background = element_rect(fill = "transparent",colour = NA)) + 
            theme(axis.text=element_text(size=12), axis.title=element_text(size=13))+
            theme(legend.title = element_text(size=13), #change legend title font size
                  legend.text = element_text(size=10)))
theme_set(theme_classic()+
            theme(panel.border = element_blank(),
          legend.key = element_blank(),
          axis.ticks = element_blank(),
          panel.grid = element_blank(),
          panel.grid.minor = element_blank(), 
          panel.grid.major = element_blank(),
          panel.background = element_blank(),
          legend.background = element_blank(),
          plot.background = element_rect(fill = "transparent",colour = NA)) + 
            theme(axis.text=element_text(size=15), axis.title=element_text(size=17))+
            theme(strip.text.x = element_text(size = 15), strip.text.y = element_text(size = 15)) +
            theme(plot.title = element_text(size = 20),plot.subtitle = element_text(size = 17)) +
            theme(legend.title = element_text(size=15), #change legend title font size
                  legend.text = element_text(size=15)) + 
  theme(plot.title = element_text(hjust = 0.5)))
```


### Figure 1 - B

```{r}
load("/storage/holab/linxy/DCATS/simulation/toy_simulation.RData")
prop_cond1M = numb_cond1M[2:nrow(numb_cond1M),]/2000
prop_cond2M = numb_cond2M[2:nrow(numb_cond2M),]/2000
prop_cond1BCM = numb_cond1BCM[2:nrow(numb_cond1BCM),]/2000
prop_cond2BCM = numb_cond2BCM[2:nrow(numb_cond2BCM),]/2000
prop_true1M = true_cond1M[2:nrow(true_cond1M),]/2000
prop_true2M = true_cond2M[2:nrow(true_cond2M),]/2000
cond1DF = rbind(prop_cond1M, prop_cond1BCM, prop_cond2M, prop_cond2BCM, prop_true1M, prop_true2M) %>% as.data.frame() %>% 
  mutate(state = c(rep("Observed cell proportions", nrow(numb_cond1M)-1), rep("Corrected cell proportions", nrow(numb_cond1BCM)-1), rep("Observed cell proportions", nrow(numb_cond2M)-1), rep("Corrected cell proportions", nrow(numb_cond2BCM)-1), rep("True cell proportions", nrow(prop_true1M) + nrow(prop_true2M))),
         cond = c(rep("Condition 1", nrow(numb_cond1M) + nrow(numb_cond1BCM) -2), rep("Condition 2", nrow(numb_cond2M) + nrow(numb_cond2BCM) -2), rep("Condition 1", nrow(prop_true1M)), rep("Condition 2", nrow(prop_true2M)))) %>% 
  mutate(state = factor(state, levels = c("True cell proportions", "Observed cell proportions", "Corrected cell proportions")))
colnames(cond1DF) = c("A", "B", "C", "state", "cond")
```

```{r}
## this is the figure we want to use in the main text
sim_prop = cond1DF %>% 
  filter(state != "Corrected cell proportions") %>% 
  pivot_longer(A:C, names_to = "cellType", values_to = "proportion") %>% 
  ggplot(aes(cellType, proportion)) + 
  geom_boxplot(aes(col = cellType)) +
  facet_grid(cond~state) + 
  theme(legend.position = "none") +
  xlab("cell type")
sim_prop
```

```{r}
## the supplementary figure
sim1_prop = cond1DF %>% 
  pivot_longer(A:C, names_to = "cellType", values_to = "proportion") %>% 
  ggplot(aes(cellType, proportion)) + 
  geom_boxplot(aes(col = cellType)) +
  facet_grid(cond~state) + 
  theme(legend.position = "none") +
  xlab("cell type")
sim1_prop
ggsave("./plot/sim1_proportion.png", bg = "transparent")
```

### Figure 1 - C

```{r}
evaluationDF = data.frame()
len = length(simulationDF_list)
#for (j in 1:8){
for (j in 1:len){
  simulationDF = simulationDF_list[[j]] %>% na.omit()
  method = colnames(simulationDF)[3:dim(simulationDF)[2]]
  numb_mthd = length(method)
  sensitivity = rep(NA, numb_mthd)
  specificity = rep(NA, numb_mthd)
  mcc = rep(NA, numb_mthd)
  auc = rep(NA, numb_mthd)
  prauc = rep(NA, numb_mthd)
  F1 = rep(NA, numb_mthd)
  
  truth = simulationDF$truth
  for (i in 3:dim(simulationDF)[2]){
    pred = simulationDF[, i]
    pred_res = ifelse(pred < 0.05, "P", "N")
    TP <- sum(pred_res=="P"&truth=="P")
    TN <- sum(pred_res=="N"&truth=="N")
    FP <- sum(pred_res=="P"&truth=="N")
    FN <- sum(pred_res=="N"&truth=="P")
    truthN = TN + FP
    truthP = FN + TP
    predP = TP + FP
    predN = FN + TN
    sensitivity[i-2] = TP/truthP
    specificity[i-2] = TN/truthN
    mcc[i-2] = (TP*TN-FP*FN)/sqrt(predP*truthP*truthN*predN)
    auc[i-2] = getROC(truth, pred)$auc
    prauc[i-2] = getPRC(truth, pred)$prauc
    F1 = 2*TP/(2*TP+FP+FN)
    
    res = data.frame(trial = as.character(j), method = method, sensitivity = sensitivity, specificity = specificity, mcc = mcc, auc = auc, prauc = prauc, F1 = F1)
    }
    evaluationDF = rbind(evaluationDF, res)
}
```

```{r}
sim1_main = evaluationDF %>% 
  tidyr::separate(method, c("method", "condition"), sep = "_") %>% 
  dplyr::select(-sensitivity, -specificity) %>% 
  filter(condition != "test") %>% 
  filter(method %in% c("estPhi", "fisher")) %>% 
  mutate(method = ifelse(method == "fisher", "Fisher", method),
         method = ifelse(method == "estPhi", "DCATS", method)) %>% 
  ggplot(aes(x = condition, y = mcc, color = method)) + geom_boxplot() +
  scale_color_manual(values=c("#F8766D", "#00BFC4")) +
  theme(axis.title.x = element_blank()) + 
  theme(legend.position = "top") +
  ylab("MCC") +
  scale_x_discrete(labels=c("null" = "No bias correction", "withBC" = "With bias correction")) 
sim1_main
ggsave("./plot/sim1_main.png", bg = "transparent")
```

Simulation 1: Supplementary Plots

```{r}
evaluationDF %>% 
  tidyr::separate(method, c("method", "condition"), sep = "_") %>% 
  dplyr::select(-sensitivity, -specificity) %>% 
  filter(condition != "test") %>% 
  mutate(method = ifelse(method == "fisher", "Fisher", method),
         method = factor(method, levels = c("wtoPhi", "adjPhi", "avrgPhi", "estPhi", "Fisher"))) %>% 
  ggplot(aes(x = condition, y = auc, color = method)) + geom_boxplot() +
  theme(axis.title.x = element_blank()) + 
  ylab("AUC") +
  scale_x_discrete(labels=c("null" = "No bias correction", "withBC" = "With bias correction"))
ggsave("./plot/sim1_allAUC.png", bg = "transparent")
```

```{r}
evaluationDF %>% 
  tidyr::separate(method, c("method", "condition"), sep = "_") %>% 
  filter(condition != "test") %>% 
  dplyr::select(-auc, -prauc) %>% 
  pivot_longer(
    sensitivity:mcc,
    names_to = "statistics", 
    values_to = "value") %>% 
  mutate(condition = ifelse(condition == "null", "noBC", condition)) %>%
  ggplot(aes(x = condition, y = value, color = method)) + geom_boxplot() + facet_grid(.~statistics) +
  ylab("Value") +
    theme(axis.title.x = element_blank())
ggsave("./plot/sim1_statistics.png", bg = "transparent")
```

```{r}
supA1 = phiDF %>% 
  ggplot(aes(x = eachPhi1, y = estPhi1)) +
  geom_point(aes(col = cluster), size = 0.5) + 
  theme(legend.position = "none") +
  xlab("eachPhi (no BC)") +
  ylab("estPhi (no BC)")
supA2 = phiDF %>% 
  ggplot(aes(x = eachPhi2, y = estPhi2)) +
  geom_point(aes(col = cluster), size = 0.5) +
  theme(legend.position = "none") +
  xlab("eachPhi (with BC)") +
  ylab("estPhi (with BC)")
supA3 = phiDF %>% 
  ggplot(aes(x = eachPhi1, y = eachPhi2)) +
  geom_point(aes(col = cluster), size = 0.5) +
  theme(legend.position = "none") +
  xlab("eachPhi (no BC)") +
  ylab("eachPhi (with BC)")
supA4 = phiDF %>% 
  ggplot(aes(x = estPhi1, y = estPhi2)) +
  geom_point(aes(col = cluster), size = 0.5) +
  theme(legend.position = "bottom") +
  xlab("estPhi (no BC)") +
  ylab("estPhi (with BC)")
(supA1 + supA2)/(supA3 + supA4)
ggsave("./plot/sim1_phi.png", bg = "transparent", width = 8, height = 5)
```

```{r}
sim_prop + sim1_main
ggsave("./plot/version2/Fig1.png", bg = "transparent", width = 12, height = 6)
ggsave("./plot/version2/Fig1.pdf", bg = "transparent", width = 12, height = 6)
```