Additional file 2: Supplementary Table 2–3. This file contains the list of cell markers in each of scTyper.db (Table S2) and CellMarker DB (Table S3) and detailed information such as identifier, study name, species, cell type, gene symbol, and PMID.

Single cell RNA-sequencing (scRNA-seq) provides gene expression profiles of individual cells from complex samples, facilitating the detection of cell type-specific marker genes. In scRNA-seq experiments with multiple donors, the population level variation brings an extra layer of complexity in cell type-specific gene detection, for example, they may not appear in all donors. Motivated by this observation, we develop a statistical model named scCTS to identify cell type-specific genes from population-level scRNA-seq data. Extensive data analyses demonstrate that the proposed method identifies more biologically meaningful cell type-specific genes compared to traditional methods.

Additional file 4: Supplementary Data. An example report summary document of scTyper.

Example single-cell RNA sequencing dataset with marker genes for testing cell type annotation

Arsenic exposure via drinking water is a serious environmental health concern. Epidemiological studies suggest a strong association between prenatal arsenic exposure and subsequent childhood respiratory infections, as well as morbidity from respiratory diseases in adulthood, long after systemic clearance of arsenic. We investigated the impact of exclusive prenatal arsenic exposure on the inflammatory immune response and respiratory health after an adult influenza A (IAV) lung infection. C57BL/6J mice were exposed to 100 ppb sodium arsenite in utero, and subsequently infected with IAV (H1N1) after maturation to adulthood. Assessment of lung tissue and bronchoalveolar lavage fluid (BALF) at various time points post IAV infection reveals greater lung damage and inflammation in arsenic exposed mice versus control mice. Single-cell RNA sequencing analysis of immune cells harvested from IAV infected lungs suggests that the enhanced inflammatory response is mediated by dysregulation of innate immune function of monocyte derived macrophages, neutrophils, NK cells, and alveolar macrophages. Our results suggest that prenatal arsenic exposure results in lasting effects on the adult host innate immune response to IAV infection, long after exposure to arsenic, leading to greater immunopathology. This study provides the first direct evidence that exclusive prenatal exposure to arsenic in drinking water causes predisposition to a hyperinflammatory response to IAV infection in adult mice, which is associated with significant lung damage.

Methods Whole lung homogenate preparation for single cell RNA sequencing (scRNA-seq).

Lungs were perfused with PBS via the right ventricle, harvested, and mechanically disassociated prior to straining through 70- and 30-µm filters to obtain a single-cell suspension. Dead cells were removed (annexin V EasySep kit, StemCell Technologies, Vancouver, Canada), and samples were enriched for cells of hematopoetic origin by magnetic separation using anti-CD45-conjugated microbeads (Miltenyi, Auburn, CA). Single-cell suspensions of 6 samples were loaded on a Chromium Single Cell system (10X Genomics) to generate barcoded single-cell gel beads in emulsion, and scRNA-seq libraries were prepared using Single Cell 3’ Version 2 chemistry. Libraries were multiplexed and sequenced on 4 lanes of a Nextseq 500 sequencer (Illumina) with 3 sequencing runs. Demultiplexing and barcode processing of raw sequencing data was conducted using Cell Ranger v. 3.0.1 (10X Genomics; Dartmouth Genomics Shared Resource Core). Reads were aligned to mouse (GRCm38) and influenza A virus (A/PR8/34, genome build GCF_000865725.1) genomes to generate unique molecular index (UMI) count matrices. Gene expression data have been deposited in the NCBI GEO database and are available at accession # GSE142047.

Preprocessing of single cell RNA sequencing (scRNA-seq) data

Count matrices produced using Cell Ranger were analyzed in the R statistical working environment (version 3.6.1). Preliminary visualization and quality analysis were conducted using scran (v 1.14.3, Lun et al., 2016) and Scater (v. 1.14.1, McCarthy et al., 2017) to identify thresholds for cell quality and feature filtering. Sample matrices were imported into Seurat (v. 3.1.1, Stuart., et al., 2019) and the percentage of mitochondrial, hemoglobin, and influenza A viral transcripts calculated per cell. Cells with < 1000 or > 20,000 unique molecular identifiers (UMIs: low quality and doublets), fewer than 300 features (low quality), greater than 10% of reads mapped to mitochondrial genes (dying) or greater than 1% of reads mapped to hemoglobin genes (red blood cells) were filtered from further analysis. Total cells per sample after filtering ranged from 1895-2482, no significant difference in the number of cells was observed in arsenic vs. control. Data were then normalized using SCTransform (Hafemeister et al., 2019) and variable features identified for each sample. Integration anchors between samples were identified using canonical correlation analysis (CCA) and mutual nearest neighbors (MNNs), as implemented in Seurat V3 (Stuart., et al., 2019) and used to integrate samples into a shared space for further comparison. This process enables identification of shared populations of cells between samples, even in the presence of technical or biological differences, while also allowing for non-overlapping populations that are unique to individual samples.

Clustering and reference-based cell identity labeling of single immune cells from IAV-infected lung with scRNA-seq

