Building the mega single-cell transcriptome ocular meta-atlas

A plethora of single-cell transcriptome studies in the retina

The retina contains a multitude of cell types that, in total, are responsible for turning light information into signal for the brain to interpret as vision. Briefly, the photoreceptors (rods and cones) are responsible for capturing the photons. The retinal pigmented epithelium (RPE) behind the photoreceptors physically supports the rods and cones by processing by-products of the visual cycle. Müller glia serve as support cells for the neurons. The retinal bipolar cells transmit the electrical signal from the photoreceptors to the retinal ganglion cells. Horizontal cells (HCs) and amacrine cells (ACs) regulate and help interpret signals from the photoreceptors. The signal is relayed via the retinal ganglion projections through the optic stalk to the brain (for review see [1]). Since 2000 many groups have investigated gene expression in small numbers of individual cells of the retina [2–5].

The recent introduction of lower cost and high-throughput single-cell sequencing technology has led to an explosion of research across many fields. As of early 2021, >40 million cells have been sequenced across >1,200 studies and the mean size of studies starting in 2020 is >100,000 cells [6]. The retina was used as the source tissue in one of the earliest works in the high-throughput single-cell transcriptomics field [7]. As of late 2020, >20 published studies, cumulatively containing >1 million cells, have used single-cell technology to profile cell type–specific gene expression patterns, cell fate trajectory, tissue and cell differentiation, and disease perturbation across multiple mammalian species [7–29].

While the gene–cell count tables are generally made available in repositories such as the Gene Expression Omnibus (GEO), there are no requirements to uniformly process the data. This means that the count tables cannot be used in cross-study comparisons because even small differences in the computational pipeline (e.g., aligner, transcriptome reference) create study-specific effects. This issue can be addressed only by re-quantifying the data in a uniform computational environment. Fortunately, owing to the continued development of computationally lightweight gene quantification tools in the single-cell space (e.g., kallisto bustools, alevin-fry), re-quantification does not require massive compute and time resources [30, 31].

Still, even after re-quantification under identical computational conditions there remain study-specific batch effects due to the diversity in single-cell technologies used and variation in tissue handling and processing across scientific groups. The single-cell community have recognized that removal of these technical (also referred to as batch) effects is a critical issue and have independently developed many tools, although it remains unclear which tools and parameters are optimal for a particular dataset [32–43]. Reinforcing this point, a couple of groups have quantified performance of multiple methods in a consistent framework across several test datasets, finding no consistently best approach [44, 45].

The projectable meta-atlas

We propose that by re-processing publicly available raw single-cell transcriptome data in a consistent bioinformatic framework and optimally using batch correction tools we can create a meta-atlas of retina single-cell transcriptomes. Because there are thousands of possible permutations of single-cell tools, references, and parameter choices, we create our meta-atlas (which we refer to as the single-cell Eye in a Disk [scEiaD]) by benchmarking integration outcomes across multiple important single-cell RNA-sequencing (scRNA-seq) processing parameters (e.g., batch removal method, number of hyper-variable genes [HVGs], clustering resolution). The benchmarking system that we developed uses a wide variety of metrics that combine in the R package scPOP (single-cell Pick Optimal Parameters). The scEiaD will be of utility to 2 communities. First, the ocular community can both search scEiaD for gene expression across many dimensions (e.g., cluster, cell type, study) and project their own single-cell data onto scEiaD for comparison and rich automatic cell labeling. Second, the computational community can use this very large, well-curated dataset to test algorithms for compute efficiency and performance in a diverse environment. Because we believe that data reuse is a powerful and efficient approach to facilitate discovery, we provide our meta-atlas code-base and the meta-atlas in several data formats and propose general guidelines to optimally create custom meta-atlases.

We identify 33 ocular scRNA datasets across 3 species

