DOI: 10.5061/dryad.t4b8gtj34
title: Large-scale integration of single-cell transcriptomic data captures transitional progenitor states in mouse skeletal muscle regeneration

description: Skeletal muscle repair is driven by the coordinated self-renewal and
 fusion of myogenic stem and progenitor cells. Single-cell gene expression
 analyses of myogenesis have been hampered by the poor sampling of rare and
 transient cell states that are critical for muscle repair, and do not
 inform the spatial context that is important for myogenic differentiation.
 Here, we demonstrate how large-scale integration of single-cell and
 spatial transcriptomic data can overcome these limitations. We created a
 single-cell transcriptomic dataset of mouse skeletal muscle by
 integration, consensus annotation, and analysis of 23 newly collected
 scRNAseq datasets and 88 publicly available single-cell (scRNAseq) and
 single-nucleus (snRNAseq) RNA-sequencing datasets. The resulting dataset
 includes more than 365,000 cells and spans a wide range of ages, injury,
 and repair conditions. Together, these data enabled identification of the
 predominant cell types in skeletal muscle, and resolved cell subtypes,
 including endothelial subtypes distinguished by vessel-type of origin,
 fibro/adipogenic progenitors defined by functional roles, and many
 distinct immune populations. The representation of different experimental
 conditions and the depth of transcriptome coverage enabled robust
 profiling of sparsely expressed genes. We built a densely sampled
 transcriptomic model of myogenesis, from stem cell quiescence to myofiber
 maturation and identified rare, transitional states of progenitor
 commitment and fusion that are poorly represented in individual datasets.
 We performed spatial RNA sequencing of mouse muscle at three time points
 after injury and used the integrated dataset as a reference to achieve a
 high-resolution, local deconvolution of cell subtypes. We also used the
 integrated dataset to explore ligand-receptor co-expression patterns and
 identify dynamic cell-cell interactions in muscle injury response. We
 provide a public web tool to enable interactive exploration and
 visualization of the data. Our work supports the utility of large-scale
 integration of single-cell transcriptomic data as a tool for biological
 discovery.
