Advances in microfluidics and sequencing technologies have allowed us to monitor biological systems at single-cell resolution. This comprehensive decomposition of complex tissues holds enormous potential in developmental biology and clinical research. However, the ever-increasing number of cells, technical noise, and high dropout rate pose significant computational challenges in scRNA-seq analysis.
To address this problem, we introduce single-cell Decomposition using Hierarchical Autoencoder (scDHA), which can efficiently detach noise from informative biological signals. The scDHA pipeline consists of two core modules. The first module is a non-negative kernel autoencoder that provides a non-negative, part-based representation of the data. Based on the weight distribution of the encoder, scDHA removes genes or components that have insignificant contributions to the representation. The second module is a Stacked Bayesian Self-learning Network that is built upon the Variational Autoencoder to project the data onto a low dimensional space. Using this informative and compact representation, many analyses can be performed with high accuracy and tractable time complexity (mostly linear or lower complexity).
In one joint framework, the scDHA software package conducts cell segregation through unsupervised learning, dimension reduction and visualization, cell classification, and time-trajectory inference.
scDHA can be installed fron GitHub or CRAN using below instruction. scDHA depends on the
torch package to build and train the autoencoders. When scDHA package is loaded, it will check for the availability of
torch package can be used to install
C++ libtorch, which is necessary for neural network computation.
#Install devtools: utils::install.packages('devtools') #Install the package from GitHub: devtools::install_github('duct317/scDHA') #With manual and vignette: devtools::install_github('duct317/scDHA', build_manual = T, build_vignettes = T) #Or from CRAN: install.packages("scDHA") #When the package is loaded, it will check for C++ libtorch library(scDHA) #libtorch can be installed using: torch::install_torch()
if (!requireNamespace("mclust", quietly = TRUE)) install.packages("mclust")
library(scDHA) #Load example data (Goolam dataset) data("Goolam") #Get data matrix and label data <- t(Goolam$data); label <- as.character(Goolam$label) #Log transform the data data <- log2(data + 1)
#Generate clustering result, the input matrix has rows as samples and columns as genes result <- scDHA(data, seed = 1) #The clustering result can be found here cluster <- result$cluster #Calculate adjusted Rand Index using mclust package ari <- round(mclust::adjustedRandIndex(cluster,label), 2) print(paste0("ARI = ", ari))
#Generate 2D representation, the input is the output from scDHA function result <- scDHA.vis(result, seed = 1) #Plot the representation of the dataset, different colors represent different cell types plot(result$pred, col=factor(label), xlab = "scDHA1", ylab = "scDHA2")
#Cell stage order in Goolam dataset cell.stages <- c("2cell", "4cell", "8cell", "16cell", "blast") #Generate pseudo-time for each cell, the input is the output from scDHA function result <- scDHA.pt(result, start.point = 1, seed = 1) #Calculate R-squared value representing correlation between inferred pseudo-time and cell stage order r2 <- round(cor(result$pt, as.numeric(factor(label, levels = cell.stages)))^2, digits = 2) #Plot pseudo-temporal ordering of cells in Goolam dataset plot(result$pt, factor(label, levels = cell.stages), xlab= "Pseudo Time", ylab = "Cell Stages", xaxt="n", yaxt="n") axis(2, at=1:5,labels=cell.stages, las=2) text(x = 1, y = 4.5, labels = paste0("R2 = ", r2))
#Split data into training and testing sets set.seed(1) idx <- sample.int(nrow(data), size = round(nrow(data)*0.75)) train.x <- data[idx, ]; train.y <- label[idx] test.x <- data[-idx, ]; test.y <- label[-idx] #Predict the labels of cells in testing set, the input matrices have rows as samples and columns as genes prediction <- scDHA.class(train = train.x, train.label = train.y, test = test.x, seed = 1) #Calculate accuracy of the predictions acc <- round(sum(test.y == prediction)/length(test.y), 2) print(paste0("Accuracy = ", acc))