The first step in building a meta-atlas is identifying studies from which to draw the data. We identified ocular scRNA-seq studies by querying PubMed, the SRA, and the European Nucleotide Archive (ENA) for the inclusive terms “retina,” “single cell,” “scRNA,” “ocular,” “eye,” and “transcriptome.” We then hand-filtered the results to keep only ocular and normal (non-perturbed or mutagenized) data from scRNA-seq technology. In December 2020 we identified 33 deposited datasets that have been published in 27 publications (Fig. 1). To provide a non-ocular reference we also downloaded the raw sequence data from the Tabula Muris project for re-processing. In cases where the fastq file from the SRA was not processed properly (always 10x v2 or v3), we acquired the 10x bam files (via SRA or personal correspondence) and re-extracted the fastq. After downloading all the data we had 11 TB across 2,470 fastq file sets.

Transcriptome quantification across multiple technologies

Droplet- and well-based scRNA-seq technologies require different quantification approaches because the former have Unique Molecular Identifiers and multiple cells are quantified within a single file. We wrote a Snakemake-based pipeline (SnakeQUANT, see Methods) to quantify and merge both droplet- and well-based technologies into a single matrix for downstream processing. For well-based data we perform both gene- and isoform-level quantification; for droplet-based technologies we quantify both exonic and intronic gene-level expression to facilitate calculation of RNA velocity. In total we quantify 6.7e+10 molecules, finding 1.1e+10 unique molecules with a mean pseudoalignment rate of 66.5%. Across the 766,615 cells (post–quality control) we have a mean of 3,410 RNA counts across 978 unique genes (see Supplemental File kallisto_stats.tsv and splicing_stats.tsv for more details).

A total of 1,204,269 cells before quality control

Gene-level counts were quantified with the kallisto bustools pseudo-aligner for both the droplet- and well-based samples. After empty droplet removal, we had 1,204,269 cells. We then removed cells that had >10% mitochondrial reads across all gene counts or <200 unique genes quantified or wsere identified as an in silico doublet (see Methods). After these quality control (for review see [46]) steps we were left with 766,615 cells (Supplemental Fig. 1).

A core objective of many scRNA-based studies is labeling the cell types. Because this information is crucial to assess dataset integration and provide an accurate reference for user querying, we extracted individual cell labels with a combination of inspecting the GEO website, supplemental information from the publication, web resources (e.g., a web app was created for the article), and personal correspondence. After normalizing cell type name nomenclature, we obtained labels for 375,966 cells across 33 cell types (Supplemental Table 1).

Running 11 tools in a Snakemake-based system

Disentangling the technical and biological effects when integrating multiple datasets is crucial. We define "batch" as each unique biological sample and assume that each study is ≥1 unique sample. We studied the metadata and methods of each study to identify the unique biological samples. Within the current scEiaD dataset we identified 86 batches across 33 deposited datasets in 26 published articles.

A wide variety of methods have been written for scRNA-seq integration. Because we were uncertain which would perform the best, we ran 11 tools with a commonly used set of key parameters such as number of HVGs, number of latent dimensional space, and the number of nearest neighbors for the Louvain clustering algorithm. The Snakemake system was used to automate the running of the wide variety of tools. In total 5,591 jobs were run to quantify gene expression and build the unified Seurat objects (SnakeQUANT) and 2,446 jobs were run to assess integration performance (SnakePOP).

The 2 key metrics that have to be balanced to optimize integration performance are cell type or cluster purity (where different cell types or clusters should be homogenous) and batch mixing (the same cell types should be similar across independent studies). While these can be visually assessed by looking at marker gene expression across the 2D UMAP projection, it is more rigorous and scalable to quantify these diametrically opposed characteristics.

scPOP wraps several different methods for measuring integration performance

Multiple methods have been proposed to quantitatively evaluate batch correction. Some of these metrics evaluate the concordance between sets of labels, while others compute distances between the individual data points of a given set of labels. While any one of these methods can be useful, we propose that calculating and evaluating them in tandem provides greater accuracy for dataset integration. We developed the R package scPOP, a lightweight, low-dependency R package that brings together the Local Simpson Index (LISI), Average Silhouette Width (ASW), Adjusted Rand Index (ARI), and Normalized Mutual Information (NMI) metrics from the R packages Harmony, kBET, and aricode, respectively. The LISI and ASW were used to measure batch mixing (where lower is better), cell type mixing (higher is better), and cluster mixing (higher is better). NMI and ARI were used to assess the consistency of cell type to cluster assignment (where 1 is perfect correspondence between cluster and cell type).

