simulation_extreme.Rmd

---
title: 'Simulation: extreme cases'
output: html_document
---

```{r}
library(Seurat)
library(SeuratWrappers)
library(speckle)
library(DCATS)
library(ggplot2)
library(tidyverse)
library(diffcyt)
library(MCMCpack)
#library(scdney)
library(scDC)
library(tidymodels)  ## for DCATS
library(miloR)
library(SingleCellExperiment)
```

## Extreme case1

The first case is all DA cell types increase.

```{r}
## simulation setting
cluster_num = 8  # numbers of clusters
setresolu = 0.6
rep1 = 3
rep2 = 3
simulation_times = 30
concentration = 100
sample_size1 = 3000
sample_size2 = 3000
cell_pool = "splatter"  # cell pool used, can be selected from ("SPAR", "splatter", "realWorld")
more_negative = ""
countC1 = rep(8, 8)
countC2 = c(rep(8,3),16, rep(8,3),16)
probC1 = countC1/sum(countC1)
probC2 = countC2/sum(countC2)
truthRes = c("N", "N", "N", "P", "N", "N", "N", "P")
print(probC1)
print(probC2)
file_name = "/storage/holab/linxy/DCATS/simulation/current/two_increase.RData"
```

## extreme case2

```{r}
## simulation setting
cluster_num = 8  # numbers of clusters
setresolu = 0.6
rep1 = 3
rep2 = 3
simulation_times = 30
concentration = 100
sample_size1 = 3000
sample_size2 = 3000
cell_pool = "splatter"  # cell pool used, can be selected from ("SPAR", "splatter", "realWorld")
more_negative = ""
#countC1 = c(rep(8,7),2)
#countC2 = c(rep(8,7), 40)
countC1 = c(2, rep(8,7))
countC2 = c(40, rep(8,7))
probC1 = countC1/sum(countC1)
probC2 = countC2/sum(countC2)
#truthRes = c("N", "N", "N", "N", "N", "N", "N", "P")
truthRes = c("P", "N", "N", "N", "N", "N", "N", "N")
print(probC1) %>% round(3)
print(probC2) %>% round(3)
file_name = "/storage/holab/linxy/DCATS/simulation/current/one_increase.RData"
```

## extreme case 3

```{r}
## simulation setting
cluster_num = 8  # numbers of clusters
setresolu = 0.6
rep1 = 3
rep2 = 3
simulation_times = 30
concentration = 100
sample_size1 = 3000
sample_size2 = 3000
cell_pool = "splatter"  # cell pool used, can be selected from ("SPAR", "splatter", "realWorld")
more_negative = ""
countC1 = rep(8, 8)
countC2 = c(16, 8, 8, 8, 12, rep(8,3))
probC1 = countC1/sum(countC1)
probC2 = countC2/sum(countC2)
truthRes = c("P", "N", "N", "N", "P", "N", "N", "N")
print(probC1)
print(probC2)
file_name = "/storage/holab/linxy/DCATS/simulation/current/two_increase_notEqual.RData"
```

```{r}
## load data and function
#load("/storage/holab/linxy/DCATS/cells_pool8_splatter.RData")
load("/storage/holab/linxy/DCATS/cells_pool8_splatter_large.RData") ## for extreme case2
source("functionsV2.r")
source("glm.R")     # for scDC
options(future.globals.maxSize = 15000 * 1024^2)
```

### Simulate count only