description: <strong>Mice.</strong> The Cornell University
 Institutional Animal Care and Use Committee (IACUC) approved all animal
 protocols, and experiments were performed in compliance with its
 institutional guidelines. Adult C57BL/6J mice (mus musculus) were obtained
 from Jackson Laboratories (#000664; Bar Harbor, ME) and were used at 4-7
 months of age. Aged C57BL/6J mice were obtained from the National
 Institute of Aging (NIA) Rodent Aging Colony and were used at 20 months of
 age. For new scRNAseq experiments, female mice were used in each
 experiment. <strong>Mouse injuries and
 single-cell isolation.</strong> To induce muscle injury, both
 tibialis anterior (TA) muscles of old (20 months) C57BL/6J mice were
 injected with 10 µl of notexin (10 µg/ml; Latoxan; France). At 0, 1, 2,
 3.5, 5, or 7 days post-injury (dpi), mice were sacrificed and TA muscles
 were collected and processed independently to generate single-cell
 suspensions. Muscles were digested with 8 mg/ml Collagenase D (Roche;
 Switzerland) and 10 U/ml Dispase II (Roche; Switzerland), followed by
 manual dissociation to generate cell suspensions. Cell suspensions were
 sequentially filtered through 100 and 40 μm filters (Corning Cellgro
 #431752 and #431750) to remove debris. Erythrocytes were removed through
 incubation in erythrocyte lysis buffer (IBI Scientific
 #89135-030). <strong>Single-cell RNA-sequencing
 library preparation.</strong> After digestion, single-cell
 suspensions were washed and resuspended in 0.04% BSA in PBS at a
 concentration of 106 cells/ml. Cells were counted manually with a
 hemocytometer to determine their concentration. Single-cell RNA-sequencing
 libraries were prepared using the Chromium Single Cell 3’ reagent kit v3
 (10x Genomics, PN-1000075; Pleasanton, CA) following the manufacturer’s
 protocol. Cells were diluted into the Chromium Single Cell A Chip to yield
 a recovery of 6,000 single-cell transcriptomes. After preparation,
 libraries were sequenced using on a NextSeq 500 (Illumina; San Diego, CA)
 using 75 cycle high output kits (Index 1 = 8, Read 1 = 26, and Read
 2 = 58). Details on estimated sequencing saturation and the number of
 reads per sample are shown in Sup. Data 1.
 <strong>Spatial RNA sequencing library
 preparation.</strong> Tibialis anterior muscles of adult (5 mo)
 C57BL6/J mice were injected with 10µl notexin (10 µg/ml) at 2, 5, and 7
 days prior to collection. Upon collection, tibialis anterior muscles were
 isolated, embedded in OCT, and frozen fresh in liquid nitrogen. Spatially
 tagged cDNA libraries were built using the Visium Spatial Gene Expression
 3’ Library Construction v1 Kit (10x Genomics, PN-1000187; Pleasanton, CA)
 (Fig. S7). Optimal tissue permeabilization time for 10 µm thick sections
 was found to be 15 minutes using the 10x Genomics Visium Tissue
 Optimization Kit (PN-1000193). H&E stained tissue sections were
 imaged using Zeiss PALM MicroBeam laser capture microdissection system and
 the images were stitched and processed using Fiji ImageJ software. cDNA
 libraries were sequenced on an Illumina NextSeq 500 using 150 cycle high
 output kits (Read 1=28bp, Read 2=120bp, Index 1=10bp, and Index 2=10bp).
 Frames around the capture area on the Visium slide were aligned manually
 and spots covering the tissue were selected using Loop Browser v4.0.0
 software (10x Genomics). Sequencing data was then aligned to the mouse
 reference genome (mm10) using the spaceranger v1.0.0 pipeline to generate
 a feature-by-spot-barcode expression matrix (10x Genomics).
 <strong>Download and alignment of single-cell RNA
 sequencing data.</strong> For all samples available via SRA,
 parallel-fastq-dump (github.com/rvalieris/parallel-fastq-dump) was used to
 download raw .fastq files. Samples which were only available as .bam files
 were converted to .fastq format using bamtofastq from 10x Genomics (github.com/10XGenomics/bamtofastq). Raw reads were aligned to the mm10 reference using cellranger (v3.1.0). <strong>Preprocessing and batch correction of single-cell RNA sequencing datasets. </strong>First, ambient RNA signal was removed using the default SoupX (v1.4.5) workflow (autoEstCounts and adjustCounts; github.com/constantAmateur/SoupX). Samples were then preprocessed using the standard Seurat (v3.2.1) workflow (NormalizeData, ScaleData, FindVariableFeatures, RunPCA, FindNeighbors, FindClusters, and RunUMAP; github.com/satijalab/seurat). Cells with fewer than 750 features, fewer than 1000 transcripts, or more than 30% of unique transcripts derived from mitochondrial genes were removed. After preprocessing, DoubletFinder (v2.0) was used to identify putative doublets in each dataset, individually. BCmvn optimization was used for PK parameterization. Estimated doublet rates were computed by fitting the total number of cells after quality filtering to a linear regression of the expected doublet rates published in the 10x Chromium handbook. Estimated homotypic doublet rates were also accounted for using the modelHomotypic function. The default PN value (0.25) was used. Putative doublets were then removed from each individual dataset. After preprocessing and quality filtering, we merged the datasets and performed batch-correction with three tools, independently- Harmony (github.com/immunogenomics/harmony) (v1.0), Scanorama (github.com/brianhie/scanorama) (v1.3), and BBKNN (github.com/Teichlab/bbknn) (v1.3.12). We then used Seurat to process the integrated data. After initial integration, we removed the noisy cluster and re-integrated the data using each of the three batch-correction tools. <strong>Cell type annotation.</strong> Cell types were determined for each integration method independently. For Harmony and Scanorama, dimensions accounting for 95% of the total variance were used to generate SNN graphs (Seurat::FindNeighbors). Louvain clustering was then performed on the output graphs (including the corrected graph output by BBKNN) using Seurat::FindClusters. A clustering resolution of 1.2 was used for Harmony (25 initial clusters), BBKNN (28 initial clusters), and Scanorama (38 initial clusters). Cell types were determined based on expression of canonical genes (Fig. S3). Clusters which had similar canonical marker gene expression patterns were merged. <strong>Pseudotime workflow.</strong> Cells were subset based on the consensus cell types between all three integration methods. Harmony embedding values from the dimensions accounting for 95% of the total variance were used for further dimensional reduction with PHATE, using phateR (v1.0.4) (github.com/KrishnaswamyLab/phateR). <strong>Deconvolution of spatial RNA sequencing spots.</strong> Spot deconvolution was performed using the deconvolution module in BayesPrism (previously known as “Tumor microEnvironment Deconvolution”, TED, v1.0; github.com/Danko-Lab/TED). First, myogenic cells were re-labeled, according to binning along the first PHATE dimension, as “Quiescent MuSCs” (bins 4-5), “Activated MuSCs” (bins 6-7), “Committed Myoblasts” (bins 8-10), and “Fusing Myoctes” (bins 11-18). Culture-associated muscle stem cells were ignored and myonuclei labels were retained as “Myonuclei (Type IIb)” and “Myonuclei (Type IIx)”. Next, highly and differentially expressed genes across the 25 groups of cells were identified with differential gene expression analysis using Seurat (FindAllMarkers, using Wilcoxon Rank Sum Test; results in Sup. Data 2). The resulting genes were filtered based on average log2-fold change (avg_logFC &gt; 1) and the percentage of cells within the cluster which express each gene (pct.expressed &gt; 0.5), yielding 1,069 genes. Mitochondrial and ribosomal protein genes were also removed from this list, in line with recommendations in the BayesPrism vignette. For each of the cell types, mean raw counts were calculated across the 1,069 genes to generate a gene expression profile for BayesPrism. Raw counts for each spot were then passed to the run.Ted function, using the “GEP” option for input.type and default parameters for the remaining inputs. Final Gibbs theta values were used as estimates for the fraction of transcripts from each spot that were derived from each of the 25 cell types. <strong>Spot co-occurrence of cell subtypes.</strong> Cell type co-occurrence was computed as before by Mantri et al. Briefly, for each spot cell types were ordered according to the computed theta values (estimated percentage of reads attributed to a cell type). Up to 10 cell types with the highest theta values were tallied for each spot. Cell types with theta values lower than 0.01 (1% of transcripts) were not counted. <strong>Ligand-receptor analysis with CellChat.</strong> Ligand-receptor analysis and visualization was performed using CellChat v1.1.0 (github.com/sqjin/CellChat). The cell type labels used were derived from the Harmony integration results using all single-cell and single-nucleus data sources. Myogenic cell types were derived from binning along the PHATE embedding (Fig. 2b). Default values were used for the parameterization of each step.
description: Additional information as well as the code used to prepare these
 data can be found at github.com/mckellardw/scMuscle. These data comprise scMuscle v1.1 Seurat/CellChat Objects (see our github repository for other version tracking information).