The Transcriptome of SH-SY5Y at Single-Cell Resolution: A CITE-Seq Data Analysis Workflow

Abstract

Cellular Indexing of Transcriptomes and Epitopes by Sequencing (CITE-seq) is a recently established multimodal single cell analysis technique combining the immunophenotyping capabilities of antibody labeling and cell sorting with the resolution of single-cell RNA sequencing (scRNA-seq). By simply adding a 12-bp nucleotide barcode to antibodies (cell hashing), CITE-seq can be used to sequence antibody-bound tags alongside the cellular mRNA, thus reducing costs of scRNA-seq by performing it at the same time on multiple barcoded samples in a single run. Here, we illustrate an ideal CITE-seq data analysis workflow by characterizing the transcriptome of SH-SY5Y neuroblastoma cell line, a widely used model to study neuronal function and differentiation. We obtained transcriptomes from a total of 2879 single cells, measuring an average of 1600 genes/cell. Along with standard scRNA-seq data handling procedures, such as quality checks and cell filtering procedures, we performed exploratory analyses to identify most stable genes to be possibly used as reference housekeeping genes in qPCR experiments. We also illustrate how to use some popular R packages to investigate cell heterogeneity in scRNA-seq data, namely Seurat, Monocle, and slalom. Both the CITE-seq dataset and the code used to analyze it are freely shared and fully reusable for future research.

Keywords: CITE-seq; gene regulatory networks; neuroblastoma; single-cell; transcriptomics; unsupervised learning.

Figure 1

Exploratory analysis of single-cell mRNA expression on SY5Y cell line. (A) SY5Y cells were left growing in flask until reaching about 90% of confluence before splitting to two flasks for library preparation and sequencing. (B) Scatterplot showing Log10 counts of HTOs. Three populations are distinguished on the basis of Hashtag counts: cells belonging to SplitA, cells belonging to SplitB, and multiplets. Proportions of UMI counts, measured genes, mtUMI and mitoRatio are reported in the table, and QC plots summarizing these metrics are shown in (C) barplot, showing the number of cells assigned to each sampling group, (D) density plot showing a similar number of UMI/cell in both SplitA and B, while multiplets contained higher UMI/cell, and (E) Box and Whisker plots showing the number of genes detected per cell in the three groups.

Figure 2

Exploratory analysis of top varying and most stable genes across the dataset. (A) Seurat standardized variance vs. average expression plot. The top 10 most varying genes are indicated. Most average expressed genes (x-axis) with lowest average variance (y-axis) are indicated in (B) SplitA, (C) SplitB, and (D) in the entire dataset. Genes showing an average variance <0.2 and >4 LogNorm Average Expression are labeled, together with commonly used HK genes, such as ACTB, B2M, and GAPDH.

Figure 3

Candidate SY5Y HK genes. Bar plots showing 13 genes selected in Figure 2D compared for LogNorm average expression (A) or LogNorm Average Expression/Variance Ratio (B). Using ratios seems to be a better choice to identify most stable genes suitable as HK. In (C), all genes showing a ratio > 15 are indicated. Forty-nine genes are suitable HK candidates.

Figure 4

Principal Component Analysis. (A) Single cells are plotted along the first two components retaining the highest variance. No difference among the two splits, or multiplets, was detectable. (B) Principal component analysis on cell cycle genes showed a clear separation of cells according to cell cycle genes expression. (C) Partial regression removed the cell cycle variance, maintaining the difference between G1 and G2M/S cells. Partial regression also allowed retaining signals separating non-cycling and cycling cells, while removing differences in cell cycle phases amongst proliferating cells (which are often uninteresting).

Figure 5

Cluster Analysis. (A) UMAP 2D projection of single cells. No detectable difference was observed between splits. (B) Cells can be assigned to four cluster-forming communities applying the Leiden algorithm. (C) Unsupervised trajectory learning. Most variant genes show expression differences following a path from cluster 3 to 2. Cluster 4 remains separated from the others. (D) Top 12 genes characterizing each community. Cluster 4 show the highest expression of TUBA1A and HSP90AA1, while S100A6 was poorly expressed in a minimal fraction of cells (<0.1). Higher PTMA, HMGB2, and H2AZ1 expression characterizes cluster 1.

Figure 6

Slalom output. (A) Graph showing the most relevant factors identified by the f-scLVM model, both annotated in Wiki Pathways (blue) or not annotated (red) (B) Most relevant genes in the WP COPPER HOMEOSTASIS component (C) Most relevant genes in the WP CELL CYCLE component, (D) Most relevant genes in the hidden05 component.