```{r}
## only simulate count
resDF = data.frame()
for (i in 1:50){
  simulation = simulator_base(rep(3000,3), rep(3000,3), probC1*100, probC2*100, diag(nrow = 8, ncol = 8))
  sim_design = data.frame(condition = c(rep("g1", rep1), rep("g2", rep2)))
  print(simulation$numb_cond1/rowSums(simulation$numb_cond1))
  print(simulation$numb_cond2/rowSums(simulation$numb_cond2))
  sim_count = rbind(simulation$numb_cond1, simulation$numb_cond2)
  phi = getPhi(sim_count, sim_design)
  estPhi_null = dcats_GLM(sim_count, sim_design, fix_phi = phi)
  ref_order = detect_reference(sim_count, sim_design, fix_phi = phi)$celltype
  print(ref_order)
  ref_res3 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = ref_order[1:3])
  ref_res2 = dcats_new(sim_count, sim_design, reference = ref_order[1:2])
  ref_res1 = dcats_new(sim_count, sim_design, reference = ref_order[1])
  org_res = dcats_GLM(sim_count, sim_design)
  truth_res3 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = c(5,6,7))
  truth_res2 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = c(5,6))
  truth_res1 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = 6)
  sub_res = data.frame(truth = truthRes, ref_pval1 = as.vector(ref_res1$LRT_pvals), ref_pval2 = as.vector(ref_res2$LRT_pvals), ref_pval3 = as.vector(ref_res3$LRT_pvals), org_pval = as.vector(org_res$LRT_pvals), truth_pval1 = as.vector(truth_res1$LRT_pvals), truth_pval2 = as.vector(truth_res2$LRT_pvals), truth_pval3 = as.vector(truth_res3$LRT_pvals))
  resDF = rbind(resDF, sub_res)
}
```

Using statistics to summary

```{r}
numb_mthd = 7
mcc = rep(NA, numb_mthd)
auc = rep(NA, numb_mthd)
prauc = rep(NA, numb_mthd)
sensitivity = rep(NA, numb_mthd)
specificity = rep(NA, numb_mthd)
precision = rep(NA, numb_mthd)
F1 = rep(NA, numb_mthd)
truth = resDF$truth
for (i in 2:(dim(resDF)[2])){
  pred = resDF[,i]
  pred_res = ifelse(pred < 0.1, "P", "N")
  auc[i-1] = getROC(resDF$truth, pred)$auc
  prauc[i-1] = getPRC(resDF$truth, pred)$prauc
  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-1] = TP/truthP
  specificity[i-1] = TN/truthN
  precision[i-1] = TP/predP
  #mcc[i-2] = sqrt(sensitivity[i]*specificity[i]*precision[i]*(TN/(TN+FN))) - sqrt((1-sensitivity[i])*(1-specificity[i])*(FN/predN)*(FP/predP))
  mcc[i-1] = (TP*TN-FP*FN)/sqrt(predP*truthP*truthN*predN)
  F1[i-1] = 2*TP/(2*TP+FP+FN)
}
res = data.frame(method = c("ref1", "ref2", "ref3", "org", "truth1", "truth2", "truth3"), mcc = mcc, auc = auc, prauc = prauc, sensitivity = sensitivity, specificity = specificity, precision = precision, F1 = F1) %>% 
  arrange(desc(auc))
```

## Add gene expression information

