explore_seurat.R

library(rlc)
library(matrixStats)
library(biomaRt)
library(Matrix)
library(umap)

#sobj <- readRDS("make_analysis_out/SN010_E115/SN010_E115_normalized.rds")
makeENSMEBLlink <- function(geneID){
  sprintf(
    "<a href='http://www.ensembl.org/Mus_musculus/Gene/Summary?db=core;g=%s' target='_blank'>%s</a>",
    geneID, geneID
  )
}

makeENSMEBLlinks <- function(geneIDs){
  out <- NULL
  for(geneID in geneIDs){
    out <- c(
      out,
      makeENSMEBLlink(geneID)
    )
  }

  return(out)
}

makeLinks <- function(geneIDs, linkFunction){
  out <- NULL
  for(geneID in geneIDs){
    out <- c(
      out,
      linkFunction(geneID)
    )
  }

  return(out)
}

makeAllenAdultLink <- function(geneSymbol){
  sprintf(
    "<a href='http://mouse.brain-map.org/search/show?page_num=0&page_size=27&no_paging=false&exact_match=false&search_term=%s&search_type=gene' target='_blank'>%s</a>",
    geneSymbol, geneSymbol
  )
}

makeAllenDevoLink <- function(geneSymbol){
  sprintf(
    "<a href='http://developingmouse.brain-map.org/search/show?page_num=0&page_size=4&no_paging=false&exact_match=false&search_term=%s&search_type=gene' target='_blank'>%s</a>",
    geneSymbol, geneSymbol
  )
}

reformatLinkVector <- function(link.mat){
  reformatted <- NULL
  for(row in 1:nrow(link.mat)){
    reformatted <- c(reformatted, link.mat[row,])
  }

  return(reformatted)
}

generateLinks <- function(df, species="mou"){
  genes <- df$external_gene_name
  geneIDs <- rownames(df)

  df$ENSEMBL <- "%s"
  df$Allen_Adult <- "%s"
  df$Allen_Devo <- "%s"


  df.html <- hwriter::hwrite(df)

  links <- reformatLinkVector(
    cbind(
      makeENSMEBLlinks(geneIDs),
      makeLinks(genes, makeAllenAdultLink),
      makeLinks(genes, makeAllenDevoLink)
    )
  )

  out.html <- do.call(sprintf, as.list(c(df.html, links)))


  return(out.html)
}