To visualize the interplay between batch mixing and cell type distinction we plot the batch mixing LISI score (which has been multiplied by −1) on the y-axis (higher is better) against the cluster LISI on the x-axis (higher is better). The best performer on both metrics will be in the top right corner (Supplemental Fig. 2A). In the same manner we plot the silhouette metric (Supplemental Fig. 2B). To merge the different scores we define , where m is a metric (LISI by batch, LISI by cluster, LISI by cell type, silhouette by batch, silhouette by cluster, silhouette by cell type, NMI, and ARI). Weight is set to 1 but can be either explicitly set or randomly chosen should it be desired to change the influence of certain metrics.

On 1 extreme we have ComBat, which merges together different batches very well but also mixes together the distinct cell types (Fig. 2A). The other extreme is not using any batch integration method, where you see very distinct groups of cells, but also each study has a distinct region in the UMAP (Fig. 2B). In our scEiaD dataset we see that methods like ComBat, Harmony, scArches, and CCA are weighed more towards batch mixing than cluster and cell type purity. Scanorama and bbknn prioritize cleanly separating the clusters. With our scEiaD meta-atlas insct, trVAE, and magic do not perform particularly well in batch mixing or cluster purity. Overall, CCA, fastMNN, and scVI perform best on this particular dataset. For further qualitative investigation of the integration performance we provide the UMAP visualizations of each integration method, normalization, and latent dimensions colored by cell type, study accession, or organism as supplementary files.

As the LISI and silhouette metrics provide independent cluster and cell type (“purity”) and batch (“mixing”) scores we looked to see whether the sumZScale scoring is highly influenced by changing the weight placed on purity or mixing by multiplying either by a multiplier of 3 (Supplemental Fig. 3). While scVI still performed well no matter the weight chosen, there were some changes when more weight was placed on batch mixing. Most notably fastMNN with 8 latent dimensions performed better than fastMNN with 30 latent dimensions. In a reversal of the fastMNN results, CCA with 30 latent dimensions performed better than CCA with 8 latent dimensions. We also bootstrap random weights across each individual metric and get similar results (Supplemental Fig. 4). Across the 11 integration methods we tested, scVI achieved the highest sumZScale for our scEiaD meta-atlas and furthermore is a desirable choice because of its very short run time (scVI can complete in <1 hour while fastMNN takes >12 hours and CCA takes >24 hours).

Different normalization methods alter integration performance

There are several normalization approaches that have been used or published. The “standard” approach that the popular analysis packages Seurat and scanpy use by default is to, per cell, divide the counts by the sum counts for the cell, multiply by a scaling factor, then log transform. This helps make the count distribution more normal, which is an assumption that many algorithms require. In contrast, the scran normalization method groups cells into pools and normalizes across the pool summed counts instead of the individual cell counts. We also use the square root (sqrt) normalization, which replaces the log transformation with a square root. Library size (libSize) normalization omits the sqrt or log transformation. Finally some methods, like scVI, directly use the raw counts for modeling the data.

As expected the libSize normalization, which omits the log or square root scaling, generally performs the worst (Fig. Cc). We see that the remaining normalization techniques alter the batch correction performance, although the outcome differs across the different methods. We also see that changing the number of latent dimensions (8 vs 30) can occasionally dramatically change performance. These results demonstrate the importance of assessing performance in a rigorous manner across many parameters.

Further optimization of scVI with grid search and projection

To find the best set of parameters for scVI in our dataset we did a grid search across key parameters: HVG, latent dimensions, and k nearest neighbors (knn). Furthermore we used a recent advance in scVI capability (≥version 0.8.0) adapted from scArches that allows one to build a reference model and query or project the new cells onto it [40]. We refer to this projectable model as “scVI-projection” and to the previous scVI model as “scVI-standard.” We built a scVI-projected model trained on human cells and then projected the mouse and macaque data onto it. We then compared scVI-projection against the previous scVI-standard across all the aforementioned parameters. Using scPOP we first saw that the scVIprojection approach generally performed better than running scVI with all of the data (Fig. 2D). We found the optimal overall parameters to be 5,000 HVG, 8 latent dimensions, and with 5 knn for the cluster finding. We also varied the UMAP projection values of nearest neighbors and minimum distance to qualitatively pick a 2D projection, selecting a minimum distance of 0.1 and 50 nearest neighbors.

