############################################################
############################################################
library(viridis)
tableau10medium = scater:::.get_palette("tableau10medium")
okabe <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
## #pal = colorRampPalette(okabe)(23)
## #barplotPalette(pal)
bupu = c("#f7fcfd", "#e0ecf4", "#bfd3e6", "#9ebcda", "#8c96c6", "#8c6bb1", "#88419d", "#810f7c", "#4d004b")
## pal = colorRampPalette(bupu)(100)

############################################################
redRows <- function(x, douniq = TRUE, sep = ","){
    if(douniq){
        xo = paste(sort(unique(x)), collapse = sep)
    } else {
        xo = paste(sort(x), collapse = sep)
    }
    return(xo)
}

## ordering stage, zone, name: works for most batches, for rest we set to "-" and add later
orderLevels <- function(calls.dt, allStages, allZones){
    if(any(grepl("_", calls.dt$name))){
        tmp = do.call(rbind, strsplit(calls.dt$name, split = "_"))
        stage = tmp[,1]
        zone = tmp[,2]
        stage.u = unique(stage)

        ## making sure name format is stage_zone
        stopifnot(all(grepl("^E|^P", stage.u[!is.na(stage.u)])))
        stopifnot(all(stage.u[!is.na(stage.u)] %in% allStages))
        zone.u = unique(zone)
        stopifnot(all(zone.u[!is.na(zone.u)] %in% allZones))
        
        calls.dt$stage = stage
        calls.dt$zone = zone
        
        orderByStage <- function(this_stages, what2order){
            oo = c()
            ll = levels(factor(what2order))
            for(i in 1:length(this_stages)){
                oo = c(oo, grep(this_stages[i], ll))
            }
            return(oo)
        }
        ## which stages present in data
        this_stages = allStages[allStages %in% stage]
        
        ## order name
        ll = levels(factor(calls.dt$name))
        oo = orderByStage(this_stages, calls.dt$name)
        calls.dt$name = factor(calls.dt$name, levels = ll[oo])
        
        ## order stage
        ll = levels(factor(calls.dt$stage))
        oo = orderByStage(this_stages, calls.dt$stage)
        calls.dt$stage = factor(calls.dt$stage, levels = ll[oo])
        
        ## order zone
        calls.dt$zone = factor(calls.dt$zone, levels = c("LA", "LV", "RA", "RV"))
        
    } else {
        stage = "-"
        zone = "-"
        calls.dt$stage = stage
        calls.dt$zone = zone
        ## order name
        ll = levels(factor(calls.dt$name))
        calls.dt$name = factor(calls.dt$name, levels = ll)
        
    }
    
    return(calls.dt)
}

## simple convert seu to sce
my.convertToSCE <- function(this_seu){
    this_sce = SingleCellExperiment(assays = list(counts = this_seu@assays$RNA@counts,
                                                  logcounts = this_seu@assays$RNA@data,
                                                  scale.data = this_seu@assays$RNA@scale.data
                                                  ),
                                    reducedDims = SimpleList(pca = this_seu@reductions$pca@cell.embeddings,
                                                             tsne = this_seu@reductions$tsne@cell.embeddings,
                                                             umap = this_seu@reductions$umap@cell.embeddings,
                                                             harmony = this_seu@reductions$harmony@cell.embeddings))
    colData(this_sce) = DataFrame(this_seu@meta.data)
    return(this_sce)
}

## gprofiler-related
run_gProf <- function(this.ensg){
    gostres <- gost(query = this.ensg,
                    organism = "mmusculus",
                    ordered_query = FALSE, 
                    multi_query = FALSE,
                    significant = FALSE,
                    exclude_iea = FALSE, 
                    measure_underrepresentation = FALSE,
                    evcodes = TRUE, 
                    user_threshold = 0.05,
                    correction_method = "g_SCS", 
                    domain_scope = "annotated",
                    custom_bg = NULL, 
                    numeric_ns = "",
                    sources = c("GO:BP"),
                    as_short_link = FALSE)
    return(gostres)
}

