单细胞数据高级分析之初步降维和聚类|Dimensionalityreduction|Clustering

# cluster cells and remove contaminating populations
cat('Doing initial clustering\n')
cl <- cluster.the.data.simple(cm, expr, 9, seed=42)
md$init.cluster <- cl$clustering
# detection rate per cluster for some marker genes
goi <- c('Igfbp7', 'Col4a1', 'Neurod2', 'Neurod6')
det.rates <- apply(cm[goi, ] > 0, 1, function(x) aggregate(x, by=list(cl$clustering), FUN=mean)$x)
contam.clusters <- sort(unique(cl$clustering))[apply(det.rates > 1/3, 1, any)]
use.cells <- !(cl$clustering %in% contam.clusters)
cat('Of the', ncol(cm), 'cells', sum(!use.cells), 'are determined to be part of a contaminating cell population.\n')
cm <- cm[, use.cells]
expr <- expr[, use.cells]
md <- md[use.cells, ]

# for clustering
# ev.red.th: relative eigenvalue cutoff of 2%
dim.red <- function(expr, max.dim, ev.red.th, plot.title=NA, do.scale.result=FALSE) {
  cat('Dimensionality reduction via diffusion maps using', nrow(expr), 'genes and', ncol(expr), 'cells\n')
  if (sum(is.na(expr)) > 0) {
    dmat <- 1 - cor(expr, use = 'pairwise.complete.obs')
  } else {
    # distance 0 <=> correlation 1
    dmat <- 1 - cor(expr)
  }
  
  max.dim <- min(max.dim, nrow(dmat)/2)
  dmap <- diffuse(dmat, neigen=max.dim, maxdim=max.dim)
  ev <- dmap$eigenvals
  # relative eigenvalue cutoff of 2%, something like PCA
  ev.red <- ev/sum(ev)
  evdim <- rev(which(ev.red > ev.red.th))[1]
  
  if (is.character(plot.title)) {
    # Eigenvalues, We observe a substantial eigenvalue drop-off after the initial components, demonstrating that the majority of the variance is captured in the first few dimensions.
    plot(ev, ylim=c(0, max(ev)), main = plot.title)
    abline(v=evdim + 0.5, col='blue')
  }
  
  evdim <- max(2, evdim, na.rm=TRUE)
  cat('Using', evdim, 'significant DM coordinates\n')
  
  colnames(dmap$X) <- paste0('DMC', 1:ncol(dmap$X))
  res <- dmap$X[, 1:evdim]
  if (do.scale.result) {
    res <- scale(dmap$X[, 1:evdim])
  } 
  return(res)
}

# jaccard similarity
# rows in 'mat' are cells
jacc.sim <- function(mat, k) {
  # generate a sparse nearest neighbor matrix
  nn.indices <- get.knn(mat, k)$nn.index
  j <- as.numeric(t(nn.indices))
  i <- ((1:length(j))-1) %/% k + 1
  nn.mat <- sparseMatrix(i=i, j=j, x=1)
  rm(nn.indices, i, j)
  # turn nn matrix into SNN matrix and then into Jaccard similarity
  snn <- nn.mat %*% t(nn.mat)
  snn.summary <- summary(snn)
  snn <- sparseMatrix(i=snn.summary$i, j=snn.summary$j, x=snn.summary$x/(2*k-snn.summary$x))
  rm(snn.summary)
  return(snn)
}


cluster.the.data.simple <- function(cm, expr, k, sel.g=NA, min.mean=0.001, 
                                    min.cells=3, z.th=1, ev.red.th=0.02, seed=NULL, 
                                    max.dim=50) {
  if (all(is.na(sel.g))) {
    # no genes specified, use most variable genes
    # filter min.cells and min.mean
    # cm only for filtering
    goi <- rownames(expr)[apply(cm[rownames(expr), ]>0, 1, sum) >= min.cells & apply(cm[rownames(expr), ], 1, mean) >= min.mean]
    # gene sum
    sspr <- apply(expr[goi, ]^2, 1, sum)
    # scale the expression of all genes, only select the gene above z.th
    # need to plot the hist
    sel.g <- goi[scale(sqrt(sspr)) > z.th]
  }
  cat(sprintf('Selected %d variable genes\n', length(sel.g)))
  sel.g <- intersect(sel.g, rownames(expr))
  cat(sprintf('%d of these are in expression matrix.\n', length(sel.g)))
  
  if (is.numeric(seed)) {
    set.seed(seed)
  }
  
  dm <- dim.red(expr[sel.g, ], max.dim, ev.red.th, do.scale.result = TRUE)
  
  sim.mat <- jacc.sim(dm, k)
  
  gr <- graph_from_adjacency_matrix(sim.mat, mode='undirected', weighted=TRUE, diag=FALSE)
  cl <- as.numeric(membership(cluster_louvain(gr)))
  
  results <- list()
  results$dm <- dm
  results$clustering <- cl
  results$sel.g <- sel.g
  results$sim.mat <- sim.mat
  results$gr <- gr
  cat('Clustering table\n')
  print(table(results$clustering))
  return(results)
}