Raw data generated by a sequencing machine are processed to have read counts or number of molecules. The count information can be represented in a matrix which shows the gene expression across cells (Fig. 1A). To exclude low-quality information from scRNAseq, a series of quality controls (QC) are generally required. QC is employed to ensure that the quality of scRNAseq data is sufficient for subsequent downstream analysis. Common practices include removing cells with a low count depth or few detected genes. Some transcriptomic information can be from multiple cells (doublets). These unwanted doublets can be removed using doublet detection tools including Scrublet19 and DoubletFinder.20

Once the matrix is obtained, normalization can be considered to compensate for cells with different numbers of barcodes or read depth. Normalization scales count data to a relative expression abundance. In addition, the batch effect should be removed when considering scRNAseq data from several sources. The batch effect is a common source of technical variation that can arise from various sources, including disparate cell dissociation protocol, library preparation, and sequencing platforms. Comprehensive benchmarking tests were applied for normalization and batch effect correction for the dataset for droplet-based scRNAseq data21 showing that pooling-based size factor estimation by Scran22 is one of the best working normalization approaches. The Seurat package is equipped with batch effect correction,17 but trying diverse approaches after considering the underlying biological processes is also recommended for batch effect correction.

Clustering is a useful approach to identify various types of cells from scRNAseq. Hierarchical clustering has been widely used for bulk-cell RNAseq analysis. Butcher et al.11 used hierarchical clustering for 270 cells and identified 3 groups of cells. As the number of cells increases, the need for new approaches to handle high-dimensional scRNAseq data arises. Various algorithms for processing multi-dimensional scRNAseq data have been developed. Generally, these clustering methods rely on dimension reduction algorithms to avoid unnecessary noise (features). Dimensionality reduction approaches, such as principal components analysis (PCA), have been widely used to analyze both bulk RNA-seq data as well as scRNAseq data.1723 For instance, PCA allows the conversion of a higher dimensional dataset into a lower, often 2 or 3, dimensional dataset with more important and uncorrelated variables (dimensions) called principal components. Subsequently, cells can be clustered and envisioned in 2- or 3-dimensional (2D or 3D) space. Similarly, t-distributed stochastic neighbor embedding (tSNE) has also been successfully used to visualize cells in a reduced space.17 Moreover, tSNE is used to visualize cells in 2D or 3D space while reflecting true distances in the original space as far as possible so that cells of a particular cell type tend to be located nearby in 2D space. Finally, clustering is performed on cells by grouping them in the reduced space. Fig. 1B shows the tSNE-based analysis using the dataset by Lin et al.16 to understand macrophage heterogeneity during atherosclerosis progression and repression. From the tSNE plot, we can identify sub-clusters with V-erbA-related protein 2 (Ear2) expression (corresponds to Rentnla^hiEar2^hi macrophage) and interferon (IFN) regulatory factor 7 (type 1 IFN signature) (Fig. 1B). In a similar manner, cell types are assigned based on a priori known marker genes. Cochain et al.12 identified cells, such as B cells (using Cd79a, Cd79b, Ly6d, and Mzb1), C-X-C chemokine receptor type (CXCR)6+ T cells (using CXCR6, Icos, Cd3g, and Il7R) and natural killer cells (using Klrb1c, Ncr1, Klra8, and Klrc1), based on the associated marker genes.

A number of algorithms have also been exclusively proposed for scRNAseq data analysis, including SC324 and SIMLR.25 CellBIC was designed to identify small cell subpopulations without losing information by dimension reduction.26 GiniClust has also been proposed to identify rare cell population8. Recent advances allow for ultra-fast clustering of more than 1 million cells.27

One of the downstream applications of scRNAseq analysis is the comparison of cell compositions. For instance, Cochain et al.12 compared the number of cells from the control and the diseased aortas for each of the 3 clusters of macrophages and found that TREM2+ macrophages were almost exclusively observed in the cells from diseased aortas. In addition, the same quantitation could also provide an estimation of cell composition of bulk-cell RNAseq. This approach may be particularly useful when samples are collected from a different section of tissue. If scRNAseq is provided for a section (so that cell subpopulations are obtained), the cell composition of another section can be estimated from the bulk RNAseq using computational deconvolution based on scRNAseq28 (Fig. 2). Winkel et al.13 used CIBERSORT29 to perform deconvolution of cells using bulk-RNA-seq from the media, adventitia, lesion and adventitia + ATLO.

When cells are represented in a lower dimensional space, those with similar transcriptomes will be located nearby on a plot, e.g. using tSNE. When cells are collected in different time stamps during differentiation, mature cells will be located far from progenitors, and cells being differentiated will be located in the middle. The path that links the cells can be regarded as a “pseudo” time9 (Fig. 3). This allows for longitudinal analysis of gene expression (e.g. development). Pseudo-time can be used to model transcriptomic changes during the development of atherosclerosis. Gene expressions can be analyzed along pseudo-time. For instance, the expression level of elastin deceases during direct cardiomyocyte conversion, while the expression level of troponin I1, slow skeletal type increases (Fig. 3). Furthermore, Lin et al.16 used the pseudo-time analysis along the fate-mapping during atherosclerosis progression and regression. This analysis found 53 genes significantly correlated with pseudo-time score, including CXCR4 and Ctsd. Monocle has been used for pseudo-time analysis.9 TSCAN combines clustering with pseudo-time analysis.30 Partition-based graph abstraction could be useful when complex trajectories are expected.31

Reverse engineering reconstructs gene regulatory networks from gene expression information.32 It usually requires a large amount of expression data. By providing transcriptomic information for each single cell, scRNAseq can be a good resource for reconstructing the regulatory networks. Pseudo-time has also been used to identify potential downstream target genes9 (Fig. 3). Software tools such as SCODE were developed to reconstruct gene regulatory networks from scRNAseq data.

Another major limitation of current transcriptomic analysis workflow is that once the cells are isolated from tissue for scRNAseq, the cell location and orientation information is lost. To restore approximate location information, tissues can be mechanically sampled from different spots. For instance, Winkel et al.,13 used the spatial information by comparing cells from whole-atherosclerotic aortas versus aortic leukocytes. Another strategy has also been introduced in which barcoding the native tissue location has been proposed.33 This approach dissects the histological section with a grid, and each spot is barcoded to provide position information. Currently, a grid contains 10–30 cells depending on the tissue. The RNAs along with the position information are sequenced together. Spatial transcriptomics is now part of 10× Genomics (https://spatialtranscriptomics.com/). Another approach called MERFISH uses fluorescence in situ hybridization to provide spatial information in order to map the identity and location of specific cell types.34 Statistical test can be applied to identify genes whose spatial patterns are significantly different.35

The expression of ligand-receptor pairs can be utilized to study cell communication.36 Specifically, Gu et al.14 investigated the ligand-receptor interaction between cell subtypes by considering the transcriptomic levels of ligands and their corresponding receptors. In cells from the adventitia, this computational prediction showed the importance of mesenchyme populations for maintaining adventitial homeostasis. For instance, cells expressing Cd34 and Cav1 interacted with Sell and Icam1 expressed by inflammatory macrophages, respectively, which may potentially modulate leukocyte influx to the adventitia. This ligand-receptor pair analysis predicted the manner in which resident mesenchyme cells interacted with and attracted immune cells in vivo.