```{r}
## add gene expression 
set.seed(123)
true_countDF = data.frame()
seurat_countDF = data.frame()
knn_matrixL = list()
true_matrixL = list()
svm_matrixL = list()
reference_orderL = list()
clusterRes = data.frame()
timeDF = data.frame()
for (idx in 1:simulation_times) {
  print(str_c("idx: ", as.character(idx)))
  timeL = rep(NA, 8) # DCATS w/wto similarity matrix, fisher, speckle, diffcyt, scDC
  ## simulation
  simulation = simulator_noInt(totals1 = runif(rep1, sample_size1, sample_size2), totals2 = runif(rep2, sample_size1, sample_size2), probC1, probC2, setresolu, sim_mat)
  if (is.na(simulation)[1] == TRUE){
    next
  }
  true_count = rbind(cond1 = simulation$numb_cond1, cond2 = simulation$numb_cond2) %>% as.data.frame() %>% mutate(sim = idx)
  true_countDF = rbind(true_countDF, true_count)
  numb_cond1 = simulation$dfRes %>% 
    filter(condition == "Cond1") %>%
    group_by(clusterRes, batch) %>% 
    summarise(n = n()) %>% 
    pivot_wider(names_from = "clusterRes", values_from = "n") %>% 
    column_to_rownames(var = "batch")
  numb_cond2 = simulation$dfRes %>% 
    filter(condition == "Cond2") %>%
    group_by(clusterRes, batch) %>% 
    summarise(n = n()) %>% 
    pivot_wider(names_from = "clusterRes", values_from = "n") %>% 
    column_to_rownames(var = "batch")
  numb_cond1[is.na(numb_cond1)] <- 0
  numb_cond2[is.na(numb_cond2)] <- 0
  conf.mat <- table(simulation$trueLabels, Idents(simulation$integratedSamples))
  true.conf <- conf.mat/rowSums(conf.mat)
  library(clue)
  order2 = as.vector(solve_LSAP(true.conf, maximum = TRUE))
  L_sum1 = colSums(numb_cond1)
  L_sum2 = colSums(numb_cond2)
  decider = sum(L_sum1 > 10)+sum(L_sum2 > 10)
  if (decider != cluster_num*2){
    next
  } # if one cluster has least than 10 cells in one condition, then won't continue
  seurat_count = rbind(numb_cond1[,order2], numb_cond2[,order2]) %>% mutate(sim = idx)
  seurat_countDF = rbind(seurat_countDF, seurat_count)
  ## Get matrices
  ### KNN matrix
  graphs = simulation$integratedSamples@graphs$RNA_snn
  labels = Idents(simulation$integratedSamples)
  knn_mat = knn_simMat(graphs, labels)
  order1 = colnames(numb_cond1)
  simil_matK = knn_mat[order1, order1]
  knn_matrixL[[idx]] = simil_matK
  ### True matrix
  simil_matT = true.conf[,order2]
  true_matrixL[[idx]] = simil_matT
  ### Uniform matrix
  simil_matU = create_simMat(cluster_num, confuse_rate=0.1)
  ### svm&rf
  mlDF = simulation$integratedSamples@reductions$pca@cell.embeddings %>% 
    as.data.frame() %>% 
    tibble::rownames_to_column("cellID") %>% 
    merge(simulation$dfRes, by = "cellID") %>% 
    dplyr::select(-cellID)
  set.seed(123)
  simil_matSVM = svm_simMat(mlDF)
  svm_matrixL[[idx]] = simil_matSVM
  
  ## Test
  ### DCATS---betabinLRT
  sim_count = rbind(numb_cond1, numb_cond2) %>% as.matrix()
  sim_design = data.frame(condition = c(rep("g1", rep1), rep("g2", rep2)))
  t1start = Sys.time()
  betabin_null = dcats_GLM(sim_count, sim_design)
  timeL[1] = Sys.time() - t1start
  ###  DCATS---betabinLRT with Bias correction from similarity matrix
  phi = getPhi(sim_count, sim_design)
  t2start = Sys.time()
  estPhi_null = dcats_GLM(sim_count, sim_design, fix_phi = phi)
  timeL[2] = Sys.time() - t2start
  ## KNN
  wtoPhi_emK = dcats_GLM(sim_count, sim_design, simil_matK)
  t3start = Sys.time()
  estPhi_emK = dcats_GLM(sim_count, sim_design, simil_matK, fix_phi = phi)
  timeL[3] = Sys.time() - t3start
  ## Uniform
  wtoPhi_emU = dcats_GLM(sim_count, sim_design, simil_matU)
  estPhi_emU = dcats_GLM(sim_count, sim_design, simil_matU, fix_phi = phi)
  ## True
  wtoPhi_emT = dcats_GLM(sim_count, sim_design, simil_matT)
  estPhi_emT = dcats_GLM(sim_count, sim_design, simil_matT, fix_phi = phi)
  ## svm
  wtoPhi_emSVM = dcats_GLM(sim_count, sim_design, simil_matSVM)
  estPhi_emSVM = dcats_GLM(sim_count, sim_design, simil_matSVM, fix_phi = phi)
  ## add reference
  ref_order = detect_reference(sim_count, sim_design, fix_phi = phi)$celltype
  knnref_order = detect_reference(sim_count, sim_design, simil_matK, fix_phi = phi)$celltype
  svmref_order = detect_reference(sim_count, sim_design, simil_matSVM, fix_phi = phi)$celltype
  
  reference_orderL[[idx]] = data.frame(cluster = colnames(numb_cond1[,order2]), originCluster = 1:8, truth = truthRes) %>% 
    merge(data.frame(cluster = ref_order, ref_orderidx = 1:8)) %>% 
    merge(data.frame(cluster = knnref_order, knnref_orderidx = 1:8)) %>% 
    merge(data.frame(cluster = svmref_order, svmref_orderidx = 1:8))
    
  print(reference_orderL[[idx]])
  emSVM_refSVM3 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = svmref_order[1:3])
  emSVM_refSVM2 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = svmref_order[1:2])
  emK_refK3 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = knnref_order[1:3])
  emK_refK2 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = knnref_order[1:2])
  emSVM_org3 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = ref_order[1:3])
  emSVM_org2 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = ref_order[1:2])
  emK_org3 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = ref_order[1:3])
  emK_org2 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = ref_order[1:2])
  emSVM_truth3 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = colnames(numb_cond1[,order2])[c(5, 6,7)])
  emSVM_truth2 = dcats_new(sim_count, sim_design, simil_matSVM, fix_phi = phi, reference = colnames(numb_cond1[,order2])[c(6,7)])
  emK_truth3 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = colnames(numb_cond1[,order2])[c(5, 6,7)])
  emK_truth2 = dcats_new(sim_count, sim_design, simil_matK, fix_phi = phi, reference = colnames(numb_cond1[,order2])[c(6,7)])
  
  ## Fisher's exact test
  t4start = Sys.time()
  fisher_pvals = getFisher(as.matrix(numb_cond1), as.matrix(numb_cond2))
  timeL[4] = Sys.time() - t4start
  ## speckle
  t5start = Sys.time()
  speckleRes = propeller(clusters = simulation$dfRes$clusterRes, sample = simulation$dfRes$batch, group = simulation$dfRes$condition) %>% 
    dplyr::rename(cluster = BaselineProp.clusters, speckle_pvals = FDR) %>%
    dplyr::select(cluster, speckle_pvals)
  timeL[5] = Sys.time() - t5start
  speckleRes = propeller(clusters = simulation$dfRes$clusterRes, sample = simulation$dfRes$batch, group = simulation$dfRes$condition) %>% 
    dplyr::rename(cluster = BaselineProp.clusters, speckle_pvals = FDR) %>%
    dplyr::select(cluster, speckle_pvals)
  ## diffcyt
  t6start = Sys.time()
  diffcytP = getDiffcyt(numb_cond1, numb_cond2, simulation$dfRes)
  timeL[6] = Sys.time() - t6start
  ## scDC
  t7start = Sys.time()
  res_scDC <- scDC_noClustering(cellTypes = simulation$dfRes$clusterRes, simulation$dfRes$batch, calCI = TRUE, calCI_method = c("percentile", "BCa", "multinom"),nboot = 1000, verbose = FALSE)
  res_GLM <- fitGLM(res_scDC, c(rep("cond1",rep1*cluster_num),rep("cond2",rep2*cluster_num)), pairwise = FALSE, fixed_only = TRUE, verbose = FALSE)
  timeL[7] = Sys.time() - t7start
  scDCRes_temp = summary(res_GLM$pool_res_fixed)
  scDCRes = scDCRes_temp[c(cluster_num+1,(dim(scDCRes_temp)[1]-cluster_num+2):dim(scDCRes_temp)[1]),]
  ## milo
  t8start = Sys.time()
  milo_res = getMilo(simulation$integratedSamples)
  timeL[8] = Sys.time() - t8start
  miloDF_temp = milo_res$da_results %>% 
    filter(ident_fraction > 0.8) %>%
    mutate(direct = ifelse( logFC > 0  & SpatialFDR < 0.1, "milo_more", "neutral"),
           direct = ifelse( logFC < 0  & SpatialFDR < 0.1, "milo_less", direct)) %>% 
    group_by(ident, direct) %>% 
    summarise(n = n()) %>% 
    pivot_wider(names_from = direct, values_from = n) %>% 
    replace(is.na(.), 0) %>% 
    dplyr::rename(cluster = ident)
  if ("milo_less" %!in% colnames(miloDF_temp)) {miloDF_temp$milo_less = rep(0, nrow(miloDF_temp))}
  if ("milo_more" %!in% colnames(miloDF_temp)) {miloDF_temp$milo_more = rep(0, nrow(miloDF_temp))}
  miloDF = miloDF_temp %>% 
    dplyr::select(cluster, milo_less, neutral, milo_more)
  ## results
  cluster_map = data.frame(cluster = colnames(numb_cond1[,order2]), truth = truthRes) %>% 
    arrange(cluster)
  clusterDF = cluster_map %>% 
    mutate(betabin_null_pvals = as.vector(betabin_null$LRT_pvals),
           estPhi_null_pvals = as.vector(estPhi_null$LRT_pvals),
           wtoPhi_emSVM_pvals = as.vector(wtoPhi_emSVM$LRT_pvals),
           estPhi_emSVM_pvals = as.vector(estPhi_emSVM$LRT_pvals),
           wtoPhi_emT_pvals = as.vector(wtoPhi_emT$LRT_pvals),
           estPhi_emT_pvals = as.vector(estPhi_emT$LRT_pvals),
           wtoPhi_emU_pvals = as.vector(wtoPhi_emU$LRT_pvals),
           estPhi_emU_pvals = as.vector(estPhi_emU$LRT_pvals),
           wtoPhi_emK_pvals = as.vector(wtoPhi_emK$LRT_pvals),
           estPhi_emK_pvals = as.vector(estPhi_emK$LRT_pvals),
           emSVM_refSVM3_pvals = as.vector(emSVM_refSVM3$LRT_pvals),
           emSVM_refSVM2_pvals = as.vector(emSVM_refSVM2$LRT_pvals),
           emK_refK3_pvals = as.vector(emK_refK3$LRT_pvals),
           emK_refK2_pvals = as.vector(emK_refK2$LRT_pvals),
           emSVM_org3_pvals = as.vector(emSVM_org3$LRT_pvals),
           emSVM_org2_pvals = as.vector(emSVM_org2$LRT_pvals),
           emK_org3_pvals = as.vector(emK_org3$LRT_pvals),
           emK_org2_pvals = as.vector(emK_org2$LRT_pvals),
           #emSVM_truth3_pvals = as.vector(emSVM_truth3$LRT_pvals),
           emSVM_truth2_pvals = as.vector(emSVM_truth2$LRT_pvals),
           #emK_truth3_pvals = as.vector(emK_truth3$LRT_pvals),
           emK_truth2_pvals = as.vector(emK_truth2$LRT_pvals),
           fisher_pvals = fisher_pvals,
           scDC_pvals = scDCRes$p.value) %>%
    merge(speckleRes, by = "cluster") %>% 
    merge(diffcytP, by = "cluster") %>% 
    left_join(miloDF, by = "cluster") %>% 
    replace(is.na(.), 0)
  clusterRes = rbind(clusterRes, clusterDF)
  timeDF = rbind(timeDF, timeL)
}

resL = list(clusterRes = clusterRes, timeDF = timeDF, svm_matrixL = svm_matrixL, true_matrixL = true_matrixL, knn_matrixL = knn_matrixL, reference_orderL = reference_orderL, seurat_countDF = seurat_countDF, true_countDF = true_countDF)
#save(resL, file = file_name)
```