High-accuracy xgboost ML model built to label unknown cell types across all technologies

To further study cell type–specific expression patterning we needed to label the 414,106 unlabeled cells. Traditionally this is done by clustering the cells, then using cell type–specific markers to label the clusters. However, because we had hundreds of thousands of expert-labeled cells across 17 publications (Fig. 1A) we built a xgboost-based machine learning (ML) model that used two-thirds of the labeled cells as a training set (see Methods for more details) to train a cell type predictor for scEiaD input. The trained model was used to predict the cell type assignments for all cells in scEiaD. In this manner we both label most cells (cells that cannot be assigned a cell type with a probability >0.5 were left unlabeled) and correct a small number of probable mislabels in the truth set (Fig. 3A).

As a brief case study of our ML performance, we look at the cell type labels assigned to the study by Shekhar et al., where they used SMART-seq on a retinal bipolar cell–enriched population of cells [21]. Even though our ML algorithm was trained only on droplet-based data, our algorithm labels most of this dataset as retinal bipolar cells, with the next most common cell types being AC. The same result was found by Shekhkar et al. (Supplemental Fig. 5). To more generally evaluate performance we use the precision recall (PR) curve, which visualizes the ability of the model to precisely label known cells at a given confidence. The area under the PR curve (AUC) summarizes the effectiveness of the model across different cell types, with 1 being the highest performance. The xgboost model can predict rods, bipolar cells, and Müller glia with near perfect performance (Supplemental Fig. 6). Most of the remaining cell types can be predicted with an AUC >0.9 (Supplemental Table 2). Several of the precursor cell types (photoreceptor, AC/HC, and neurogenic) were labeled with lower confidence (Supplemental Table 3). The next most common labels for these cell types were either other precursor cells (e.g., 100 AC/HC precursors were labeled as neurogenic) or the adjacent terminal cell type (e.g., 254 photoreceptor precursors were labeled as cones). Other cell types that were challenging for the model to predict were the artery, choriocapillaris, and vein. Artery, choriocapillaris, and vein are constructed from endothelial cells, and we find that for all 3 of these, endothelial was the second most common label. Overall, our xgboost-based ML model shows strong accuracy across the major cell types of the retina (overall AUC of 0.98).

ML cell type labels result in high study diversity for each cell type

After ML projection of cell type labels from the original 375,966 labels onto a total of 758,278 cells we have substantially improved the number of studies per cell type. For example, we went from 8 human studies with labeled Müller glia to 14 after labeling (Supplemental Table 4). Overall we go from a mean of 4 studies per human cell type to 8 and from 2 studies per mouse cell type to 8 after transferring the cell type labels.

As another check on the quality of the cell type assignments, we ran the cell and cluster-independent haystack gene search and pairwise differential expression tests between the predicted cell types (see Methods for further details) [47]. We show that the 5 most differentially expressed genes for each of the major retina cell types are consistent with known retinal cell markers (Fig. 3A). As a simple metric to identify known and unknown genes relating to the cell type–specific expression we search PubMed for the number of publications with 2 searches per gene. We expect most of the genes identified to be known in the literature. The first search is the more precise “gene AND cell type” (e.g., “PDE6H AND Cones”) and the second search is the more inclusive “gene AND retina” (e.g., “PDE6H AND Retina”). Of the 50 genes (in Fig. 3A), 37 had ≥1 citation in the gene–cell type search (Supplemental File celltype_markers.tsv) and 45 had ≥1 citation in the inclusive search. The 50 genes had a mean of 46 studies (with the inclusive gene by “retina” search). In contrast, 100 randomly chosen genes had a mean of 2 (Wilcoxon test P < 1.44 × 10⁻¹⁷).