clean_gProf <- function(gostres){
    gores.dt = data.table(gostres$result)
    gores.dt[,"evidence_codes" := NULL]
    gores.dt[, "neglog10pval" := -log10(p_value)]
    ## limiting to terms that are not too generic
    ## gores.dt = gores.dt[term_size <= 400, c(3,4,6,11, 15)]
    gores1.dt = gores.dt[term_size <= 400, c(3,4,5,6,9,11, 16)]
    toret.dt = gores1.dt[, c("term_id", "term_name",
                             "p_value", "neglog10pval",
                             "term_size", "query_size",
                             "intersection_size")]
    return(toret.dt)
}


getMusOrtho <- function(this.hu.gg){
    human <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
    mouse <- useMart(biomart="ensembl", dataset="mmusculus_gene_ensembl")
    biomaRt:::martHost(human)
    biomaRt:::martHost(mouse)
    
    genesmus <- getLDS(attributes = "hgnc_symbol",
                       filters = "hgnc_symbol",
                       values = this.hu.gg,
                       mart = human,
                       attributesL = "mgi_symbol",
                       martL = mouse,
                       uniqueRows = TRUE)
    genesmus2 <- rename(genesmus,
                        "human" = "HGNC.symbol",
                        "mouse" = "MGI.symbol") %>%
        group_by(human) %>%
        mutate(nmouse = n()) %>%
        group_by(mouse) %>%
        mutate(nhuman = n()) %>%
        filter(nmouse == 1, nhuman == 1) %>%
        dplyr::select(-nmouse, -nhuman)
    gg = as.vector(genesmus2$mouse)
    ## needed to add this:
    ## https://github.com/grimbough/biomaRt/issues/31
    httr::set_config(httr::config(ssl_verifypeer = FALSE))
    ens <- getBM(attributes = c("ensembl_gene_id",
                                "mgi_symbol",
                                "chromosome_name"),
                 filter = "mgi_symbol",
                 values = gg,
                 mart = mouse,
                 useCache = FALSE) ## needed to add this
    ## https://support.bioconductor.org/p/p132704/
    genesmus3 <- inner_join(genesmus2, ens, by=c("mouse"="mgi_symbol")) %>%
        dplyr::filter(!grepl("PATCH", chromosome_name)) %>%
        dplyr::filter(!grepl("CHR_MM", chromosome_name)) %>%
        dplyr::select(-chromosome_name)
    genesmus3.dt = data.table(genesmus3)
    setnames(genesmus3.dt, "human", "Gene")
    length(unique(genesmus3.dt$Gene)) == length(genesmus3.dt$Gene) ## FALSE
    genesmus3 = data.frame(genesmus3.dt)
    
    genesmus3.dt[mouse %in% genesmus3.dt[which(duplicated(genesmus3.dt$mouse))]$mouse]
    ## since in all cases its the ensembl id that is different, but gene name is same in hu and mus
    genesmus3.dt = unique(genesmus3.dt[, c("Gene", "mouse")])
    setnames(genesmus3.dt, c("hu_gene", "mus_gene"))
    return(genesmus3.dt)
}


splitN <- function(N, nmax = 1000000){
    nfiles = ceiling(N/nmax)
    nlims = (0: (nfiles-1)) * nmax
    if(nfiles > 1){
        cat("Splitting file into #files = ", nfiles, " \n")
    } else {
        cat("File is small enough to not be split \n")
    }
    out.list = list()
    for(i in 1:nfiles){
        ## cat(i, nlims[i]+1, min(nlims[i+1],N, na.rm=T), "\n")
        ## rng = (nlims[i]+1):min(nlims[i+1],N, na.rm=T)
        this.rng = c((nlims[i]+1), min(nlims[i+1],N, na.rm=T))
        out.list[[i]] = this.rng
    }
    return(out.list)
}