Using statistics to summary

```{r}
## calculate statistics of "three_increase.RData" results
#load(file_name)
#clusterRes = resL$clusterRes
methods = colnames(clusterRes)[colnames(clusterRes) %>% str_detect("_pvals")] %>% str_remove("_pvals")
#methods = colnames(clusterRes)[colnames(clusterRes) %>% str_detect("_pvals")] %>% str_remove("_pvals")
methods = c(methods, "milo")
numb_mthd = length(methods)
mcc = rep(NA, numb_mthd)
auc = rep(NA, numb_mthd)
prauc = rep(NA, numb_mthd)
sensitivity = rep(NA, numb_mthd)
specificity = rep(NA, numb_mthd)
precision = rep(NA, numb_mthd)
F1 = rep(NA, numb_mthd)
truth = clusterRes$truth
clusterRes = clusterRes %>% 
  mutate(milo_main = ifelse(milo_less > milo_more, milo_less, milo_more),
         milo_pct = milo_main/(milo_less + neutral + milo_more),
         milo_diff = abs(milo_less - milo_more)) %>% 
  dplyr::select(-milo_less, -milo_more, - neutral, -milo_main)

for (i in 3:(dim(clusterRes)[2]-1)){
  pred = clusterRes[,i]
  if (colnames(clusterRes)[i] == "milo_pct") {
    pred_res = ifelse(clusterRes$milo_diff > 0, "P", "N")
    auc[i-2] = getROC(clusterRes$truth, 1-pred)$auc
    prauc[i-2] = getPRC(clusterRes$truth, 1-pred)$prauc}
  else {
    pred_res = ifelse(pred < 0.1, "P", "N")
    auc[i-2] = getROC(clusterRes$truth, pred)$auc
    prauc[i-2] = getPRC(clusterRes$truth, pred)$prauc}
  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
  precision[i-2] = TP/predP
  #mcc[i-2] = sqrt(sensitivity[i]*specificity[i]*precision[i]*(TN/(TN+FN))) - sqrt((1-sensitivity[i])*(1-specificity[i])*(FN/predN)*(FP/predP))
  mcc[i-2] = (TP*TN-FP*FN)/sqrt(predP*truthP*truthN*predN)
  F1[i-2] = 2*TP/(2*TP+FP+FN)
}
res = data.frame(methods = methods, mcc = mcc, auc = auc, prauc = prauc, sensitivity = sensitivity, specificity = specificity, precision = precision, F1 = F1) %>% 
  arrange(desc(mcc))
  #filter(methods == "wtoPhi_fullK" | methods == "diffcyt" | methods == "speckle" | methods == "fisher" | methods == "scDC")
print(res)
```