The scVI-based scEiaD UMAP projection blends batches and species while separating cell types

The 2D UMAP projection of the scVI-calculated batch-corrected 8 latent dimensional space blends the 33 studies together while also maintaining distinct space for the 31 unique cell types. We also see good mixing across all the droplet- and well-based single-cell technologies (Supplemental Fig. 7). We see the neurogenic and progenitor populations from which the retinal cell types are derived near the center of the UMAP visualization. The photoreceptor precursors are adjacent to the neurogenic population and, as demonstrated by the RNA velocity dynamics, flow into the rods and cones. The AC and HC precursors (AC/HC) likewise flow from the neurogenic center into the mature ACs and HCs.

The photoreceptors (cones and rods) of the retina, which are responsible for color and low-light vision, respectively, are near each other in the UMAP space. The major remaining retina cell types, by proportion in the mammalian eye, are the Müller glia, which are a glial cell type that help support the neurons of the retina. Next we have the neural cell types that transmit and help interpret the signals from the photoreceptors before they leave the retina via the optic stalk: the ACs, retinal ganglia, HCs, and bipolar cells. All of these cells are in well-separated spaces in the UMAP. Finally we see across species that the major cell types overlap each other (Fig. 3D). The macaque retina cells are not present in the precursor/neurogenic center of the UMAP as expected because only fully developed tissues were sampled in these data.

While by eye the UMAP 2D projection generally blends together the 3 different species (macaque, mouse, human) in the UMAP visualization, we more rigorously tested this by changing the inputs to our xgboost cell type predictor ML system. We trained 2 new models: 1 with only human data and 1 with only mouse data. We then applied this model to the other 2 species to see whether the model trained on 1 species was generally transferable to another species. Again, both models (human only and mouse only) proved to be highly accurate at predicting cell types in the other species with the terminal cell types (Supplemental Fig. 8).

Projection of outside data onto scEiaD demonstrates similarities and differences of induced pluripotent stem cells and organoids to primary cells

The hard work of creating this resource can be leveraged and extended by the wider community with a few relatively simple steps and modest compute requirements. Very briefly, if outside groups quantify their mouse or human scRNA with kallisto (bustools) and the same references (see Methods), they can overlay their data on top of scEiaD by (i) installing scVI (version 0.9.0 or higher), (ii) downloading our 13-MB scVI model, and (iii) following the Jupyter notebook on Google colab that we provide as a live demo [69]. We demonstrate the power of this approach in 2 ways.

First, a wild-type retinal organoid dataset from Kallman et al. was projected onto scEiaD [48]. In the Google colab notebook we show how a subset of the reads from Kallman et al. SRR12130660 can be processed from the raw reads to a UMAP visualization in <10 minutes. Indeed we see how organoid-based retinal cell types can be detected overlapping primary cells for RPCs, photoreceptors, ACs, and retinal ganglion cells (Supplemental Fig. 9).

Finally, we demonstrate how the RPE derived from the Bharti induced pluripotent stem cell (iPSC) differentiation process express canonical RPE markers of TTR and RPE65 but end up in a slightly different position in the UMAP projection (Fig. 4A). Differential expression analysis (see Methods) identified that vimentin is nearly exclusively expressed in the iPSC-based RPE (Fig. 4B). This result is not surprising because Hunt and Davis demonstrated that cultured and proliferating RPE express high levels of vimentin [49] and this same trend is seen in bulk RNA-seq from stem cell RPE and primary RPE (Supplemental Fig. 10) [50].

Quantification of gene counts

Gene quantification is handled by the SnakeQUANT snakefile. First, we generated multiple quantification indices to facilitate calculation of RNA velocity. For each droplet technology, we generated a separate set of transcript sequences that contain both exonic and intronic sequences using the “get_velocity_files” function from the R package BUSpaRse. The reference transcriptome annotations used to build sets of transcript sequences were the Gencode “gencode.vM25.annotation.gtf.gz” and “gencode.v35.annotation.gtf.gz” for mouse and human, respectively [51, 52]. Because the Macaca fascicularis genome is less well annotated we used the Ensembl release 101 genome and transcriptome annotation “Macaca_mulatta.Mmul_101.gtf.gz” [53]. A single set of exonic transcript sequences was created for each species for all well techonologies. A kallisto quantification index was generated for each of these fastqs using kallisto index (0.46.2) [54]. Well-based samples were quantified using kallisto quant. For the relatively few single-ended samples, the parameters “–single -l 200 -s 30” were used for kallisto quant. Otherwise the “–bias” flag was added. For the droplet-based samples, we adapted the bustools workflow for generating spliced and unspliced count matrices (0.39.4) [55].

