marker heatmap for shreejoy

lazyShreejoy.R

devtools::install_github('oganm/neuroExpressoAnalysis')
library(neuroExpressoAnalysis)
devtools::install_github('oganm/ogbox')
library(ogbox)
devtools::install_github('oganm/markerGeneProfile')
library(markerGeneProfile)

library(scales)
library(gplots)
library(stringr)
library(magrittr)
library(dplyr)
library(viridis)
library(reshape2)
library(pheatmap)

orderMarkers = function(genesDir =  '/home/omancarci/brainGenesManuscript/analysis/01.SelectGenes/'){
    
    genes = neuroExpressoAnalysis::mouseMarkerGenesCombined
    
    genes %<>% lapply(function(x){
        x = x[!grepl(pattern = '(?!^Pyramidal$)Pyra',x = names(x),perl = TRUE)]
        return(x)
    })
    
    
    genes %<>% lapply(function(x){x[!grepl(pattern='Microglia_',names(x))]})
    geneNames = names(genes)
    
    genesMicroarray = lapply(1:len(genes), function(i){
        nm = names(genes[[i]])
        nm = nm[nm %in% (n_expressoSamples %>% select(CellTypes,PyramidalDeep) %>% unlist %>% unique)]
        lapply(1:len(nm), function(j){
            print(paste(i,j))
            component1 = names(genes[i])
            component2 = names(genes[[i]][nm[j]])
            component1 = paste0(component1,'_')
            if (component1=='All_'){
                component1=''
            }
            if(component2 == 'Pyramidal'){
                scores = read.table(paste0(genesDir,'/Quick/', component1,'CellTypes/',component2))
            } else{
                scores = read.table(paste0(genesDir,'/Quick/', component1,'PyramidalDeep/',component2))
            }
            (scores %>% filter(V1 %in% genes[[i]][[nm[j]]]))[,]$V1 %>% unlist %>% as.char
        }) ->out
        names(out) = nm
        return(out)
    })
    
    names(genesMicroarray) = geneNames
    genesMicroarray = lapply(genesMicroarray, function(x){
        lapply(x,trimNAs)
    })
    
    genesSingleCell  = lapply(1:len(genes$Cortex), function(i){
        print(names(genes$Cortex[i]))
        if(names(genes$Cortex[i]) == 'Pyramidal'){
            scores = read.table(paste0('analysis/01.SelectGenes/QuickJustSingleCell//CellTypes/',names(genes$Cortex[i])))
        } else{
            scores = read.table(paste0('analysis/01.SelectGenes/QuickJustSingleCell//PyramidalDeep/',names(genes$Cortex[i])))
        }
        assertthat::are_equal(all(genes$Cortex[[i]] %in% scores$V1),TRUE)
        (scores %>% filter(V1 %in% genes$Cortex[[i]]))[,]$V1 %>% unlist %>% as.char
    })
    names(genesSingleCell) = names(genes$Cortex)
    genesSingleCell %<>% lapply(trimNAs)
    
    genesMicroarray %<>% lapply(function(x){
        names(x) = publishableNameDictionary$ShinyNames[match(names(x),publishableNameDictionary$PyramidalDeep)]
        return(x)
    })
    names(genesSingleCell) = publishableNameDictionary$ShinyNames[match(names(genesSingleCell),publishableNameDictionary$PyramidalDeep)]
    
    list(microarray = genesMicroarray,singleCell = genesSingleCell)
    
}



lazyShreejoy = function(expression, # expression data with rownames and gene names
                        n = 5, # plot top 5
                        outDir = 'output', # output directory
                        regions = NULL, # which regions to make
                        geneTransform =NULL , # transform function. eg. function(x){homologene::mouse2human(x)$humanGene}
                        genesDir = '/home/omancarci/brainGenesManuscript/analysis/01.SelectGenes/', # needs Quick and QuickJustSingleCell directories
                        annotation_col = NULL, # for annotating samples, see pheatmap documentation. a data frame with rownames as your colnames and columns as the things you want to annotate with
                        annotation_colors = list(), # for annotating samples, see pheatmap documentation. a named list of named vectors. names corresponds to column names of annotation_col)
                        gaps_col = NULL){ # for adding column gaps. see pheatmap documentation. vector of integers.
                        
    dir.create(outDir,showWarnings = FALSE)
    
    
    order = cellOrder %>% translatePublishable()
    
    orderedMarkers = orderMarkers(genesDir)
    
    if(!is.null(geneTransform)){
        orderedMarkers$microarray %<>% lapply(function(x){x %>% lapply(geneTransform)})
        orderedMarkers$singleCell %<>% lapply(geneTransform)
    }
    
    orderedMarkers$microarray %<>% lapply(function(x){x %>% lapply(function(y){y[y %in% rownames(expression)][1:n] %>% trimNAs()})})
    orderedMarkers$singleCell %<>% lapply(function(x){x[x %in% rownames(expression)][1:n] %>% trimNAs()})
    
    
    regionsToDo = names(orderedMarkers$microarray)
    if(!is.null(regions)){
        regionsToDo = regions
    }
    
    null = list(NULL)
    names(null) = NA
    
    for(x in regionsToDo){
        if(x == 'Cortex'){
            genes = names(orderedMarkers$microarray$Cortex) %>% sapply(function(x){
                c( orderedMarkers$microarray$Cortex[[x]],orderedMarkers$singleCell[[x]]) %>% unique
            })
            
            genes = genes[order] %>% trimElement(null)
            
        }  else{
            genes = orderedMarkers$microarray[[x]] %>% {.[order]} %>% trimElement(null)
        }
        
        relExp = expression[genes %>% unlist,]
        
        geneCellTypes = repIndiv(names(genes), genes %>% sapply(len))
        
        geneCols = toColor(geneCellTypes, cellColors())
        anotRow = data.frame('Specific Genes' = geneCellTypes, check.names=F)
        rownames(anotRow) = rownames(relExp)
        
        png(paste0(outDir,'/',x,'.png'),height=1900,width=2000)
        pheatmap(relExp,
                 gaps_row= anotRow$`Specific Genes` %>% duplicated %>% not %>% which %>% {.-1},
                 gaps_col = gaps_col,
                 fontsize = 30,
                 color=viridis(20),
                 cluster_rows=F,
                 cluster_cols=F,
                 annotation_row = anotRow,
                 annotation_col = annotation_col,
                 annotation_colors = c(list('Specific Genes' = geneCols$palette[order][order %in% names(genes)]
                 ),annotation_colors),
                 show_rownames=TRUE ,
                 show_colnames=FALSE
        )
        dev.off()
    }
}