```{r}
orderDF = data.frame()
for (idx in 1:simulation_times){
  temp = reference_orderL[[idx]]
  orderDF = rbind(orderDF, temp)
}
```

```{r}
orderDF %>% 
  pivot_longer(ref_orderidx:svmref_orderidx, names_to = "method", values_to = "idx") %>% 
  mutate(originCluster = as.factor(originCluster)) %>% 
  ggplot(aes(x = originCluster, y = idx, color = method)) +
  geom_boxplot()
```




```{r}
seurat_countDF %>% 
  rownames_to_column("condition") %>% 
  mutate(condition = str_sub(condition, 1,5)) %>% 
  pivot_longer(D:B, names_to = "cellType", values_to = "count") %>% 
  mutate(proportion = count/3000) %>% 
  ggplot(aes(cellType, proportion)) +
  geom_boxplot(aes(col = condition)) +
  theme(legend.position = "top") + 
  xlab("Cell type")
```

```{r}
true_countDF %>% 
  group_by(sim) %>% 
  mutate(condition = c(rep("cond1", 3), rep("cond2", 3))) %>% 
  #rownames_to_column("condition") %>% 
  #mutate(condition = str_sub(condition, 1,5)) %>% 
  pivot_longer(V1:V8, names_to = "cellType", values_to = "count") %>% 
  mutate(proportion = count/3000) %>% 
  ggplot(aes(cellType, proportion)) +
  geom_boxplot(aes(col = condition)) +
  theme(legend.position = "top") + 
  xlab("Cell type") +
  ggtitle("True count")
```

