R/findGenes.R

In zhupinglab/Aodav: Advanced omics data analysis and visualization

#' @title markerSelect
#' @description select most cell-type representative genes from differential gene expression matrix
#' @param DEGs differentially expressed gene matrix obtained from pwDEG, conservedMarkers, or
#' Seurat function 'FindAllMarkers'.
#' @param obj Seurat object
#' @param adjP cutoff of adjusted p value to filter genes
#' @param fc fold change to filter genes
#' @param only.pos selection of cluster highly expressed genes.
#' @param gene.num number of genes to select for each cluster
#' @param conserved indication of input matix is obtained from function 'conservedMarkers'
#' @param order sort genes by clusteirng their expression profiles
#' @examples
#' \dontrun{
#' markerSelect(DEGs, obj)
#' markerSelect(DEGs, obj, conserved = T) # this is for output from function 'conservedMarkers'
#' }
#' @export
#' @import Seurat
markerSelect <- function(
  DEGs,
  obj,
  adjP = 0.05,
  fc = log(2),
  only.pos = T,
  gene.num = NULL,
  conserved = F,
  order = F
){
  if(isTRUE(conserved)){
    p_val_adj <- 5
    avg_logFC <- 2
  }	else {
    p_val_adj <- "p_val_adj"
    avg_logFC <- "avg_logFC"
  }
  topMarkers <- unique(unlist(lapply(unique(DEGs$cluster), function(x){
    this.data <- DEGs[DEGs$cluster == x & DEGs[,p_val_adj] <= adjP & abs(DEGs[,avg_logFC]) >= fc,]
    if (isTRUE(only.pos)){
      this.data <- DEGs[DEGs$cluster == x & DEGs[,p_val_adj] <= adjP & DEGs[,avg_logFC] >= fc,]
    }
    this.data <- this.data[ order(this.data[,avg_logFC], decreasing=T),]
    if(is.null(gene.num)){
      return(this.data$gene)
    }
    return(head(this.data$gene, n=gene.num))
  })))
  if(isTRUE(order)){
    topMarkers <- topMarkers[hclust(as.dist(1-cor(t(obj@scale.data[topMarkers,]))), method="ward.D2")$order]
  }
  return(topMarkers)
}

#' @title pcaGenes
#' @description use pca method to extract genes with most contributions in first three dimensions. The pca genes could be used as gene list for downstream dimention reduction.
#' @param objs A list of Seurat objects
#' @param dim.use PCA dimensions used to select genes. default 1:3.
#' @param gene.num Number of genes selected on each PC.
#' @examples
#' \dontrun{
#' pcaGenes(obj)
#' }
#' @export
#' @import Seurat
pcaGenes <- function(objs, dim.use=1:3, gene.num=100){
  unique(unlist(lapply(objs, function(obj){
    this.pca <- prcomp(t(as.matrix(obj@data)))
    pc.dim <- dim.use
    pc.gene.num <- gene.num
    this.pca.genes <- unique(unlist(lapply(pc.dim, function(x){
      head(row.names(this.pca$rotation[order(abs(this.pca$rotation[,x]), decreasing=T),]), n=pc.gene.num)
    })))
    return(this.pca.genes)
  })))
}

#' @title pwDEG
#' @description find differentially expressed genes thoroughly for any pairs of clusters.
#' @param obj Seurat object
#' @examples
#' \dontrun{
#' pwDEG(obj)
#' }
#' @export
#' @import Seurat
pwDEG <- function(obj){
  DEGs <- c()
  for(i in 1:length(levels(obj@ident))){
    for(j in (i + 1):length(levels(obj@ident))){
      if (j > i && j <= length(levels(obj@ident))){
        print(paste0(levels(obj@ident)[i], " vs ", levels(obj@ident)[j]))
        this.marker <- FindMarkers(obj, ident.1=levels(obj@ident)[i], ident.2=levels(obj@ident)[j])
        this.marker$gene <- row.names(this.marker)
        this.marker$cluster <- paste0(levels(obj@ident)[i], "_",levels(obj@ident)[j])
        DEGs <- rbind(DEGs, this.marker)
      }
    }
  }
  return(DEGs)
}

#' @title conservedMarkers
#' @description find conserved markers for each cluster
#' @param obj Seurat object
#' @param only.pos selection of cluster highly expressed genes.
#' @param groupBy group cells by a specified feature.
#' @examples
#' \dontrun{
#' conservedMarkers(obj, groupBy)
#' }
#' @export
#' @import Seurat
conservedMarkers <- function(obj, groupBy, only.pos=T){
  conserved.markers <- c()
  for(i in levels(obj@ident)){
    print(paste0("Processing cluster c", i, " ..."))
    this.marker <- FindConservedMarkers(
      obj,
      ident.1 = i,
      grouping.var = groupBy,
      print.bar = FALSE, only.pos = only.pos
    )
    this.marker <- this.marker[order(this.marker[,2] + this.marker[,7], decreasing =T), ]
    this.marker$gene <- row.names(this.marker)
    this.marker$cluster <- i
    conserved.markers <- rbind(conserved.markers, this.marker)
  }
  return(conserved.markers)
}