load("data/cancer_cell_line/sincell_with_class_5cl.RData")
require(Seurat)
require(ggplot2)
require(tidyr)
require(dplyr)
require(randomForest)
require(GGally)
require(e1071)
require(ggh4x)
require(gtable)
require(reshape2)
require(patchwork)
source("script/20241231.R")

################################# Seurat #################################
obj <- CreateSeuratObject(counts = sce_sc_10x_5cl_qc@assays$data$counts, project = "celline",
                          meta.data = data.frame(sce_sc_10x_5cl_qc@colData))
obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj)
obj <- ScaleData(obj)
obj <- RunPCA(obj)
obj <- RunUMAP(obj, dims = 1:10)
Idents(obj) <- obj$cell_line
genes <- VariableFeatures(obj)

# PBMC
# obj <- readRDS("result/obj.pbmc3k.expression.rds")
# genes <- VariableFeatures(obj)

# Build the null data matrix
perm.data <- obj@assays$RNA$data[genes, ] %>%
  apply(.,1,sample) %>% t
rownames(perm.data) <- paste("NULL", rownames(perm.data), sep = "-")
new.obj <- CreateSeuratObject(counts = rbind(obj@assays$RNA$data[genes,], perm.data), project = "celline",
                              meta.data = data.frame(sce_sc_10x_5cl_qc@colData))
LayerData(new.obj, "data") <- LayerData(new.obj, "counts")
null.genes <- rownames(perm.data)
expand.genes <- rownames(new.obj)

################################# NNet #################################
new.obj <- prepare.seurat(new.obj, genes = expand.genes)
new.obj <- prepare.graph(new.obj)
new.obj <- prepare.reg(new.obj, predictors = expand.genes,
                       responses = expand.genes)

################################# Regression #################################
response.list <- tapply(genes, cut(1:length(genes), breaks = 200), list)

# Calculate power and FDR for one set of responses
one.iter <- function(responses){
  predictors <- paste("NULL",responses,sep="-")
  new.obj <- run.nn.reg(new.obj, responses, return.p.val = T, predictors = predictors)

  # Calculate power and FDR for one response
  evaluation <- function(p){
    null.p <- paste("NULL", p, sep = "-")
    p.val <- new.obj@misc$mod$p.val[p,p,]
    null.p.val <- new.obj@misc$mod$p.val[p,null.p,]
    expressed <- new.obj@assays$RNA$data[p,] != 0
    plot.df <- data.frame(expressed = expressed, p.val = p.val, null.p.val = null.p.val)
    plot.df %>% group_by(expressed) %>%
      dplyr::summarise(power = mean(p.val > 0.95), fdr =  mean(null.p.val > 0.95)) %>%
      dplyr::select(-1)%>%
      unlist(use.names = T)
  }

  sapply(responses, evaluation)
}

# Test
#one.iter(response.list[[3]])

result <- lapply(response.list, one.iter)


# Save results
# saveRDS(result, "result/benchmark3.result.rds")
# result <- readRDS("result/benchmark2.result.rds")

################################# Plot permutation #################################
# Some results are not in the correct format
result <- readRDS("result/benchmark3.result.rds")
bad.result <- which(sapply(result, is.list))
result[bad.result] <- lapply(result[bad.result], function(x){
  sapply(x, function(x){
    if(length(x) < 4){
      new.x <- rep(0,4)
      names(new.x) <- c("power1", "power2", "fdr1", "fdr2")
      new.x[c(2,4)] <- x
      new.x
    }else{
      x
    }
  })
})
result <- do.call(cbind, result)
result <- t(result) %>% data.frame() %>% mutate(F1 = 2*(power2 * (1-fdr2))/(1+power2-fdr2))
result <- melt(t(result))
result <- result %>% mutate(expressed = Var1 %in% c("power2","fdr2"),
                            type = ifelse(Var1 %in% c("power1","power2"), "Sensitivity", "1-Specificity"))
# result$expressed <- result$Var1 %in% c("power2","fdr2")
# result$type <- result$Var1 %in% c("power1","power2")

p1 <- ggplot(result) + geom_violin(aes(Var1, value, fill = expressed)) + ylim(0,1) +
  geom_jitter(aes(Var1, value), size = 0.8)+
  facet_wrap(~type, scales = "free_x")+
  scale_fill_manual("Response expressed?", values = c("darkblue","darkred")) +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.ticks = element_blank(), legend.position = "none") + xlab("Response expressed?") + ylab("Rate") +
  scale_x_discrete(labels = c("Not expressed", "Expressed", "Not expressed", "Expressed"))


result <- subset(result, Var1 == "F1")
mean.f1 <- round(mean(result$value),2)
p2 <- ggplot() + geom_violin(aes(y= "F1", x = result$value), fill = "darkred")+
  geom_jitter(aes(y= "F1", x = result$value),size = 0.8) + xlim(0,1) +
  xlab("F1 score") + ylab("")+
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.ticks = element_blank(), legend.position = "none") +
  ggtitle("Permutation based pruning", paste("F1 mean = ",mean.f1 ))
p1/p2 +
  plot_layout(heights = c(3, 1))