```{r}
print(reference_orderL[[idx]])
head(clusterRes)
```

## Test effect size

```{r}
## simulation setting
cluster_num = 8  # numbers of clusters
setresolu = 0.6
rep1 = 3
rep2 = 3
simulation_times = 30
concentration = 100
sample_size1 = 3000
sample_size2 = 3000
cell_pool = "splatter"  # cell pool used, can be selected from ("SPAR", "splatter", "realWorld")
more_negative = ""
countC1 = rep(8, 8)
#countC2 = c(rep(10,8),12,12)
probC1 = countC1/sum(countC1)
#probC2 = countC2/sum(countC2)
truthRes = c("N", "N", "N", "N", "N", "N", "P", "P")
#print(probC2)
file_name = "/Users/linxy29/Documents/Data/DCATS/simulation/current/two_increase_toy.RData"
```

```{r}
## load data and function
load("/Users/linxy29/Documents/Data/DCATS/cells_pool8_splatter.RData")
#load("/storage/holab/linxy/DCATS/cells_pool8_splatter_large.RData") ## for extreme case2
source("functionsV2.r")
source("glm.R")     # for scDC
options(future.globals.maxSize = 15000 * 1024^2)
```

### Simulate count only

```{r}
## only simulate count
resL = list()
valueV = seq(10, 20, 2)
for (v in 1:length(valueV)){
  countC2 = c(rep(8,6),valueV[v],valueV[v])
  probC2 = countC2/sum(countC2)
  print(probC2)
  for (i in 1:50){
  simulation = simulator_base(rep(3000,3), rep(3000,3), probC1*100, probC2*100, diag(nrow = 8, ncol = 8))
  sim_design = data.frame(condition = c(rep("g1", rep1), rep("g2", rep2)))
  #print(simulation$numb_cond1/rowSums(simulation$numb_cond1))
  #print(simulation$numb_cond2/rowSums(simulation$numb_cond2))
  sim_count = rbind(simulation$numb_cond1, simulation$numb_cond2)
  phi = getPhi(sim_count, sim_design)
  estPhi_null = dcats_GLM(sim_count, sim_design, fix_phi = phi)
  ref_order = detect_reference(sim_count, sim_design, fix_phi = phi)$celltype
  #print(ref_order)
  ref_res3 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = ref_order[1:3])
  ref_res2 = dcats_new(sim_count, sim_design, reference = ref_order[1:2])
  ref_res1 = dcats_new(sim_count, sim_design, reference = ref_order[1])
  org_res = dcats_GLM(sim_count, sim_design)
  truth_res3 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = c(4,5,6))
  truth_res2 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = c(5,6))
  truth_res1 = dcats_new(sim_count, sim_design, fix_phi = phi, reference = 6)
  resDF = data.frame(truth = truthRes, ref_pval1 = as.vector(ref_res1$LRT_pvals), ref_pval2 = as.vector(ref_res2$LRT_pvals), ref_pval3 = as.vector(ref_res3$LRT_pvals), org_pval = as.vector(org_res$LRT_pvals), truth_pval1 = as.vector(truth_res1$LRT_pvals), truth_pval2 = as.vector(truth_res2$LRT_pvals), truth_pval3 = as.vector(truth_res3$LRT_pvals))
  }
  resL[[v]] = resDF
}

```