Principal components were identified from the integrated dataset and were used for Uniform Manifold Approximation and Projection (UMAP) visualization of the data in two-dimensional space. A shared-nearest-neighbor (SNN) graph was constructed using default parameters, and clusters identified using the SLM algorithm in Seurat at a range of resolutions (0.2-2). The first 30 principal components were used to identify 22 cell clusters ranging in size from 25 to 2310 cells. Gene markers for clusters were identified with the findMarkers function in scran. To label individual cells with cell type identities, we used the singleR package (v. 3.1.1) to compare gene expression profiles of individual cells with expression data from curated, FACS-sorted leukocyte samples in the Immgen compendium (Aran D. et al., 2019; Heng et al., 2008). We manually updated the Immgen reference annotation with 263 sample group labels for fine-grain analysis and 25 CD45+ cell type identities based on markers used to sort Immgen samples (Guilliams et al., 2014). The reference annotation is provided in Table S2, cells that were not labeled confidently after label pruning were assigned “Unknown”.

Differential gene expression by immune cells

Differential gene expression within individual cell types was performed by pooling raw count data from cells of each cell type on a per-sample basis to create a pseudo-bulk count table for each cell type. Differential expression analysis was only performed on cell types that were sufficiently represented (>10 cells) in each sample. In droplet-based scRNA-seq, ambient RNA from lysed cells is incorporated into droplets, and can result in spurious identification of these genes in cell types where they aren’t actually expressed. We therefore used a method developed by Young and Behjati (Young et al., 2018) to estimate the contribution of ambient RNA for each gene, and identified genes in each cell type that were estimated to be > 25% ambient-derived. These genes were excluded from analysis in a cell-type specific manner. Genes expressed in less than 5 percent of cells were also excluded from analysis. Differential expression analysis was then performed in Limma (limma-voom with quality weights) following a standard protocol for bulk RNA-seq (Law et al., 2014). Significant genes were identified using MA/QC criteria of P < .05, log2FC >1.

Analysis of arsenic effect on immune cell gene expression by scRNA-seq.

Sample-wide effects of arsenic on gene expression were identified by pooling raw count data from all cells per sample to create a count table for pseudo-bulk gene expression analysis. Genes with less than 20 counts in any sample, or less than 60 total counts were excluded from analysis. Differential expression analysis was performed using limma-voom as described above.

Currently most methods take manual strategies to annotate cell types after clustering the single-cell RNA sequencing (scRNA-seq) data. Such methods are labor-intensive and heavily rely on user expertise, which may lead to inconsistent results. We present SCSA, an automatic tool to annotate cell types from scRNA-seq data, based on a score annotation model combining differentially expressed genes (DEGs) and confidence levels of cell markers from both known and user-defined information. Evaluation on real scRNA-seq datasets from different sources with other methods shows that SCSA is able to assign the cells into the correct types at a fully automated mode with a desirable precision.

Single cell RNA-sequencing dataset of peripheral blood mononuclear cells (pbmc: T, B, NK and monocytes) extracted from two healthy donors.

Cells labeled as C26 come from a 30 years old female and cells labeled as C27 come from a 53 years old male. Cells have been isolated from blood using ficoll. Samples were sequenced using standard 3' v3 chemistry protocols by 10x genomics. Cellranger v4.0.0 was used for the processing, and reads were aligned to the ensembl GRCg38 human genome (GRCg38_r98-ensembl_Sept2019). QC metrics were calculated on the count matrix generated by cellranger (filtered_feature_bc_matrix). Cells with less than 3 genes per cells, less than 500 reads per cell and more than 20% of mithocondrial genes were discarded.

