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)
source("script/20240910.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

################################# NNet #################################
obj <- prepare.seurat(obj)
obj <- prepare.graph(obj)
obj <- prepare.reg(obj, predictors = VariableFeatures(obj),
                   responses = VariableFeatures(obj))
pcs <- obj@misc$setting$pcs
loadings <- obj@misc$setting$loadings
ids <- ncol(obj)
genes <- VariableFeatures(obj)
data <- obj@assays$RNA$scale.data[genes,]

################################# Regression function #################################
calc.der <- function(mod, pred.fun, npcs){
  d <- seq(-0.3,0.3, length.out = 15)
  eps <- min(d[d>0])
  x.new <- matrix(0, ncol = npcs, nrow = length(d))
  f.prime <- c()
  for(i in 1:ncol(x.new)){
    x.perm <- x.new
    x.perm[,i] <-  x.perm[,i] + d
    yhat.perm <- pred.fun(x.perm, mod)
    w.plus <-  exp(-1*abs(d-eps)/eps)
    p1 <- sum(w.plus*yhat.perm)/sum(w.plus)
    w.minus <-  exp(-1*abs(d+eps)/eps)
    p2 <- sum(w.minus*yhat.perm)/sum(w.minus)
    f.prime[i] <- (p1-p2)/(2*eps)
  }
  f.prime
}

# Ridge
ridge.reg <- function(y, x, lambda = 5){
  qr.mod <- rbind(x,diag(sqrt(lambda), ncol(x))) %>% qr
  v <- qr.qty(qr.mod, c(y, rep(0, ncol(x))))
  b <- backsolve(qr.R(qr.mod), v)
  list(par.derv = b, mod = b)
}
ridge.pred.fun <- function(x, mod) x %*% mod

# Random Forest
rf.reg <- function(y, x){
  rf.mod <- randomForest(x, as.vector(y), ntree = 100)
  b <- calc.der(rf.mod, rf.pred.fun, ncol(x))
  list(par.derv = b, mod = rf.mod)
}
rf.pred.fun <- function(x, mod) predict(mod, newdata = x)

# SVM
svm.reg <- function(y, x){
  svm.mod <- svm(y ~ ., data = x, scale = FALSE, kernel = "radial")
  b <- calc.der(svm.mod, svm.pred.fun, ncol(x))
  list(par.derv = b, mod = svm.mod)
}
svm.pred.fun <- function(x, mod) predict(mod, x)

################################# Run experiment #################################
experiment <- function(){
  p <- sample(genes, size = 1)
  q <- sample(genes, size = 1)
  n <- sample(ids, size = 1)

  y <- obj@misc$setting$lra[,p]
  neighbours <- obj@misc$setting$nn.idx[n,]

  y <- scale(y[neighbours])
  x <- scale(pcs[neighbours,])
  center.x <- attr(x, "scaled:center")
  scale.x <- attr(x, "scaled:scale")
  var.q <- var(data[q,neighbours])

  perm.x <- data[genes,neighbours]
  perm.x[q,] <- sample(perm.x[q,])
  perm.x <- crossprod(perm.x, loadings)
  perm.x <- scale(perm.x, center = center.x, scale = scale.x)

  calc.imp <- function(reg.fun, pred.fun, suffix){
    reg.out <- reg.fun(y,x)
    par.derv <- reg.out$par.derv
    imp <- 2 * var.q * ((loadings[q,]/scale.x) %*% par.derv)^2
    perm.imp <- mean((pred.fun(x, reg.out$mod)-pred.fun(perm.x, reg.out$mod))^2)
    df <- data.frame(sqrt(imp), sqrt(perm.imp))
    colnames(df) <- paste(suffix, c("sqrt.imp", "sqrt.perm.imp"), sep = ".")
    df
  }

  ridge.imp <- calc.imp(ridge.reg, ridge.pred.fun, "ridge")
  rf.imp <- calc.imp(rf.reg, rf.pred.fun, "rf")
  svm.imp <- calc.imp(svm.reg, svm.pred.fun, "svm")

  bind_cols(ridge.imp, rf.imp, svm.imp)
}
result <- sapply(1:1000, function(x) experiment(), simplify = F) %>% do.call(rbind,.)
colnames(result) <- c("Ridge approx", "Ridge perm", "RF approx", "RF perm", "SVM approx", "SVM perm")

################################# Make biplot #################################
cor_func <- function(data, mapping, method, symbol, ...){
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)

  corr <- cor(x, y, method=method, use='complete.obs')

  colFn <- colorRampPalette(c("brown1", "white", "dodgerblue"),
                            interpolate ='spline')
  fill <- rev(colFn(100))[findInterval(corr, seq(-1, 1, length = 100))]

  ggally_text(
    label = paste(symbol, as.character(round(corr, 2))),
    mapping = aes(),
    xP = 0.5, yP = 0.5,
    color = 'black',
    ...) +
    theme_void() +
    theme(panel.background = element_rect(fill = fill))
}

p <- ggpairs(result, diag  = NULL,
             upper = list(continuous = wrap(cor_func,
                                            method = 'spearman', symbol = expression('\u03C1 ='), size = 5)),
             lower = list(continuous = function(data, mapping, ...) {
               ggally_points(data = data, mapping = mapping, size = 0.5) + theme_bw() + geom_abline(linetype = "dashed")
             }
             )
)+
  theme(text = element_text(size = 13, face = "bold"), axis.text.x = element_text(angle = 90),
        strip.text = element_text(size = 13))


p11 <- ggally_text("Ridge approx", aes(size = 3, fontface = "bold"), color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkred"))
p22 <- ggally_text("Ridge perm", aes(size = 3, fontface = "bold"), color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkred"))
p[1,1] <- p11
p[2,2] <- p22

p33 <- ggally_text("RF approx", aes(size = 3, fontface = "bold") , color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkblue"))
p44 <- ggally_text("RF perm", aes(size = 3 , fontface = "bold"), color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkblue"))
p[3,3] <- p33
p[4,4] <- p44

p55 <- ggally_text("SVM approx", aes(size = 3 , fontface = "bold"), color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkgreen"))
p66 <- ggally_text("SVM perm", aes(size = 3 , fontface = "bold"), color = "white") +
  theme_classic() + theme(panel.background = element_rect(fill = "darkgreen"))
p[5,5] <- p55
p[6,6] <- p66


p
