scRICA: An R package for multiple-sample single-cell RNA-seq data integrative comparative analysis

Operation

scRICA categorizes the workflow of integrative and comparative analysis using multiple-sample scRNA-seq data into four steps ( Figure 1A): ‘Step 1’ for pre-processing and quality control; ‘Step 2’ for multi-sample integration; ‘Step 3’ for visualization by attribute groups; and ‘Step 4’ for downstream analysis including differential expression (DE) analysis, pseudo-time trajectory analysis, and cell clusters identification. Everything is implemented in R as an R package with various functions implementations, additionally it also offers command line execution option to make the implementations easier on High-Performance Computing (HPC) systems. Instructions for installing prerequisite packages are available at https://github.com/yan-cri/scRICA. Users can follow these instructions to set up the appropriate analysis environment in R.

Figure 1. scRICA schematic overview.

A. scRICA workflow analysis steps; B. scRICA analysis results outline; C&D. Cellular composition percentage stacking bar plots by different sample’s attribute groups; E&F. Marker genes dot plots with respect to different sample’s attribute groups; G&H. Heat map of gene of interest with respect to the selected sample’s attribute groups and cell types.

Workflow implementation

‘Step 1’ performs pre-processing and quality control checks for each individual sample listed in the input metadata table. This step is initiated by the function ‘processQC()’, with main parameters inputting from a user-provided metadata table (details in input and output structures). ‘Step 1’ performs: 1) counts importation into a SeuratObject; 2) doublet/multiplet detection and removal using DoubletDecon⁶; 3) quality control results visualization; 4) mitochondrial and ribosomal contents summarization; and 5) cell numbers summarization. It has additional options including: 1) ‘genomeSpecies’, which allows users to select the species of the reference genome; 2) ‘mtPerCutoff’, the cut-off for mitochondrial content values. Values above the cut-off indicate cells with low viability; 3) ‘extraFilter’, which allows users to further eliminate certain cells from the analysis. These cells need to be specified in a separate file with respect to each sample, and the full path of these files should be provided in the metadata table column ‘filterFname’; and 4) ‘multi-omics ’, which specifies the format of the input count matrices.

‘Step 2’ integrates all samples listed in the input metadata table using two popular integration methods: the CCA and RPCA algorithms in Seurat.^7,8 This step is performed on cells that pass the quality control from ‘Step 1’ via the function ‘getClusterMarkers()’. In ‘Step 2’, the number of gene features used for integration can be specified by the option ‘nfeatures’. Additionally, users can exclude ribosomal genes from the gene set by setting ‘ribo.removal = TRUE’. By default, the identified cell clusters are listed as numbers, (i.e., 0, 1, 2, etc.). scRICA allows users to further annotate cell clusters via the option ‘newAnnotationRscriptName’ for downstream analysis.

‘Step 3’ provides various visualization techniques to inspect cellular compositions and gene expression patterns of different cell types, with respect to sample attributes specified in the option ‘expCondCheck’. For example, we can analyze fallopian tube data from the Human Cell Atlas by anatomic sites (isthmus, ampulla, fimbriae, abbreviated as I, A and F) by specifying the anatomic site as the desired sample attribute. Four main functions are available: 1) Function ‘getClusterSummaryReplot()’ generates box plots to inspect cellular compositions across attribute groups. For example, users can visualize cellular compositions across anatomic sites (I, A and F) ( Figure 1C) and across patients (D1-7) ( Figure 1D). 2) Function ‘getGoiDotplot()’ generates dot plots for a pre-defined set of gene markers, inputting through a separate file. This plot can visualize normalized and center scaled gene expression values across attribute groups, with dot color representing the center scaled average gene expression of cells with each group and dot size representing the percentage of cells with non-zero expression. Gene annotations can be specified and displayed on top as a legend bar (as shown in Figure 1E and 1F). 3) Function ‘getGoiHeatmap()’ generates heatmaps for a pre-defined set of gene markers, which can be specified either in a separate file or through the option ‘geneNames’. This heatmap can visualize gene expression at the single-cell level ( Figure 1G) or the donor level ( Figure 1H) across groups, providing complementary information for the dot plot that display scaled average values. 4) Function ‘getScatterPlot()’ generates scatter plots for any two selected cell types or attribute groups via the option ‘selectedGroups’.

‘Step 4’ can perform downstream analysis including: 1) differential expression (DE) analysis by the function ‘getclusterExpCondDe()’; 2) over-expressed gene identification via the function ‘getExpCondClusterMarkers()’; 3) pseudo-time trajectory analysis via the function ‘getExpCondClusterPseudotime()’; and 4) sub-cluster identification via the function ‘getHippoRes()’. All analyses can be applied to any specified cell types, any specified attribute groups, or a combination of cell types and attributes. DE analysis can be conducted at single-cell level or pseudo-bulk level. The options for single-cell based DE include wilcox, MAST and t-test. For pseudo-bulk-based DE, expression levels are averaged and normalized for each donor with respect to the total number of cells, followed with DE test using DESeq2.⁹ The pseudo-time trajectory analysis in scRICA is performed with slingshot.¹⁰ Users can select any subset of cell types from specified attribute group for pseudo-time analysis. Sub-clustering identification is conducted by HIPPO,¹¹ which implements a zero-inflation test that runs iteratively to select heterogeneous features in order to refine and identify biologically important sub-cluster cell types.

Input and output structures

scRICA’s input includes a metadata table and a processed count matrix. The processed count matrix can be in different formats, including a processed count matrix saved in a folder, which is directly generated from the 10X genomic Cell Ranger analysis tool (https://support.10xgenomics.com/single-cell-gene-expression/software/overview/welcome), or a count matrix stored in text or hd5 format. The input metadata table includes two mandatory columns: the first column specifies all sample names, and the second column (‘path’) specifies the full path of the count matrix for each corresponding sample in a row. Optional columns include: 1) ‘doubletsRmMethod’, which specifies the doublets removal algorithm metroid, detroid,⁶ or both (designated by OL); 2) ‘filterFname’, which provides the full path of Excel files that specify cells to be removed from the analysis; and 3) ‘expCond*’, which allows the inputting of multiple sample attributes for each sample. These sample attributes can be inherited through the entire analysis workflow by setting the option ‘expCondCheck’ properly.

scRICA outputs all analysis results to a main folder named ‘scRICA_results’ by default or it can be specified by user via the option ‘resDirName’. Results from each step are organized into corresponding sub-folders ( Figure 1B). The integration analysis results are saved as an RDS object. All visualizations are saved as PDF files. QC summary, DE analysis, and trajectory analysis are saved as both text files and Rdata objects, making them easily accessible for additional analysis.