The processing steps was performed with the R package Seurat (https://satijalab.org/seurat/), including sample integration, data normalisation and scaling, dimensional reduction, and clustering. SCTransform method was adopted for the normalisation and scaling steps. The clustered cells were manually annotated using known cell type markers.

Files content:

- raw_dataset.csv: raw gene counts

- normalized_dataset.csv: normalized gene counts (single cell matrix)

- cell_types.csv: cell types identified from annotated cell clusters

- cell_types_macro.csv: cell macro types

- UMAP_coordinates.csv: 2d cell coordinates computed with UMAP algorithm in Seurat

Single-cell RNA sequencing has revolutionized the study of cellular heterogeneity, yet accurate cell type annotation remains a significant challenge. Inconsistent labels, technological variability, and limitations in transferring annotations from reference datasets hinder precise annotation. This study presents a novel approach for accurate cell type annotation in scRNA-seq data using unique marker gene sets. By manually curating cell type names and markers from 280 publications, we verified marker expression profiles across these datasets and unified nomenclatures to consistently identify 166 cell types and subtypes. Our customized algorithm, which builds on the AUCell method, achieves accurate cell labeling at single-cell resolution and surpasses the performance of reference-based tools like Azimuth, especially in distinguishing closely related subtypes. To enhance accessibility and practical utility for researchers, we have also developed a user-friendly application that automates the cell typing process, enabling efficient verification and supporting comprehensive downstream analyses. The desktop application can be accessed at https://omnibusx.com.

Using single-cell RNA sequencing, we provide a comprehensive map of lymph node resident lymphatic endothelial cells, identifying subpopulations, new markers and functions.

By leveraging single-cell transcriptome and T cell receptor (TCR) sequencing, we aimed to track the transcriptional signatures of CAR T cell clonotypes throughout the course of treatment and furthermore identify molecular patterns leading to potent CAR T cell cytotoxicity. The data presented in this study encompass blood and bone marrow samples from patients ≤ 21 years of age with relapsed or refractory B-cell acute lymphoblastic leukemia (B-ALL) participating in the SJCAR19 phase I/II clinical trial (NCT03573700). In brief, patients enrolled in the clinical trial received either 1 x 10^6 (dose level 1) or 3 x 10^6 (dose level 2) per kilogram of body weight following successful generation of autologous CAR T cell products and lymphodepleting chemotherapy. Peripheral blood was drawn from each participant every week until week 4 post-infusion, at week 6 or 8, and month 3 or 6 if feasible. At week 4 post-infusion, blood marrow was also collected from participants. Total T cells (CD3+) were sorted from each post-infusion sample, as well as the pre-infusion CAR T cell products, and processed through 10x Genomics’ single-cell gene expression and V(D)J sequencing platform using the standard protocol. We identified a unique and unexpected transcriptional signature in a subset of pre-infusion CAR T cells that shared TCRs with post-infusion cytotoxic effector CAR T cells. Functional validation of cells with even a subset of these pre-effector markers demonstrated their immediate cytotoxic potential and resistance to exhaustion. Methods Cells were processed using the Chromium Single Cell V(D)J 5' reagents (10X Genomics). T cell receptor V(D)J cDNA was enriched using the Chromium Single Cell V(D)J Enrichment kit for Human T cells. Corresponding libraries were sequenced on the Illumina NovaSeq platform. Sequencing data were processed using CelLRanger v3.1.0 (10X Genomics) with the GRCh38 reference (v3.0.0) modified to include the first 825 nucleotide bases of the CD19-CAR transcript. The resulting gene expression matrices were aggregated, with read depth normalization based on the number of mapped reads. TCR sequences were processed with version 3.1.0 of the GRCh38 V(D)J reference. Aggregated gene expression matrices were analyzed using Seurat (Hao et al, Cell 2021). Cells with fewer than 300 detected genes, more than 4,999 detected genes, with at least 10% of their expression owed to mitochondrial genes, or with no detected CD19-CAR UMIs (unique molecular identifiers) were excluded from downstream analyses. TCR lineages were integrated with gene expression data using shared cellular barcodes. Additional analyses are described in the corresponding manuscript.

Genomic methods are used increasingly to interrogate the individual cells that compose specific tissues. However, current methods for single cell isolation struggle to phenotypically differentiate specific cells in a heterogeneous population and rely primarily on the use of fluorescent markers. Many cellular phenotypes of interest are too complex to be measured by this approach, making it difficult to connect genotype and phenotype at the level of individual cells. Here we demonstrate that microraft arrays, which are arrays containing thousands of individual cell culture sites, can be used to select single cells based on a variety of phenotypes, such as cell surface markers, cell proliferation and drug response. We then show that a common genomic procedure, RNA-seq, can be readily adapted to the single cells isolated from these rafts. We show that data generated using microrafts and our modified RNA-seq protocol compared favorably with the Fluidigm C1. We then used microraft arrays to select pancreatic cancer cells that proliferate in spite of cytotoxic drug treatment. Our single cell RNA-seq data identified several expected and novel gene expression changes associated with early drug resistance.

This RDS file contains processed single-cell RNA sequencing (scRNA-seq) data comparing immune cell populations from germ-free (GF) and specific-pathogen-free (SPF) mice. The dataset includes:Samples: Peripheral blood (PB) and bone marrow (BM) from GF and SPF miceCell Counts:Raw: 21,827 cells (PB) and 19,940 cells (BM)Quality-filtered: 18,344 high-quality cells (PB) and 16,537 high-quality cells (BM)Gene Coverage: Median 1,426 genes per cell (PB) and 1,391 genes per cell (BM)Cell Classifications: 18 major cell identities further divided into 25 subpopulationsAnnotation: Cells identified using established marker genes for blood cells

We performed CODEX (co-detection by indexing) multiplexed imaging on 64 sections of the human intestine (~16 mm2) from 8 donors (B004, B005, B006, B008, B009, B010, B011, and B012) using a panel of 57 oligonucleotide-barcoded antibodies. Subsequently, images underwent standard CODEX image processing (tile stitching, drift compensation, cycle concatenation, background subtraction, deconvolution, and determination of best focal plane), single cell segmentation, and column marker z-normalization by tissue. The outputs of this process were data frames of 2.6 million cells with 57 antibody fluorescence values quantified from each marker. Each cell has its cell type, cellular neighborhood, community of neighborhooods, and tissue unit defined with x, y coordinates representing pixel location in the original image. This is from a total of 25 cell types, 20 multicellular neighborhoods, 10 communities of neighborhoods, and 3 tissue segments that could be used to understand the cellular interactions, composition, and structure of the human intestine from the duodenum to the sigmoid colon and understand differences between different areas of the intestine. This data could be used as a healthy baseline to compare other single-cell datasets of the human intestine, particularly multiplexed imaging ones. The overall structure of the datasets is individual cells segmented out in each row. Columns MUC2 through CD161 are the markers used for clustering the cell types. These are the columns that are the values of the antibody staining the target protein within the tissue quantified at the single-cell level. This value is the per cell/area averaged fluorescent intensity that has subsequently been z normalized along each column as described above. OLFM4 through MUC6 were captured in the quantification but not used within the clustering of cell types. Other columns are explained in the table in the Usage Notes section below. Along with this main data table, there is also a donor metadata table that links the donor ids to clinical metadata such as: age, sex, race, BMI, history of diabetes, history of cancer, history of hypertension, and history of gastorintestinal disease. The raw imaging data can be found at (https://portal.hubmapconsortium.org/). We have created a landing page with links to all the raw dataset IDs and the HuBMAP ID for this Collection is HBM692.JRZB.356 and the DOI is:10.35079/HBM692.JRZB.356. This can be used to also pair it with the matched snRNAseq and snATACseq for each section of tissue. Methods For a detailed description of each of the steps of protocols and processes to obtain this data see the detailed materials and methods in the associated manuscript. Briefly, intestine pieces from 8 different sites across the small intestine and colon were taken and frozen in OCT. These were assembled into an array of 4 tissues, cut into 7 um slices, and stained with a panel of 54 CODEX DNA-oligonucleotide barcoded antibodies. Tissues were imaged with a Keyence microscope at 20x objective and then processed using image stitching, drift compensation, deconvolution, and cycle concatenation. Processed data were then segmented using CellVisionSegmenter, a neural network R-CNN-based single-cell segmentation algorithm. Cell type analysis was completed on B004, 5, and 6 by z normalization of protein markers used for clustering and then overclustered using leiden-based clustering. The cell type labels were verified by looking back at the original image. Cell type labels were transferred to other donors using STELLAR framework for annotating spatially resolved single-cell data as described in detail in the companion Nature Methods manuscript. With set cell type labels we performed neighborhood analysis by clustering windows of the 10 nearest neighbors around a given cell and were named based on cell type enrichment and location in the tissue. Similarly, communities of neighborhoods were determined by taking the 100 nearest neighbors with the neighborhood labels and clustered. Finally, tissue segments were determined through multiple rounds of clustering the 300 nearest neighbors of the community labels of each cell. Broad categories for cell types, neighborhoods, and communities were expert annotated based on epithelial, immune, or other stromal compartments.

allGOnKEGG_June2020.gmt: Tab delimited file describing gene sets from KEGG and GO pathways

Control1_Deprez_exprMatrix.tsv: COntrol dataset 1

Contro2_Madissoon_lung.cellxgene.h5ad: Control dataset 2

COVID19COMORBID_DATABASE_UDDINLAB.xls: Comorbid gene dataset. The genes which are upregulated in severe COVID19-comorbid conditions

Final_LUNGMarker_DB_UDDINLAB.xlsx: Database consisting of canonical markers for cells associated to human lung region constructed by combining cell marker databases and thorough literature

All files starting with GSM: Matrix(h5 file) and annotation file for bronchoalveolar lavage fluid (BALF) from 6 severe and 3 moderate COVID-19 patients and 3 healthy control (From Liao et al., 2020: PMID: 32398875)

Human_cell_markers.txt: Markers from CellMarker database

PanglaoDB_markers_27_Mar_2020.tsv: Markers from PanglaoDB database

All files starting with sample: Sample Seurat processed file (sample.rds); Sample file for executing gene enrichment analysis; Sample file for lung marker database

All files starting with Validation: Validation files used in the paper.

A comprehensive cellular anatomy of normal human bladder is vital to address the cellular origins of benign bladder disease and bladder cancer. The physiological function and pathological changes of bladder are associated with its cell type. To investigate the classification and underlying function of bladder cells, we conducted single-cell RNA sequencing (scRNA-seq) of 12,423 cells from human bladder and 12,884 cells from mouse bladder.Here, we show a single-cell transcriptomic map of the human and mouse bladders, including 16 clusters of human bladder cells and 15 clusters of mouse bladder cells. Homology and heterogeneity of human and mouse were compared and analyzed by our study, and indicated that there were both conservative and heterogeneous aspects of human and mouse evolution. We also discovered two novel types of human bladder cells; one type is ADRA2A+ and HRH2+ interstitial cells which may be associated with nerve conduction and allergic reactions; the other type is TNNT1+ epithelial cells that may be different from other epithelial cells in biological function and verify in rat and mouse tissues. Collectively, this analysis provides a resource for the discovery of novel cell type, specific cell markers, signaling receptors, and genes dataset that will help us to study the relationship between cell types and diseases. In this website, we show the average expression of every gene in every cell cluster of human and mouse bladder. At the same time, we shared differentially expressed genes (DEGs) of every cell cluster in human or mouse bladder.

Innate lymphoid cells (ILCs) are enriched at mucosal surfaces where they respond rapidly to environmental stimuli and contribute to both tissue inflammation and healing. To gain insight into the role of ILCs in the pathology and recovery from COVID-19 infection, we employed a multi-omic approach consisting of Abseq and targeted mRNA sequencing to respectively probe the surface marker expression, transcriptional profile and heterogeneity of ILCs in peripheral blood of patients with COVID-19 compared with healthy controls. We found that the frequency of ILC1 and ILC2 cells was significantly increased in COVID-19 patients. Moreover, all ILC subsets displayed a significantly higher frequency of CD69-expressing cells, indicating a heightened state of activation. ILC2s from COVID-19 patients had the highest number of significantly differentially expressed (DE) genes. The most notable genes DE in COVID-19 vs healthy participants included a) genes associated with responses to virus infections and b) genes that support ILC self-proliferation, activation and homeostasis. In addition, differential gene regulatory network analysis revealed ILC-specific regulons and their interactions driving the differential gene expression in each ILC. Overall, this study provides mechanistic insights into the characteristics of ILC subsets activated during COVID-19 infection. Methods Study participants, blood draws and processing Participants were recruited as described previously from adults who had a positive SARS-COV-2 RT-PCR test at Stanford Health Care (NCT04373148). Collection of Covid samples occurred between May to December 2020. The cohort used in this study consisted of asymptomatic (n=2), mild (n=17), and moderate (n=3) COVID-19 infections, some of whom developed long term COVID-19 (n=15). The clinical case severities at the time of diagnosis were defined as asymptomatic, moderate or mild according to the guidelines released by NIH. Long term (LT) COVID was defined as symptoms occurring 30 or more days after infection, consistent with CDC guidelines. Some participants in our study continued to have LT COVID symptoms 90 days after diagnosis (n=12). Exclusion criteria for COVID sample study were NIH severity diagnosis of severe or critical at the time of positive covid test. Samples selected for this study were obtained within 76 days of positive PCR COVID-19 test date. Healthy controls were selected who had sample collection before 2020. Informed consent was obtained from all participants. All protocols were approved by the Stanford Administrative Panel on Human Subjects in Medical Research. Peripheral blood was drawn by venipuncture and using validated and published procedures, peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll-based density gradient centrifugation, frozen in aliquots and stored in liquid nitrogen at -80°C , until thawing. A summary of participant demographics is presented in Supp. Table 1.
ILC Enrichment, single cell captures for Abseq and targeted mRNAseq Participant PBMCs were thawed, and each sample stained with Sample Tag (BD #633781) at room temperature for 20 minutes. Samples were combined in healthy control or COVID-19 tubes. Cells were surface stained with a panel of fluorochrome-conjugated antibodies (Supp. Table 2) in buffer (PBS with 0.25% BSA and 1mM EDTA) for 20 minutes at room temperature prior to immunomagnetic negative selection for ILCs. Following ILC enrichment using the EasySep human Pan-ILC enrichment kit (StemCell Technologies #17975), cells from healthy and COVID-19 recovered participants were counted and normalized before combining. ILCs were sorted using a BD FACS Aria at the Stanford FACS facility prior to incubation with AbSeq oligo-linked mAbs (Supp. Table 3). Sorted cells were processed by the Stanford Human Immune Monitoring Center (HIMC) using the BD Rhapsody platform. Library was prepared using the BD Immune Response Targeting Panel (BD Kit #633750) with addition of custom gene panel reagents (Supp. Table 4) and sequenced on Illumina NovaSeq 6000 at Stanford Genomics Sequencing Center (SGSC). ILCs were identified as Lineageneg (CD3neg, CD14neg, CD34neg, CD19neg), NKG2Aneg, CD45+ and ILCs further defined as CD127+CD161+ and as subsets: ILC1 (CD117negCRTH2neg), ILC2 (CRTH2+) and ILCp (CD117+CRTH2neg) (Supp. Fig. 1). Computational data analysis The above multi-modal setup allowed paired measurements of cellular transcriptome and cell surface protein abundance. The ILC1, ILC2 and ILCp cells were manually gated based on the abundance profile of CD127, CD117, CD161 and CRTH2 (Supp. Fig. 1). Before the integrative analysis, the complete multi-modal single cell dataset containing ILC subsets was converted into single Seurat object. All the subsequent protein-level and gene-level analyses were performed using multimodal data analysis pipeline of Seurat R package version 4.0. The normalized and scaled protein abundance profile was used for estimating the integrated harmony dimensions using runHarmony function in Seurat R package (reduction= ‘apca’ and group.by.vars = ‘batch’) . The batch corrected harmony embeddings were then used for computing the Uniform Manifold Approximation and Projection (UMAP) dimensions to visualize the clusters of ILC subsets. Differential marker analysis of surface proteins, between two groups of cells (COVID-19 and Healthy cohort), from abseq panels was computed with normalized and scaled expression values using FindMarkers function from Seurat R package (test.use=’wilcox’). Similarly, differential gene expression was performed on normalized and scaled gene expression values from between two groups of cells (COVID-19 and Healthy cohort) using the FindMarkers function from Seurat R package (test.use=’MAST’ and latent.vars=’batch’). Genes with log-fold change > 0.5 and adjusted p-value < 0.05 (method: Benjamini-Hochberg) (were considered as significant for further evaluation. The resulting adjusted p-values box-plots were plotted using ggplot2 R package (version 3.4.2) after computing the number of cells expressing a given protein or gene in each sample. Pathway enrichment analysis of DE genes was performed using web-server metascape (version 3.5). The AUCells score and gene regulatory network analysis was performed using pySCENIC pipeline (version 0.12.1). Gene regulatory network was reconstructed using GRNBoost2 algorithm and the list of TFs in humans (genome version: hg38) were obtained from cisTarget database. (https://resources.aertslab.org/cistarget). Cellular enrichment (aka AUCell) analysis that measures the activity of TF or gene signatures across all single cells was performed using aucell function in pySCENIC python library. The ggplot2 R package (version 3.4.2) was used for boxplot visualization. The differential gene co-expression analysis was performed using scSFMnet R package. Circular plots were generated using the R package circlize (version 0.4.15).

A Transcriptome Database for Astrocytes, Neurons, and Oligodendrocytes: A New Resource for Understanding Brain Development and Function Understanding the cell-cell interactions that control CNS development and function has long been limited by the lack of methods to cleanly separate astrocytes, neurons, and oligodendrocytes. Here we describe the first method for the isolation and purification of developing and mature astrocytes from mouse forebrain. This method takes advantage of the expression of S100β by astrocytes. We used fluorescent activated cell sorting (FACS) to isolate EGFP positive cells from transgenic mice that express EGFP under the control of an S100β promoter. By depletion of astrocytes and oligodendrocytes we obtained purified populations of neurons, while by panning with oligodendrocyte-specific antibodies we obtained purified populations of oligodendrocytes. Using GeneChip Arrays we then created a transcriptome database of the expression levels of over 20,000 genes by gene profiling these three main CNS neural cell types at postnatal ages day 1 to 30. This database provides the first global characterization of the genes expressed by mammalian astrocytes in vivo and is the first direct comparison between the astrocyte, neuron, and oligodendrocyte transcriptomes. We demonstrate that Aldh1L1, a highly expressed astrocyte gene, is a highly specific antigenic marker for astrocytes with a substantially broader, and therefore potentially more useful, pattern of astrocyte expression than the traditional astrocyte marker GFAP. This transcriptome database of acutely isolated and highly pure populations of astrocytes, neurons and oligodendrocytes provides a resource to the neuroscience community by providing improved cell type specific markers and for better understanding of neural development, function, and disease. We acutely purified mouse astrocytes from early postnatal ages (P1) to later postnatal ages (P30), when astrocyte differentiation is morphologically complete (Bushong et al., 2004), and acutely purified mouse OL-lineage cells from stages ranging from OPCs to newly differentiated OLs to myelinating OLs. We extracted RNA from each of these highly purified, acutely isolated cell types and used GeneChip Arrays to determine the expression levels of over 20,000 genes and construct a comprehensive database of cell type specific gene expression in the mouse forebrain. Analysis of this database confirms cell type specific expression of many well characterized and functionally important genes. In addition, we have identified thousands of new cell type enriched genes, thereby providing important new information about astrocyte, OL, and neuron interactions, metabolism, development, and function. This database provides a comparison of the genome-wide transcriptional profiles of the main CNS cell types and is a resource to the neuroscience community for better understanding the development, physiology, and pathology of the CNS. Keywords: Developmental CNS Cell type comparision FACS purification of astrocytes: Dissociated forebrains from S100β-EGFP mice were resuspended in panning buffer (DBPS containing 0.02% BSA and 12.5 U/ml DNase) and sequentially incubated on the following panning plates: secondary antibody only plate to deplete microglia, O4 plate to deplete OLs, PDGFRα plate to deplete OPCs, and a second O4 plate to deplete any remaining OLs. This procedure was sufficient to deplete all OL-lineage cells from animals P8 and younger, however, in older animals that had begun to myelinate, additional depletion of OLs and myelin debris was accomplished as follows. The nonadherent cells from the last O4 dish were harvested by centrifugation, and the cells were resuspended in panning buffer containing GalC, MOG, and O1 supernatant and incubated for 15 minutes at room temperature. The cell suspension was washed and then resuspended in panning buffer containing 20 μg donkey anti-mouse APC for 15 minutes. The cells were washed and resuspended in panning buffer containing propidium iodide (PI). EGFP+ astrocytes were then purified by fluorescence activated cell sorting (FACS). Dead cells were gated out using high PI staining and forward light scatter. Astrocytes were identified based on high EGFP fluorescence and negative APC fluorescence from indirect immunostaining for OL markers GalC, MOG, and O1. Cells were sorted twice and routinely yielded >99.5% purity based on reanalysis of double sorted cells.; FACS purification of neurons: EGFP- cells were the remaining forebrain cells after microglia, OLs, and astrocytes had been removed, and were primarily composed of neurons, and to a lesser extent, endothelial cells (we estimate < 4% endothelial cells at P7 and < 20% endothelial cells at P17). EGFP- cells from S100β-EGFP dissociated forebrain were FACS purified in parallel with astrocyte purification and were sorted based on their negative EGFP fluorescence immunofluorescence. Cells were sorted twice and routinely yielded >99.9% purity. In independent preparations, the EGFP- cell population was additionally depleted of endothelial cells and pericytes by sequentially labeling with biotin-BSL1 lectin and streptavidin-APC while also labeling for OL markers as described above. Cells were sorted twice and routinely yielded >99.9% purity.; Panning purification of oligodendrocyte lineage cells: Dissociated mouse forebrains were resuspended in panning buffer. In order to deplete microglia, the single-cell suspension was sequentially panned on four BSL1 panning plates. The cell suspension was then sequentially incubated on two PDGFRα plates (to purify and deplete OPCs), one A2B5 plate (to deplete any remaining OPCs), two MOG plates (to purify and deplete myelinating OLs), and one GalC plate (to purify the remaining PDGFRα-, MOG-, OLs). The adherent cells on the first PDGFRα, MOG, and GalC plates were washed to remove all antigen-negative nonadherent cells. The cells were then lysed while still attached to the panning plate with Qiagen RLT lysis buffer, and total RNA was purified. Purified OPCs were >95% NG2 positive and 0% MOG positive. Purified Myelin OLs were 100% MOG positive, >95% MBP positive, and 0% NG2 positive. Purified GalC OLs depleted of OPCs and Myelin OLs were <10% MOG positive and ~50% weakly NG2 positive, a reflection of their recent development as early OLs.; Data normalization and analysis: Raw image files were processed using Affymetrix GCOS and the MAS 5.0 algorithm. Intensity data was normalized per chip to a target intensity TGT value of 500, and expression data and absent/present calls for individual probe sets were determined. Gene expression values were normalized and modeled across arrays using the dChip software package with invariant-set normalization and a PM model. (www.dchip.org, Li and Wong, 2001). The 29 samples were grouped into 9 sample types: Astros P7-P8, Astros P17, Astros P17-gray matter (P17g), Neurons P7, Neurons P17, Neurons-endothelial cell depleted (P7n, P17n), OPCs, GalC-OLs, and MOG-OLs. Gene filtering was performed to select probe sets that were consistently expressed in at least one cell type, where consistently expressed was defined as being called present and having a MAS 5.0 intensity level greater than 200 in at least two-thirds of the samples in the cell type. We identified 20,932 of the 45,037 probe sets that were consistently expressed in at least one of the nine cell types. The Significance Analysis of Microarrays (SAM) method (Tusher et al., 2001) was used to determine genes that were significantly differentially expressed between different cell types (see Supplemental Table S2 for SAM cell type groupings). Clustering was performed using the hclust method with complete linkage in R. Expression values were transformed for clustering by computing a mean expression value for the gene using those samples in the corresponding SAM statistical analysis, and then subtracting the mean from expression intensities. In order to preserve the log2 scale of the data, unless otherwise indicated, no normalization by variance was performed. Plots were created using the gplots package in R. The Bioconductor software package (Gentleman et al., 2004) was used throughout the expression analyses. Functional analyses were performed through the use of Ingenuity Pathways Analysis (Ingenuity® Systems, www.ingenuity.com).

bioRxiv preprint: https://www.biorxiv.org/content/10.1101/2020.04.02.022566v4

Contact: Grigorios Georgolopoulos (ggeorgol@altius.org); Jeff Vierstra (jvierstra@altius.org)

Lineage commitment and differentiation is driven by the concerted action of master transcriptional regulators at their target chromatin sites. Multiple efforts have characterized the key transcription factors (TFs) that determine the various hematopoietic lineages. However, the temporal interactions between individual TFs and their chromatin targets during differentiation and how these interactions dictate lineage commitment remains poorly understood. Here we delineate the temporal interplay between the cis- and the trans-regulatory landscape in establishing lineage commitment and differentiation in human hematopoiesis by performing a dense timecourse of chromatin accessibility (DNase I-seq), and gene expression (total and single cell RNA-seq).

All data uploaded correspond to human genome build version GRCh38.

Contents

DNase I Hotspot (DHS) metadata: Supplementary_Data_1.txt

DNase I Hotspot quantile-normalized counts: A tab-separated matrix with quantile-normalized DNase I density counts from 79,085 FDR 5% hotspots, across 12 erythroid differentiation timepoints from 3 donors, present in at least n=2 samples. Rows correspond to DHS information in Supplementary_Data_1.txt (hotspots.fdr.0.05.qnorm.counts.tsv.gz)

Column information for DNase I Hotspot quantile-normalized counts: hotspots.fdr.0.05.qnorm.counts.info.tsv

Developmentally regulated gene metadata (erythroid): Supplementary_Data_2.csv

Gene matrix of quantile-normalized FPKM values (erythroid): A tab-separated matrix with the quantile-normalized FPKM values of all detected genes, across 13 erythroid differentiation timepoints from 3 donors. (fpkm_erythroid_qnorm.tsv.gz)

Column information for the quantile-normalized FPKM gene matrix (erythroid): A tab-separated table (fpkm_erythroid_qnorm.info.tsv)

CD34+ HSPC TADs at 10kb resolution: Supplementary_Data_3.bed

Day 11 ex vivo erythroid progenitor TADs at 10kb resolution: Supplementary_Data_4.bed

Transcription factor motif enrichment per DHS cluster: Supplementary_Data_5.csv

Correlation information (links) between developmentally regulated DHS and target genes: Supplementary_Data_6.csv

Chromatin anchor loops called from 10kb resolution Hi-C data: Supplementary_Data_7.bedgraph

Developmentally regulated gene metadata (megakaryocytic): Supplementary_Data_8.csv

Gene matrix of quantile-normalized FPKM values (megakaryocytic): A tab-separated matrix with the quantile-normalized FPKM values of all detected genes, across 13 megakaryocytic differentiation timepoints from 3 donors. (fpkm_megakaryocyte_qnorm.tsv.gz)

Column information for the quantile-normalized FPKM gene matrix (megakaryocytic): A tab-separated table (fpkm_megakaryocyte_qnorm.info.tsv)

Marker (differentially expressed) genes per single cell population: Supplementary_Data_9.csv

A SCANPY h5ad Annotated DataFrame object: Annotated Data frame anndata in h5ad format including the gene-by-cell count matrix, Velocyto splicing kinetics (RNA velocity) information layer, along with obs, obsm, var, varm, and uns layers. (SCANPY_anndata_object.h5ad)

Adult zebrafish possess a strong regenerative capacity. Focusing on the heterogeneity of zebrafish optic nerve regeneration, single-cell sequencing analysis of the optic nerve of adult zebrafish two weeks after injury identified five major cell types: fibroblasts, vascular wall cells, immune cells, mature oligodendrocytes, and oligodendrocytes in the process of myelination.Experimental ProcedureDissociation of optic nerve to a single-cell suspensionOptic nerve samples from n=18 zebrafish (half of male and female) at 2 wpi were pooled to create a single-cell suspension, ensuring sufficient cell capture while minimizing inter-individual variability. Therefore, there are no biological replicates due to small tissue size of the optic nerve in zebrafish and low expression of mRNA in the optic nerve (Yu et al., 2024). Dissection of the optic nerve and then were minced with a sterile scalpel into 1 mm fragments, suspended in 5ml of digestion buffer consisting of 2 mg/mL Collagenase type II and 200U/ml DNase I in RPMI medium, and incubated in 37℃ water bath with shaking for 30 min. The suspension was passed through a 100 μm filter and centrifugated (400g, 10 min, 4℃). Pelleted cells were resuspended in red blood cell lysis buffer, incubated for 2 min, passed through a 40 μm filter, collected by centrifugation (400g, 10 min, 4℃) and resuspended in PBS containing 0.04% BSA. Cells were manually counted by Trypan blue and AO-PI (LUNA, D23001) after each centrifugation (400g, 10 min, 4℃) and resuspended. Single cells were processed using Chromium Controller (10X Genomics) according to the manufacturer’s protocol.Single-cell RNA sequencingBy using Chromium Next GEM Single Cell 3ʹ Kit v3.1and Chromium Next GEM Chip G Single Cell Kit, we performed single cell 3’gene expression profiling. The cell suspension was loaded onto the Chromium single cell controller (10x Genomics) to generate single-cell gel beads in the emulsion according to the manufacturer’s protocol. Captured cells were lysed and the released RNA were barcoded through reverse transcription in individual GEMs. Cell-barcoded 3’gene expression libraries were sequenced on an Illumina NovaSeq6000 system (Illumina, USA) by Shanghai Biochip Co., Ltd.,China.Single-cell RNA sequencing and data analysisThe raw single-cell RNA sequencing (scRNA-Seq) data were processed as described in previous paper ((Yu et al., 2024)). Low-quality cells were excluded based on the retaining criteria as RNA feature count between 300 and 5000, RNA count between 500-8000, a mitochondrial gene percentage below 20%, hemoglobin- and red blood cell–related gene percentages below 10%, and total RNA counts not exceeding the 95th percentile. A total of 3,359 cells were sequenced, and 1,341 cells were retrieved after sequencing and quality control. Clusters that did not belong to optic nerve tissue were excluded from further analysis. The detailed quality control workflow is present in Supplementary file1. Optic nerve samples were integrated using Seurat v5.1, followed by data normalizing, scaling, dimensional reduction, and clustering. Visualization was performed with t-SNE using the first 30 principal components. Cell types were annotated based on canonical marker genes, including mature oligodendrocyte (tspan2a, mag, mbpa), myelin forming oligodendrocyte (egr2b, gldn, mbpb, s100b, si:ch211-234p6.5, ndrg1a), fibroblast (col1a1a, col1a1b, dcn), immune cell (nr4a1, srgn) and mural cell (rgs5a, rbpms2a, myh11a). We note that this immune cell cluster does not express any microglial or macrophage markers but instead shows high expression of more generalized immune function–related genes. Consequently, we can only designate it as an “immune cell” cluster and cannot further subdivide it. Gene Ontology (GO) enrichment analysis was conducted to identify functional pathways and biological processes associated with each group. Proliferation scores were computed based on curated gene sets representing proliferation-related pathways

ObjectiveThis study amied to investigate the prognostic characteristics of triple negative breast cancer (TNBC) patients by analyzing B cell marker genes based on single-cell and bulk RNA sequencing.MethodsUtilizing single-cell sequencing data from TNBC patients, we examined tumor-associated B cell marker genes. Transcriptomic data from The Cancer Genome Atlas (TCGA) database were used as the foundation for predictive modeling. Independent validation set was conducted using the GSE58812 dataset. Immune cell infiltration into the tumor was assessed through various, including XCELL, TIMER, QUANTISEQ, CIBERSORT, CIBERSORT-ABS, and ssGSEA. The TIDE score was utilized to predict immunotherapy outcomes. Additional investigations were conducted on the immune checkpoint blockade gene, tumor mutational load, and the GSEA enrichment analysis.ResultsOur analysis encompassed 22,106 cells and 20,556 genes in cancerous tissue samples from four TNBC patients, resulting in the identification of 116 B cell marker genes. A B cell marker gene score (BCMG score) involving nine B cell marker genes (ZBP1, SEL1L3, CCND2, TNFRSF13C, HSPA6, PLPP5, CXCR4, GZMB, and CCDC50) was developed using TCGA transcriptomic data, revealing statistically significant differences in survival analysis (P