Intersection of gene names between mouse, macaque, and human

To facilitate comparison of gene expression across species, where possible we converted mouse and macaque gene ids and names to human ones. We downloaded a mapping of orthologous genes between human, mouse, and macaque using the Ensembl BioMart web browser in November 2020. We identified 15,759 human genes that could be directly mapped to mouse and macaque orthologs. Genes present in mouse or macaque that were not found in human were not used for HVG gene selection but were retained and used for differential gene expression.

Custom macaque reference quantification

Because we noticed that several retina marker genes (e.g., NRL and CRX) had very low expression in the macaque data, we quantified the scRNA data twice: once with the Ensembl reference and again with the same Gencode human reference used for the human data. We compared the gene-level counts for each cell and replaced the macaque gene count with the human counts if the human counts were greater than the macaque counts, and, to prevent genes with very few total counts from being used, we required the counts to be greater than the first quartile of non-zero macaque gene expression.

Removal of empty droplets and further quality control

After bustools count, we used R (3.6.2) to remove empty droplets. The BUSpaRse package was used to input the bustools counts mtx file. The DropletUtils package with the “barcodeRanks” function was used to automatically detect the inflection point in the barcode count ranks that delineates the likely empty droplets [56]. We then removed cells with percent mitochondrial reads of >10%. After merging the individual count matrices into 1 sparse matrix, we created a Seurat version 3 object and removed cells with <200 detected unique genes, and for the droplet data, >3,000 detected genes (these are more likely to be doublets) [32].

Normalization and batch effect correction

The following steps were carried out (normalization through benchmarking is handled by the SnakePOP pipeline). We tested several gene count normalization approaches because we were not certain which would produce an optimal outcome: standard (default Seurat, library size normalization, then log transform), sqrt (same, but with sqrt normalization), libSize (omit the log or sqrt normalization), scran, SCT from Seurat, and for scVI, no normalization (counts) [57]. Our R implementation of the normalization approaches, as well as how we constructed the Seurat v3 object, can be found in the supplementary file “make_seurat_obj_functions.R”.

Batch normalization under a grid search procedure

We tested scArches, bbknn, insct, magic, scVI, CCA, scanorama, harmony, fastMNN, combat, and none against 2,000 HVGs, the aforementioned different gene count normalization procedures, and both 8 and 30 outputted batch-corrected latent dimensions. The latent dimensions are the input for clustering, the 2D UMAP visualization, and the xgboost ML to transfer cell type labels to unlabeled cells. We were unable to run every method successfully with every normalization method. For example, magic could not complete with the standard or libSize normalization. Insct was only able to complete the scran normalization. We also tried the DESC, liger, and Conos batch correction methods but were unable to get them to work reliably, so they were dropped. While we attempted to use “default” parameters wherever possible, we had to deviate from this to get Seurat's CCA procedure to complete. CCA reciprocally tries to integrate all batches and even with a subset of cells we regularly got “long vector” errors. We got around this issue by setting the human batches as the “reference.” The batch correction steps implementation can be found in the supplementary file “merge_methods.R” and as a deposit in Zenodo [58].

Clustering and UMAP

Louvain-Jaccard clustering was performed against the batch-corrected latent dimensions used the Seurat implementation [59]. For the all methods benchmarking we used knn of 20. The scVI-only tuning reduced the knn parameters to 5 and 7 (where 5 gives more clusters than 7) to increase the cluster number. We also used the Leiden algorithm as implemented by PARC with a resolution of 0.6 and 0.8 (higher results in more clusters) [60, 61]. These 2 resolutions were chosen because they gave roughly the same number of clusters as the Seurat Louvain-Jaccard approach with a knn 7.