################################# Plot houskeeping #################################
# Some results are not in the correct format
result <- readRDS("result/benchmark3.result.rds")
bad.result <- which(sapply(result, is.list))
result[bad.result] <- lapply(result[bad.result], function(x){
  sapply(x, function(x){
    if(length(x) < 4){
      new.x <- rep(0,4)
      names(new.x) <- c("power1", "power2", "fdr1", "fdr2")
      new.x[c(2,4)] <- x
      new.x
    }else{
      x
    }
  })
})
result <- do.call(cbind, result)
result <- t(result) %>% data.frame() %>% mutate(F1 = 2*(power2 * (1-fdr2))/(1+power2-fdr2))
result <- melt(t(result))
result <- result %>% mutate(expressed = Var1 %in% c("power2","fdr2"),
                            type = ifelse(Var1 %in% c("power1","power2"), "Sensitivity", "1-Specificity"))
# result$expressed <- result$Var1 %in% c("power2","fdr2")
# result$type <- result$Var1 %in% c("power1","power2")

p1 <- ggplot(result) + geom_violin(aes(Var1, value, fill = expressed)) + ylim(0,1) +
  geom_jitter(aes(Var1, value), size = 0.8)+
  facet_wrap(~type, scales = "free_x")+
  scale_fill_manual("Response expressed?", values = c("darkblue","darkred")) +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.ticks = element_blank(), legend.position = "none") + xlab("Response expressed?") + ylab("Rate") +
  scale_x_discrete(labels = c("Not expressed", "Expressed", "Not expressed", "Expressed"))


result <- subset(result, Var1 == "F1")
mean.f1 <- round(mean(result$value),2)
p2 <- ggplot() + geom_violin(aes(y= "F1", x = result$value), fill = "darkred")+
  geom_jitter(aes(y= "F1", x = result$value),size = 0.8) + xlim(0,1) +
  xlab("F1 score") + ylab("")+
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.ticks = element_blank(), legend.position = "none") +
  ggtitle("Housekeeping gene based pruning", paste("F1 mean = ",mean.f1 ))
p1/p2 +
  plot_layout(heights = c(3, 1))

################################# Co-expression plot #################################

source("script/20241231.R")

p <- "SEC61G"
null.p <- paste("NULL", p, sep = "-")
new.obj <- run.nn.reg(new.obj, p, return.p.val = T)
mu <- new.obj@misc$mod$mus[p]
sigma <- new.obj@misc$mod$sigmas[p]

p1 <- FeaturePlot(obj, p, pt.size = 1, order = T)+xlab("UMAP1")+ylab("UMAP2")+
  theme_classic()+
  scale_color_gradientn("Scaled expression",colors =gg.color.spec(11))+
  theme(text = element_text(face = "bold", size = 13),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(),
        plot.title = element_text(hjust = 0.5, face = "bold.italic")
  )+ xlab("")+ylab("") +
  guides(color = guide_colorbar(theme = theme(
    legend.text = element_text(size = 0)
  )))

source("script/20240910.R")

new.obj1 <- run.nn.reg(new.obj, p, return.p.val = T)
mu1 <- new.obj1@misc$mod$mus[p]
sigma1 <- new.obj1@misc$mod$sigmas[p]

p2 <- ggplot() + geom_jitter(aes(x= "Permuted response", y= log(abs(new.obj@misc$mod$effect[p,null.p,]))), size= 0.8) +
  geom_jitter(aes(x= "Response", y= log(abs(new.obj@misc$mod$effect[p,p,]))), size= 0.8) +
  geom_errorbar(aes(x = "Null distribution", ymin = mu-1.64*sigma, ymax =mu + 1.64*sigma, color = "Permutation"),size = 2)+
  geom_errorbar(aes(x = "Null distribution", ymin = mu1-1.64*sigma1, ymax =mu1 + 1.64*sigma1, color = "Housekeeping"),size = 2)+
  theme_classic()+
  scale_color_manual("Error bar",values = c("grey40","grey80"))+
  theme(text = element_text(size = 13, face = "bold"),
        axis.text.x = element_text(angle = 15, hjust = 1),
        legend.position = "none") + xlab("Co-expression type") +
  ylab("Log absolute co-expression")
p1 + p2 + plot_layout(ncol = 2)


################################# Test 1 (not shown) #################################
result <- readRDS("result/benchmark2.result.rds")
bad.result <- which(sapply(result, is.list))
result[bad.result] <- lapply(result[bad.result], function(x){
  sapply(x, function(x){
    if(length(x) < 4){
      new.x <- rep(0,4)
      names(new.x) <- c("power1", "power2", "fdr1", "fdr2")
      new.x[c(2,4)] <- x
      new.x
    }else{
      x
    }
  })
})
result <- do.call(cbind, result)
result <- melt(result)
result <- result %>% mutate(expressed = Var1 %in% c("power2","fdr2"),
                            type = ifelse(Var1 %in% c("power1","power2"), "Sensitivity", "1-Specificity"))


result[result$Var1 == "power2",] %>% arrange(value) %>% head
result[result$Var1 == "fdr2",] %>% arrange(desc(value)) %>% head
result %>% subset(Var2 == "MAGEA10")
one.iter("SEC61G")