```{r}
statisticsDF = data.frame()
numb_mthd = 7
for (idx in 1:length(resL)){
  resDF = resL[[idx]]
  mcc = rep(NA, numb_mthd)
  auc = rep(NA, numb_mthd)
  prauc = rep(NA, numb_mthd)
  sensitivity = rep(NA, numb_mthd)
  specificity = rep(NA, numb_mthd)
  precision = rep(NA, numb_mthd)
  F1 = rep(NA, numb_mthd)
  truth = resDF$truth
  for (i in 2:(dim(resDF)[2])){
    pred = resDF[,i]
    pred_res = ifelse(pred < 0.1, "P", "N")
    auc[i-1] = getROC(resDF$truth, pred)$auc
    prauc[i-1] = getPRC(resDF$truth, pred)$prauc
    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-1] = TP/truthP
    specificity[i-1] = TN/truthN
    precision[i-1] = TP/predP
    mcc[i-1] = (TP*TN-FP*FN)/sqrt(predP*truthP*truthN*predN)
    F1[i-1] = 2*TP/(2*TP+FP+FN)
    }
  res = data.frame(method = c("ref1", "ref2", "ref3", "org", "truth1", "truth2", "truth3"), mcc = mcc, auc = auc, prauc = prauc, sensitivity = sensitivity, specificity = specificity, precision = precision, F1 = F1) %>% 
  arrange(desc(mcc)) %>% 
    mutate(diff = idx*2)
  statisticsDF = rbind(statisticsDF, res)
}
```

```{r}
save(resL, statisticsDF, file = file_name)
```