explore_seurat <- function(seurat_object_rds_path, species, singleSample = T){

  cat("Reading the data\n")
  sobj <<- readRDS(seurat_object_rds_path)

  if(tolower(species) == "mou"){
    ensembl <- useEnsembl("ensembl", host = "http://dec2017.archive.ensembl.org", dataset = "mmusculus_gene_ensembl")
  }else if(tolower(species) == "hum"){
    ensembl <- useEnsembl("ensembl", host = "http://dec2017.archive.ensembl.org", dataset = "hsapiens_gene_ensembl")
  }
  # MOU data


  # HUM data
  #

  if(singleSample){
    all.genes <- t([email protected])
    vars <- apply(all.genes, 2, var)
    sf <- rowSums(all.genes)
    means <- apply(all.genes / sf, 2, mean)
    vpm <- vars/means

    vpm.data <- data.frame(means = means, vpm = vpm)
    rownames(vpm.data) <- colnames(all.genes)
  }


  exprs <- sobj@data
  sobj.umap <- sobj@[email protected]
  hvg.genes <- scale(t(sobj@data[[email protected], ]))
  g.pca <- irlba::prcomp_irlba(t(hvg.genes), n = 25)
  rownames(g.pca$x) <- colnames(hvg.genes)
  hvg.genes.umap <- umap(g.pca$x, metric = "correlation", method = "umap-learn", min_dist = .05, n_neighbors = 10)



  selGene <- rownames(exprs)[1]
  selGenes <- colnames(hvg.genes)
  selCell <- rownames(hvg.genes)[1]
  meanMarkedExp <- rowMeans(exprs)

  cellColVec <- rowMeans(hvg.genes)

  # Getting the gene symbols from ENSEMBL


  geneInfo <- getBM(
    attributes = c(
      "ensembl_gene_id",
      "external_gene_name"
    ),
    filters = "ensembl_gene_id",
    values = rownames(exprs),
    mart = ensembl
  )



  #########################################
  ##### Starting the rlc session ##########
  #########################################
  openPage( useViewer = FALSE, layout = "table2x2" )

  lc_scatter(
    dat(
      x = sobj.umap[,1],
      y = sobj.umap[,2],
      #colourValue = log2(as.numeric(as.matrix(sn15.exprs[selGene,])) + 1),
      colourValue = log1p(cellColVec + 1),
      palette = RColorBrewer::brewer.pal( 9, "YlOrRd" ),
      size = 2,
      on_click = function(i){
        selCell <<- rownames(sobj.umap)[i]
        updateCharts("A2")
      },
      on_marked = function(){
        markedCells <<- getMarked("A1")
        if(!is.null(markedCells)){
          meanMarkedExp <<- rowMeans(exprs[,markedCells])
          updateCharts("A2")
        }else{
          meanMarkedExp <<- rowMeans(exprs)
        }

        showHighGenes( markedCells )
      }),
    "A1"
  )


  lc_scatter(
    dat(
      x = hvg.genes.umap$layout[,1],
      y = hvg.genes.umap$layout[,2],
      #colourValue = meanMarkedExp[rownames(genes15.umap$layout)],
      #palette = RColorBrewer::brewer.pal( 9, "YlOrRd" ),
      size = 2,
      on_click = function(i){
        selGene <<- rownames(hvg.genes.umap$layout)[i]
        if(!is.null(selGene)){
          cellColVec <<- hvg.genes[,selGene]
        }else{
          cellColVec <<- rowMeans(hvg.genes)
          # markGene <<- rep(0, nrow(sn15.exprs))
          # markGene[selGene == rownames(sn15.exprs)] <<- 1
        }
        updateCharts(c("A1", "C1", "C2", "C3"))
      },
      on_marked = function(){
        selGenes <<- getMarked("A2")
        if(!is.null(selGenes)){
          cellColVec <<- apply(hvg.genes[,selGenes],1,mean)
        }else{
          cellColVec <<- rowMeans(hvg.genes)
        }
        updateCharts("A1")
      }
    ),
    "A2"
  )

  if(singleSample){
    lc_scatter(
      dat(
        x = means,
        y = vpm,
        #colourValue = markGene,
        size = 2,
        logScaleX = 2,
        logScaleY = 2,
        on_click = function(i){
          genes.vec <- exprs[i,]
          cellColVec <<- genes.vec
          updateCharts("A1")
          selGene <<- rownames(exprs)[i]
          updateCharts(c("C1", "C2", "C3"))
        }
      ),
      "B1"
    )

  }


  showHighGenes <- function( marked ){

    # If no genes are marked, clear output, do nothing else
    if( is.null(markedCells) ) {
      #lc_html( dat(content= ""), "B2" )
      rlc::removeChart("B2")
      return()
    }

    rlc::removeChart("B2")
    exprsMarked <- as.matrix(exprs[ , marked ])
    exprsUnmarked <- as.matrix(exprs[ , -marked ])
    df <- data.frame(
      meanMarked = rowMeans( exprsMarked ),
      meanUnmarked = rowMeans( exprsUnmarked ),
      sdMarked = rowSds( exprsMarked ),
      sdUnmarked = rowSds( exprsUnmarked )
    )

    print(head(df))

    df$sepScore <- ( df$meanMarked - df$meanUnmarked ) / pmax( df$sdMarked + df$sdUnmarked, 0.0002 )

    df <- merge(df, geneInfo, by.x = 0, by.y = 1)

    rownames(df) <- df$Row.names
    df <- df[,-1]

    #html_table <- hwriter::hwrite( head( df[ order(-df$sepScore), ], 20 ) )

    # Now, send the html table to appear in slot B2
    output.df <<- df
    # lc_html( dat(content = head( df[ order(-df$sepScore), ], 20 )) , "B2" )

    lc_html(dat(content = generateLinks(head( df[ order(-df$sepScore), ], 20 ))), "B2")
  }


  lc_html(
    dat( content = sprintf(
      "<a href='http://www.ensembl.org/Mus_musculus/Gene/Summary?db=core;g=%s' target='_blank'>ENSEMBL entry for %s</a>",
      selGene, selGene) ),
    "C1"
  )


  lc_html(
    dat(
      content = sprintf(
        "<a href='http://mouse.brain-map.org/search/show?page_num=0&page_size=27&no_paging=false&exact_match=false&search_term=%s&search_type=gene' target='_blank'>ABA search for %s (adult)</a>",
        geneInfo[geneInfo$ensembl_gene_id == selGene,2], geneInfo[geneInfo$ensembl_gene_id == selGene,2]
      )
    ),
    "C2"
  )

  lc_html(
    dat(
      content = sprintf(
        "<a href='http://developingmouse.brain-map.org/search/show?page_num=0&page_size=4&no_paging=false&exact_match=false&search_term=%s&search_type=gene' target='_blank'>ABA search for %s (developing)</a>",
        geneInfo[geneInfo$ensembl_gene_id == selGene,2], geneInfo[geneInfo$ensembl_gene_id == selGene,2]
      )
    ),
    "C3"
  )
}

markGene <- function(exprs, gene){
  cellColVec <<- exprs[gene,]
  updateCharts("A1")
}

# http://mouse.brain-map.org/search/show?page_num=0&page_size=27&no_paging=false&exact_match=false&search_term=Vav3&search_type=gene