The UMAP visualization was calculated with the Seurat “RunUMAP” using the uwot R package [62]. We tried min.dist parameters of 0.001, 0.1, and 0.3 and tried n.neighbors across 15, 30, 100, and 500. A smaller min.dist value gives “tighter” groupings while a higher number of n.neighbors uses a larger number of near cells to calculate the global positioning.

Benchmarking and scPOP

We wrote the scPOP R package (available via the R package resource CRAN) to unify the LISI and Silhouette metrics from Harmony and kBet, respectively, along with NMI and ARI [36, 63]. LISI and Silhouette require a dense matrix, which is a problem for our data because a 766,615 cell × 8 latent dimension dense matrix cannot fit in our largest available compute node (1.5 TB). We down-sampled the dataset to ∼100,000 cells, taking care to keep all rarer cell types for the LISI and Silhouette benchmarking.

To merge these metrics into a balanced single score, we Z scale each and sum them. scPOP produces both tables and visualizations, allowing the user to quickly see both the interplay of batch mixing and cluster/cell type separation and the overall performance. If a user wishes to prioritize batch mixing or cluster/cell type separation, we let the user provide a custom batch/cluster–cell type scaling value (1 is the default).

Multi-step doublet removal

To identify probable doublets (>1 cell in a droplet) we ran DoubletDetect and scrublet and calculated the distribution of DoubletDetect and scrublet scores across all clusters and removed clusters with a score in both metrics >4 standard deviations above the mean [64, 65]. This removed another 23,457 cells, leaving 766,615 in total.

xgboost-based cell type model

To identify cell types for the 361,456 unlabeled cells we designed a custom xgboost-based cell type classifier. We took labeled data and split them into training (2/3) and test (1/3) sets, stratified by cell type. The input features used to train the model are the scVI latent dimensions, the total number of reads in each cell, the number genes detected in each cell, and the percent mitochondrial gene expression of each cell. We additionally generated features using the age of each sample to group the data into 3 developmental categories (early development, late development, and adult) and then generated a one-hot encoded feature for each category. To speed up training times, we used the gpu implementation of the xgboost algorithm from the the xgboost Python library. The model was trained using default parameters. The trained model had an overall macro and micro AUC score of 0.98 and 0.99, respectively. This model was then used to identify labels for all cells. Unlabeled data were pre-processed identically to training data and fed into the model to generate a vector of label probabilities for each cell. We selected the highest label probability for each cell and required a minimum probability of 0.5 to assign a label to a cell.

For the organism-specific xgboost ML we followed the above procedure, except that we combined early RPCs, late RPCs, and RPCs into 1 category and did not attempt to predict the non-retina cell types (e.g., fibroblasts) because there were very few labeled cells across all 3 organisms.

Marker gene identification

To identify marker genes across the CellType (predict) and cluster groups, we used the scran findmarkers (Wilcox test) along with the singleCellHaystack algorithm [47]. The scran findmarkers test runs a Wilcox test in a pairwise manner (e.g., rods vs all other cell types). It returns an overall P-value (and false discovery rate [FDR]) that assesses how good the gene is at separating the group of interest from all other cells. It also returns for each pairwise comparison an AUC score, where 1 is a perfect power to distinguish and 0 is no power. The singleCellHaystack algorithm uses a Kullback-Leibler divergence (D KL) measurement of the scVI lower dimensional space to identify genes with non-random distribution. A higher D KL score represents a gene with “specific” expression in the lower dimensional space and is used to calculate an FDR-corrected P-value against the full distribution of D KL values. We filtered to keep genes with scran FDR < 1, a mean AUC > 0.2, and a log₁₀(D KL FDR) < −10,000. No more than 50 genes for each cell type were retained (sorted by mean AUC).

Calculation of RNA velocity

