# 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 1dmat <- 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 PCAev.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 matrixnn.indices <- get.knn(mat, k)$nn.indexj <- as.numeric(t(nn.indices))i <- ((1:length(j))-1) %/% k + 1nn.mat <- sparseMatrix(i=i, j=j, x=1)rm(nn.indices, i, j)# turn nn matrix into SNN matrix and then into Jaccard similaritysnn <- 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 filteringgoi <- rownames(expr)[apply(cm[rownames(expr), ]>0, 1, sum) >= min.cells & apply(cm[rownames(expr), ], 1, mean) >= min.mean]# gene sumsspr <- apply(expr[goi, ]^2, 1, sum)# scale the expression of all genes, only select the gene above z.th# need to plot the histsel.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 <- dmresults$clustering <- clresults$sel.g <- sel.gresults$sim.mat <- sim.matresults$gr <- grcat('Clustering table\n')print(table(results$clustering))return(results)
}