getMarkerGenes function

Hi!

I have noticed that there is no function to get the marker genes from a loom file. Something like this works, but may need some check ups (e.g. check that the clustering is actually included in the loom etc):

# This function is not yet included in SCopeLoomR
getMarkersPerCluster <- function(loom, clustering, log_FC_thr=1.5, adj_pval=0.05){
  # Get clustering index
  available_clusterings <- sapply(1:length(get_global_meta_data(loom)$clusterings), function (i) get_global_meta_data(loom)$clusterings[[i]]$name) 
  clustering_idx <- which(available_clusterings == clustering)
  # Get matrices
  ra <- loom[["row_attrs"]]
  genes <- ra[['Gene']][]
  avglog_mat <- ra[[paste0('ClusterMarkers_', clustering_idx, '_avg_logFC')]][]
  pval_mat <- ra[[paste0('ClusterMarkers_', clustering_idx, '_pval')]][]
  rownames(avglog_mat) <- genes
  rownames(pval_mat) <- genes
  # Get cluster names
  cluster_names <- sapply(1:ncol(avglog_mat), function (i) get_global_meta_data(loom)$clusterings[[clustering_idx]]$clusters[[i]]$description)
  colnames(avglog_mat) <- rep('avg_logFC', ncol(avglog_mat))
  colnames(pval_mat) <- rep('adj_pval', ncol(pval_mat))
  # Create list
  marker_list <- lapply(1:ncol(avglog_mat), function (i) cbind(avglog_mat[,i, drop=FALSE], pval_mat[,i, drop=FALSE]))
  marker_list <- lapply(marker_list, function(x) x <- x[which(x[,1] >  log_FC_thr),])
  marker_list <- lapply(marker_list, function(x) x <- x[which(x[,2] < adj_pval),])
  marker_list <- lapply(marker_list, function(x) x <- x[with(x, order(-x[,'avg_logFC'], x[,'adj_pval'])),])
  names(marker_list) <- cluster_names
  return(marker_list)
} 

markers_list <- getMarkersPerCluster(loom, "Leiden resolution 0.8", log_FC_thr=1.5, adj_pval=0.05)

aertslab / SCopeLoomR

getMarkerGenes function #24