RNA velocity calculations were performed with the Velociraptor wrapping of the scVelo Python library [66]. From the anndata objects generated by our Snakemake pipeline we calculated velocity across all genes. Genes without detectable velocity were dropped. The scVI-generated latent space (instead of principal component analysis) was used to calculate first- and second-order moments. The calculated moments were used to estimate RNA velocity. Differential velocity was tested between cell types using pairwise Wilcox rank sum tests.

Limitations

The scVI model is first built on human data. The mouse and macaque data are then projected (or queried) onto it with the scVI implementation of the scArches method. While this system works very well to integrate information between these 3 species, this approach may not scale to more distantly related species. Another limitation is that discrepancies between cell type labels between different laboratories make certain transitioning cell type labels a bit imprecise. One example is how the rods and photoreceptor precursor labels partially overlap. Although we attempted to ameliorate the issue by removing cell type labels in large disagreement with the consensus, some disagreements could propagate into our machine-labeled cell type assignments. These issues may result from labeling continuous processes with discrete labels.

While scEiaD distinguishes the major cell types very well, some of the cell types contain many “subtypes”—notably the ACs have a huge variety in morphology, with a single-cell–based study identifying >60 different types of ACs in mouse [29]. At this time our batch-corrected pan-retina cell space does not precisely resolve these sub-cell types with high resolution. We are actively working to “sub-cluster” the cell types so we can robustly and reliably identify the high diversity of retinal cell types across the entire retina.

The scEiaD is a unique ocular resource that provides a highly diverse, large N dataset with a relatively small amount of compute power

We have assembled the largest ocular single-cell transcriptome database to date. The rapid advancement of algorithms to batch-correct and process data continues to reduce the computational requirements to handle huge numbers of cells. The scVI batch correction step with ∼1 million cells runs within 1 hour on a GPU and 150 GB of memory. This places this crucial step within the capabilities of a moderately powerful computer or a cloud compute node. Further downstream processing can largely be done on computers with ≥64 GB of memory and a few hundred GB of disk space. We believe that our efforts can be replicated in any other tissue/system with a large number of independent studies by a small number of computational scientists following our general approach. We provide the completed analysis as both Scanpy (h5ad) and Seurat objects (Supplemental Table 5).

Benchmarking and quantitation of integration performance is crucial for meta-atlas studies

We originally intended to use the Seurat CCA method to integrate the datasets. However, the long run-times of CCA and poor integration of our known cell types led us to benchmark more methods and parameters. After adding more integration methods we first attempted to “hand-assess” the integration results by using the UMAP 2D projection view. This proved to scale poorly, and this led us to curate some of the more useful benchmarking algorithms (NMI, ARI, silhouette, and LISI) that roughly matched our hand-assessed results into the scPOP R package. While we chose the scVI algorithm, we strongly suggest any other groups attempting a similar meta-atlas construction choose a quantifiable set of criteria so that optimal methods and parameters can be picked. We found in our analysis that substantial differences in the integration performance can come from altering the number of HVGs and the number of latent dimensions.

Transfer of cell type labels from a smaller number of studies onto the remaining cells is a powerful way to increase diversity in a meta-atlas

We first hand-curated >350,000 published cell type labels across 24 publications. With a xgboost algorithm using the latent dimensions, cell age classification (developing or matured), and the UMAP coordinates, we can very accurately label the remaining cells. Many other cell type labeling algorithms and systems exist for those groups less willing or able to tune an ML algorithm. For example, the developers of scVI also have a cell type label projection algorithm called scANVI [67]. Whichever approach is used, taking a smaller number of high-quality labels and projecting them onto the remaining cells is a powerful way to leverage community knowledge across a huge diverse dataset.

Projection allows community knowledge to be leveraged by all

Many retina atlases have been published to date. We argue that we have created the first atlas that is generally useful because (i) our dataset/atlas is several times larger than any other published set, (ii) our data are available via download in several forms at [69] and Zenodo accession 5129265, and crucially (iii) we provide a Google colab/Jupyter notebook that step-by-step demonstrates how to use scVI to project (or query) outside data onto our scEiaD with minimum compute resources. We demonstrate concretely how this can work by showing how iPSC-based RPE can be queried onto the reference dataset to demonstrate both similarities and dissimilarities in their transcriptomes.