A single-cell atlas and Shiny-based framework for murine lung injury and remodeling

Single-cell RNA sequencing has provided unprecedented insights into the cellular heterogeneity of acute and chronic pulmonary injury. Lack of data-driven analytical frameworks remains a significant barrier to reproducible and biologically robust analyses. To address this challenge, we developed a modular, interactive analysis platform for single-cell RNA sequencing data processing and visualization. As a case study, we applied the framework to an immune-enriched murine lung dataset containing samples from different batches and across four conditions: control, inflammation and remodeling, fibrosis, and fibrotic lung with ozone exposure. Manual annotations performed on individual samples were used as ground-truth labels to evaluate automated cell-type annotations generated using the Tabula Muris Senis and LungMAP reference atlases. These annotations were also utilized to benchmark the performance of data preprocessing and integration combinations on batch-effect removal and the preservation of biological identity. Within the LIANA framework, we evaluated eight inference algorithms and 18 references for cell-cell communication analysis. In our results, we outlined the framework of the analysis pipeline and interactive platform. Through manual annotation, we identified 36 distinct cell populations, highlighting the limitations of automated reference-based annotations in resolving cell states missing from underlying references - particularly injury-induced monocyte-derived macrophages and eosinophils. In our dataset, Log-normalization outperformed SCTransform in both preserving cell-type separation and batch mixing, while FastMNN provided the best balance between technical batch correction and the retention of cell subtype identities. Cell-cell communication analysis revealed that permutation-based algorithms used by CellChat and CellPhoneDB improved prediction confidence, while MouseConsensus and OmniPath resources increased interaction coverage. These findings demonstrate that preprocessing, integration, and cell annotation strategies should be selected in a context-dependent, data-driven manner rather than relying on default workflows. The annotated murine lung dataset expands the existing atlas, and the modular interactive platform offers a practical framework to improve accessibility, reproducibility, and analytical precision in pulmonary research.

Single-cell RNA sequencing (scRNA-seq) has emerged as an essential tool for characterizing cell type diversity and elucidating disease-associated mechanisms in acute and chronic pulmonary injury (Reyfman et al., 2019; Aran et al., 2019; Adams et al., 2020; Habermann et al., 2020; Wang et al., 2024). Over the past decade, rapid advancements in scRNA-seq technologies have been paralleled by the development of a wide array of computational methods for dimensionality reduction, clustering, cell-type annotation, integration, trajectory analysis, and inference of cellular communication (Dest et al., 2024; Arya et al., 2025; Molla Desta and Birhanu, 2025). Nonetheless, there is a paucity of comprehensive standardized pipelines that optimize the efficiency of these tools (Hou et al., 2020; Cheng et al., 2023).

A notable example pertains to data integration during the merging of multiple samples. Canonical Correlation Analysis (CCA), Reciprocal Principal Component Analysis (RPCA), Harmony, FastMNN (Fast Mutual Nearest Neighbors), and scVI (single-cell Variational Inference) are common methods designed to correct technical batch effects while preserving true biological variation (Butler et al., 2018; Hao et al., 2021; Korsunsky et al., 2019; Haghverdi et al., 2018; Lopez et al., 2018). However, the performance of these methods is not universally robust across all conditions (Cheng et al., 2023; Luecken et al., 2022; Tran et al., 2020). Use of different instruments, varying sequencing depths, and heterogeneity in sample dissociation and preparation across laboratories introduce complex technical noise (Luecken et al., 2022; Luecken and Theis, 2019; Tung et al., 2017). While dimension reduction techniques like Uniform Manifold Approximation and Projection (UMAP) are effective at visualizing large-scale batch effects (e.g., those arising from distinct sequencing depths), they remain qualitative diagnostics (Luecken et al., 2022; Luecken and Theis, 2019). Consequently, quantifying integration performance through objective metrics is essential to optimize data integration by minimizing technical variance without sacrificing the biological signal. With mild technical noise, determining whether integration is needed and which method to use is challenging because the analysis requires a delicate balance between alignment and preservation.

Cell-type annotation remains the most time-consuming step prior to downstream analyses, as current solutions heavily rely on reference datasets or established marker gene sets (Luecken and Theis, 2019; Abdelaal et al., 2019; Clarke et al., 2021). Although public data availability has improved significantly, high-quality reference datasets for disease-associated mouse lung tissue remain scarce. While the widely recognized Tabula Muris Senis (TMS) and LungMAP CellRefs (LungMAP) comprehensively cover the transcriptional landscape of the healthy lung across the lifespan, neither offers insights of lung injury and disease (Tabula Muris Consortium, 2020; Guo et al., 2023). This is significant because inflammatory disease is characterized by the influx of phenotypically active monocytes, which contribute directly to the inflammatory milieu and can replace the resident macrophage pool, alongside polymorphonuclear cells and T and B lymphocytes (Aran et al., 2019; Moos P. J. et al., 2024; Cembellin-Prieto et al., 2025). Prioritization of a specific cell type (e.g., alveolar or airway epithelial, fibroblast or myofibroblasts, interstitial or alveolar macrophages) influences the methodology of tissue dissociation (e.g. Dispase vs. Collagenase D, inclusion of DNase, use Fluorescence-Activated Cell Sorting) which are critical for cell viability and cell phenotype (Wang et al., 2024; Rodriguez et al., 2023; Rodríguez et al., 2025; Quintanal-Villalonga et al., 2022; Fujino et al., 2024). Thus, the choice of reference dataset and the underlying technical preparation may introduce significant bias, potentially compromising the accuracy of cell-type prediction and annotation across various methods.

Cell-cell communication (CCC) inferences provide a forecast of intercellular interactions in a technique that is inherently unable to capture post-transcriptional events. It involves leveraging computational models to integrate gene and gene-complex expression levels with curated ligand-receptor (LR) interaction resources. Heterogeneous curation of LR pairings across CCC resources (<15–20% overlap) and the paucity of species, tissue, or disease-specificity have been shown as constraints of its predictive power (Armingol et al., 2021; Jin et al., 2021; Efremova et al., 2020; Dimitrov et al., 2022). Evidence also highlights that the performance of these CCC inference tools is highly dependent on the computational method utilized, which produces conflicting “top-ranked” integrations when applied to the same dataset (Dimitrov et al., 2022). Thus, in the absence of systematic evaluation against biological ground truth or validated signaling networks, reliance on a single computational method or a resource may introduce significant bias, leading to the omission or misidentification of biological information.

Commercial platforms often provide user-friendly interfaces but frequently restrict the available methods for data normalization, integration, and CCC inference, or support direct comparisons between different algorithms (Best single; Wei et al., 2024). Among the six most popular scRNA-seq analysis platforms, only “Trailmaker” and “Partek Flow” allow selection of multiple integration algorithms, but neither platform quantifies the batch correction performance or the preservation of biological signal (Best single). Similarly, the incorporation of CCC inference tools by “Partek Flow” and “BBrowserX Pro” still lacks options for alternative resources or computational methods (CellPhoneDB; BBrowserX Pro). The code-based analysis pipelines offer greater flexibility but demand advanced programming and computational expertise, which is often a barrier for wet-bench researchers. Even for experienced programmers, manually repeating scripts during data exploration can introduce errors. Taken together, these limitations highlight the pressing need for an open-source platform that ensures a reproducible workflow, facilitates interactive data exploration, enables sharing with non-coding collaborators, and remains fully customizable for specific research requirements.

To address these critical gaps, our study introduces a high-resolution, annotated scRNA-seq dataset of a murine model of lung injury (Nureki et al., 2018). This resource extends beyond existing reference atlases by specifically incorporating multiple stages of inflammation and tissue remodeling. Utilizing this dataset as a representative case study, we demonstrate and explore key analytical steps in a data-driven scRNA-seq workflow, including cell-type annotation, normalization and transformation, data integration, and CCC inferences. In addition, providing the research community with a reusable and flexible platform, we developed an interactive R Shiny application. This platform enables researchers to explore complex lung fibrosis scRNA-seq data and perform validated analyses regardless of their programming and computational expertise. This comprehensive resource is strategically positioned to support future studies in lung biology, systems immunology, and computational method development.

Mouse model and sample collection

The samples utilized in these studies are obtained from surfactant protein (SP)-C mutant mice. Our group has previously developed and characterized the inflammatory and fibrogenic events occurring during the tamoxifen-dependent induction of the most common disease-linked mutation in the gene encoding for the SP-C (Nureki et al., 2018). Four experimental conditions (n = 2 for each) were utilized to validate our framework: air-exposed controls (Control); 2-week post-tamoxifen-induced SP-C mutant injury, consistent with the peak of inflammation and early lung remodeling (Inflam); 6-week post-tamoxifen-induced SP-C mutant injury, consistent with established fibrosis (Fib); 6-week post-tamoxifen induced fibrosis also challenged with the air pollutant ozone (Fib-O3). Following bronchoalveolar lavage collection, the perfused lung was collected and enzymatically processed (Collagenase D, 2 mg/mL), as previously described (Moos P. J. et al., 2024). Single cell suspensions were achieved through mechanical mincing of the tissue and passage through a 70 μm cell strainers. Single cell suspension and bronchoalveolar lavage were then combined and counted for the following sequencing steps.

Single cell RNA library preparation and sequencing

Samples were processed across five independent batches (Batch IDs: 24143, 20354, 21401, 22713, 23957) at the High-Throughput Genomics Shared Resource at the University of Utah (Huntsman Cancer Institute). Single-cell library preparation was performed using the 10x Genomics Chromium Single Cell 3′ GEM Library and Gel Bead Kit v3 to enable 3′ digital gene expression profiling. Suspensions containing approximately 800 to 16,000 cells per sample were loaded into a Chromium Chip B with partitioning oil, reverse transcription reagents, and gel beads containing 3.5 million unique 10x Barcodes. Gel Beads-in-emulsion (GEMs) were generated during a 9-min run on the 10x Genomics Chromium Controller, facilitating the co-encapsulation of single cells with barcoded beads. Within these droplets, poly-adenylated mRNA was reverse-transcribed, incorporating a unique 10x Cell Barcode and Unique Molecular Identifier (UMI) to each cDNA molecule for cell-specific transcript assignment. Following GEM breakage and library construction, library quality control was conducted, which included a Qubit dsDNA Assay, Agilent DNA ScreenTape Assay, and Kapa Library Quantification by qPCR to qualify and normalize the library pool. Final libraries were sequenced on an Illumina NovaSeq 6000 platform using a paired-end configuration. Sequencing was performed in blocks of at least 100 million read pairs to ensure sufficient depth for high-resolution transcriptomic analysis. Raw sequencing data (FASTQ files) were aligned to the GRCm38 (mm10) reference genome using the 10X Genomics Cell Ranger (v7.0) count pipeline.

Modular workflow and computational environment design

The interactive R Shiny platform was intentionally designed as a modular workflow to accommodate different computational environments while improving memory management and computational efficiency. Module I (individual sample quality control, filtering, and cell annotation) operates on individual samples and is designed to run on a personal computer. Module II (data integration and evaluation) and Module IV (cell-cell communication analysis) computing part contains the most computationally intensive analyses and is therefore recommended for Linux-based servers or high-performance computing (HPC) environments. To minimize computational requirements for downstream analyses, the outputs from Module II retain only the information required for subsequent analyses, enabling Module III (data analysis and visualization) and Module IV (cell-cell communication analysis) visualization part to run efficiently on personal computers. This modular design allows users to perform routine analyses locally while reserving server or HPC resources for computationally intensive integration tasks.

Data quality control (QC) and filtering

To ensure high data quality, cells were filtered based on the following criteria: mitochondrial gene expression <10%, ribosomal gene expression <50%, hemoglobin gene expression <1%, and unique feature counts >500 gene per cell. Potential doublets were identified and excluded using the DoubletFinder algorithm (McGinnis et al., 2019). To account for transcriptional heterogeneity driven by cell proliferation, were assigned cell cycle score using CellCycleScoring function in Seurat v5.0, utilizing the gene sets defined by (Hao et al., 2024; Tirosh et al., 2016).

Data transformation and scaling

To compare the impact of normalization on integration, two methods were utilized: Log-Normalization and SCTransform. Log-Normalization scales the raw counts to 10, 000 and natural-log transformed using “LogNormalize” option with NormalizeData function in Seurat (Butler et al., 2018). SCTransform applies a regularized negative binomial regression to stabilize variance and remove technical artifacts associated with sequencing depth (Chou et al., 2022; Hafemeister and Satija, 2019). For both methods, technical covariant (mitochondrial gene percentage and cell cycle score) were regressed to minimize their influence on following dimension reduction and clustering analysis.

Dimensionality reduction and clustering

Following normalization and regression of technical covariates, dimensionality reduction was performed using Principal Component Analysis (PCA) on the scaled expression matrix restricted to highly variable genes (HVGs). HVGs were identified within each dataset to enrich biologically informative features while minimizing stochastic noise from low-abundance transcripts. The top 30 principal components (PCs) were retained for downstream analyses based on variance structure and stability of clustering results across a range of dimensional cutoffs. PCA served to compress the high-dimensional transcriptomic space into a reduced representation that preserves dominant sources of biological variation while attenuating gene-level noise. This reduced space formed the basis for neighborhood graph construction and subsequent clustering. Cellular relationships were modeled using a Shared Nearest Neighbor (SNN) graph constructed in PCA space. By emphasizing shared local neighborhood structure, the SNN framework enhances robustness to technical noise and prevents artificial merging of transcriptionally distinct populations (Butler et al., 2018). Unsupervised clustering was performed in Seurat using the Leiden community detection algorithm applied to the SNN graph to identify groups of transcriptionally similar cells (Butler et al., 2018; Waltman and Van Eck, 2013). Clustering resolution was systematically tuned to balance lineage-level separation with detection of biologically meaningful subpopulations. For lineage-specific analyses, subsets were iteratively reprocessed (HVG identification, PCA, graph construction, clustering) to resolve finer cellular states. For visualization, Uniform Manifold Approximation and Projection (UMAP) was applied to the PCA space to generate two-dimensional embeddings. UMAP is a non-linear dimensionality reduction method that preserves local neighborhood relationships while maintaining the global structure of the data manifold, enabling intuitive visualization of cellular population structure in complex transcriptomic datasets (McInnes et al., 2020).

SingleR cell type prediction

We utilized SingleR to assign cell identities by calculating Spearman rank correlations between cells of our Sftpc I73T murine dataset and two established references: LungMap Mouse Lung CellRef (LungMap) and Tabula Muris Senis (TMS). LungMap is a comprehensive scRNA-seq reference build from 17 mouse lung samples that cover embryonic day (E) 16.5 to postnatal day (PND) 28 at eight developmental time points (E16 n = 2, E18 n = 2, PND1 n = 2, PND3 n = 2, PND7 n = 3, PND10 n = 2, PND14 n = 2, PND28 n = 2) (Guo et al., 2023). It encompassed two lineage levels (major cell type) and three refined cell type levels, up to 40 distinct mouse lung cell types (Guo et al., 2023). TMS is a large-resolution scRNA-seq atlas of aging in the mouse from 1 month to 30 months (1-month M n = 2, 3-month F n = 1, 3-month M n = 2, 18-month F n = 2, 18-month M n = 2, 21-month F n = 2, 30-month M n = 4) in 23 tissues, including the lung. It contains 30 different cell-type labels (Tabula Muris Consortium, 2020). In this study, we extracted droplet-based (10X Genomics) lung tissue data from the entire reference.

Manual cell type annotation

Annotation was refined through an iterative process within our R Shiny application. Initial major cell types were identified using canonical biomarkers. Subpopulations were further resolved by employing an interactive free-hand lasso tool (via Plotly) to select specific clusters for re-PCA, re-clustering, and re-UMAP until all subtypes were expertly identified and validated against differential expression results and biomarkers in the literature (Guo et al., 2023; Moos P. et al., 2024; Tabula Muris Consortium, 2018; Duong et al., 2022; Ghaedi et al., 2020; Zhang et al., 2022; Jirmo et al., 2023; Negretti et al., 2021).

Transcriptional similarity analysis

Transcriptional similarity among manually annotated cell populations was evaluated using aggregated expression profiles. Raw counts were summed across all cells within each cell population using Seurat’s AggregateExpression function and subsequently log-normalized. Pairwise similarity was calculated using Pearson correlation across genes. The resulting correlation matrix was visualized as a bi-clustered heatmap using Euclidean distance and complete-linkage hierarchical clustering. Aggregation across cells reduced the influence of stochastic cell-level dropout on the similarity estimates.

Benchmarking of integration methods

We evaluated five commonly used integration methods implemented in the Seurat v5 framework to correct for technical variation across five experimental batches and four disease conditions: Canonical Correlation Analysis (CCA), Reciprocal Principal Component Analysis (RPCA), Harmony, Fast Mutual Nearest Neighbors (FastMNN), and single-cell Variational Inference (scVI). These methods were selected as they represent widely applied integration strategies available within a unified analytical framework. CCA is an anchor-based method that aligns datasets by identifying shared correlation structures across features and projecting them into a common space (Butler et al., 2018).

RPCA extends this framework by identifying anchors through reciprocal projection between datasets, improving computational efficiency, and reducing the risk of overcorrection in large or heterogeneous datasets (Hao et al., 2021). Harmony performs batch correction through iterative adjustment of cell embeddings, enabling alignment across batches while preserving global and local data structure (Korsunsky et al., 2019). FastMNN corrects batch effects by identifying mutual nearest neighbors across datasets and aligning cells based on local similarity, which helps maintain separation of biologically distinct populations (Haghverdi et al., 2018). scVI is a deep generative modeling approach that learns a batch-corrected latent representation directly from raw count data (Lopez et al., 2018). By modeling gene expression using a probabilistic framework, it accounts for technical variability while capturing underlying biological structure (Lopez et al., 2018). The resulting latent space is used for downstream visualization and clustering.

Quantitative evaluation of data transformation and integration

The Sftpc I73T dataset was generated across five batches, with Inflam samples (n = 2) originating from a single batch and Control, Fib, and Fib_O3 samples (n = 2 each) distributed across separate batches (Supplementary Table S1). To systematically assess the impact of preprocessing and integration, we evaluated 10 workflows comprising two data transformation methods (Log-Normalization and SCTransform) and five integration methods (CCA, RPCA, Harmony, FastMNN, and scVI). Integration performance was assessed using both qualitative and quantitative criteria. UMAP visualization was used for initial inspection of batch mixing and cluster structure, while quantitative evaluation was performed using Average Silhouette Width (ASW), Local Inverse Simpson’s Index (LISI), and Graph LISI (GLISI). These metrics were selected to jointly evaluate two key objectives: (i) effective removal of technical variation (batch mixing) and (ii) preservation of biological identity (cell-type separation). Manual cell-type annotations and batch labels were used as ground truth for evaluation. LISI-based metrics quantified local neighborhood diversity, with batch LISI (iLISI) assessing the extent of mixing across batches and cell-type LISI (cLISI) evaluating the purity of cell-type neighborhoods (Luecken et al., 2022). Graph LISI extended this framework to the Shared Nearest Neighbor (SNN) graph, providing a topology-aware measure of integration quality that is less sensitive to distortions from low-dimensional embeddings (Luecken et al., 2022). ASW was used as a complementary metric to evaluate both cell-type separation and batch mixing (Luecken et al., 2022). It quantifies the relative distance of each cell to its assigned group compared to neighboring groups, providing a measure of cluster cohesion and separation. For cell-type labels, higher ASW values indicate well-separated and cohesive clusters, reflecting strong preservation of biological identity. In contrast, for batch labels, lower ASW values indicate effective mixing, as cells are not preferentially grouped by batch. Thus, while ASW captures the same two evaluation objectives as LISI, its interpretation depends on the label context and the directionality of the score. To enable comparison across metrics with different scales, we applied a median-based ranking approach, assigning relative scores (1–10) to each workflow for each metric. These rankings were then aggregated to identify the optimal combination of transformation and integration methods (Table 2).

Cell-cell communication (CCC) inference

To evaluate the coverage of intercellular communication inference, we employed the LIANA framework to conduct a systematic comparison across multiple computational methods and ligand-receptor (LR) resources. The methodological comparison included eight distinct algorithms: Connectome, logFC, NATMI, SingleCellSignalR (SCA), CellPhoneDB, CytoTalk, CellChat, and iTALK. Simultaneously, we assessed the content and overlap of 18 independent LR databases, including Cellinker, OmniPath, Consensus, MouseConsensus, CellTalkDB, LRdb, ConnectomeDB2020, Baccin2019, iTALK, Ramilowski2015, EMBRACE, CellChatDB, HPMR, ICELLNET, Guide2Pharma, and Kirouac2010. To compare inference methods independently of the underlying resources, each method was evaluated across all 18 LR resources, and duplicate interactions identified across resources were collapsed before summarization. Conversely, to compare LR resources independently of the inference methods, each resource was evaluated using all eight methods, and duplicate interactions identified across methods were collapsed. The analysis focused on quantifying the structural differences in LR pair composition, specifically evaluating the diversity of identified ligands, receptors, and their respective multimeric complexes across frameworks. Method comparisons were summarized by the total predicted interactions, whereas resource comparisons additionally evaluated the total and unique ligands, receptors, ligand subunits, and receptor subunits. The intersections of unique and shared interactions, ligands, receptors, ligand subunits, and receptor subunits were visualized using UpSet plots. To highlight the primary drivers of variability or consensus among the tools, visualization was restricted to the top 20 intersection groups.

Ortholog mapping and resource integration

For human-based resources, mouse genes were mapped to human orthologs. To leverage the comprehensive human-based ligand-receptor (LR) resources provided within the LIANA framework, we performed a systematic ortholog conversion of the murine gene expression data. We restricted our analysis to one-to-one orthologs between Mus musculus and Homo sapiens to ensure high confidence signaling predictions and avoid ambiguities associated with many-to-many gene relationships. Gene symbols were mapped using the internal mapping dictionaries provided by the OmniPath and LIANA R packages. During this conversion process, a total of 16,578 mouse genes were excluded due to a lack of a strict one-to-one human match. The resulting filtered expression matrix was then utilized as input for 17 human-centric resources. For the MouseConsensus resource, the analysis was performed using native murine gene symbols, identifying a core set of 96 CCC-related genes specifically involved in the mouse lung fibrosis context.

Software and environment

All analyses were conducted in R (v4.3.3) using RStudio (2025.5). The interactive platform was developed using Shiny. Key R packages included: Seurat (v5.3.1), SeuratObject (v5.3.0), SeuratWrapper (v0.3.0), SingleR (v2.10.0), DoubletFinder (v2.0.4), LIANA (v1.5.0), CellChat (v2.1.2), ggplot2 (v4.0.0), plotly (v4.10.4), glmGamPoi (1.18.0), DT (v0.35), tidyverse (v2.0.0).

Single cell RNA sequencing analysis pipeline development

We established a standardized scRNA-seq analysis pipeline consisting of five sequential stages: (1) Quality Control (QC) and Filtering, (2) Cell-Type Annotation, (3) Data Integration, (4) Integration Evaluation, and (5) Downstream Analysis and Visualization (Figure 1A).

To process individual samples, the initial step applies stringent filters to exclude low-quality cells based on the number of detected features (nFeature), the total number of transcripts (nCount), and the percentage of reads mapping to mitochondrial, ribosomal, and hemoglobin genes (Mitochondrial%, Ribosomal%, Hemoglobin%). Doublet identification and removal are performed using the DoubletFinder algorithm (McGinnis et al., 2019). Next, cell types are annotated through either an automated or manual process. SingleR, an automated annotation tool embedded, assigns identities by calculating the highest correlation between individual single-cell transcriptomes and reference datasets. The manual process requires a combination of data processing methods, including PCA, clustering, and UMAP; cells are labeled using known biomarkers (Aran et al., 2019).

Once annotated, processed data from individual samples are merged and batch-corrected through the Data Integration and Integration Evaluation stages of the pipeline. The final Downstream Analysis and Visualization stage of the pipeline is designed to address core biological questions through differential expression (DE) analysis, pathway enrichment analysis, CCC inferences, supported by visualization techniques for both gene expression levels and high-level analytical results.

Interactive analysis platform development

To standardize this pipeline and ensure it is accessible, reproducible, and time-efficient, we developed an interactive R Shiny application. The application is organized into four functional modules: (I) Individual sample QC, Filtration, and Cell Annotation, (II) Data Integration and Evaluation, (III) Data Analysis and Visualization, and (IV) Cell-Cell Communication (Figure 1B). The first module performs initial QC filtering and doublet removal, followed by data preprocessing, PCA, clustering, and UMAP generation. An interactive UMAP interface allows users to zoom, pan, overlay biomarker expression levels, or SingleR prediction results. It further enables manual selection of cell subsets for refined reanalysis and subtype annotation (e.g., clustering, UMAP generation, DE analysis).

Processed data from individual samples are merged in the second module, where the user can choose from multiple integration methods based on the specific requirements for batch correction. To guide the selection of the optimal method, the platform provides an interactive UMAP display, overlaying batch identity and integrated coordinates, alongside quantified evaluation scores for batch-mixing and cell-type separation. The third module facilitates deep data exploration and figure generation. It offers a range of analytical tools and customization options for producing high-resolution figures suitable for reports, grants, and publications. The final module is dedicated to CCC inferences. The App features both CellChat v2, the primary visualization and analysis tool, alongside the LIANA framework, which incorporates eight computational methods and 18 interaction resources (Dimitrov et al., 2022; Jin et al., 2025).

Sftpc I73T lung reference dataset construction and characterization

We characterized a comprehensive reference dataset using a previously established murine model of inflammatory fibrogenic exacerbations, driven by a clinically relevant mutation in the Surfactant Protein C gene (Sftpc I73T) (Nureki et al., 2018; Cheminant et al., 2025). The study investigated four distinct conditions: steady state (Control), inflammation and early remodeling occurring 2 weeks post injury (Inflam), established fibrosis at 6 weeks (Fib), and fibrosis with ozone exposure (Fib_O3) (Supplementary Table S1). Our analysis annotated 36 distinct cellular identities encompassing both unique cell types and distinct activation states of common populations (Figures 2A,B; Table 1; Supplementary Table S2). Cell identities were confirmed by evaluating top unique gene signatures for each population (Figure 2C). Cellular profiling revealed 22 annotated subtypes of immune cells alongside a sparse recovery of alveolar and ciliated epithelial cells, likely stemming from the mild digestive profile of Collagenase D (Table 1).

Mapping of Alveolar Macrophages (AM), Monocyte-Derived Macrophages (MoDM), and Interstitial Macrophages (IM) into a FastMNN-integrated UMAP space revealed evident population shifts across conditions (Figure 2B). As expected, MoDMs were largely absent in steady-state Control but exhibited a progressive increase through the Inflam and Fib stages, reaching peak abundance in the Fib_O3 group (Figure 2B). Moreover, the UMAP embedding highlighted the transient appearance of blood-derived eosinophils in the Inflam study group, underscoring their specific recruitment during the early remodeling phase (Moos P. J. et al., 2024).

Bi-clustered heatmap showed strong positive transcriptional similarity among myeloid (macrophages, monocytes, and dendritic cells) and lymphoid populations (B, T, natural killer cells, innate lymphoid cells) (Figure 2D). In contrast, plasma cells segregated independently of the broader B cell cluster. This separation reflected their specialized secretory transcriptome dedicated to antibody production (Luecken and Theis, 2019). Similarly, neutrophils and eosinophils formed a distinct cluster sequestered from other immune populations. This independent grouping is attributed to the unique granulocyte gene signature and low RNA content inherent to these cell types (Figure 2D) (Zilion et al., 2019). Among the non-immune compartments, Endothelial Progenitor Cells (EPCs) clustered more closely with stromal populations than with mature Endothelial compartments (Figure 2D). This likely reflects a shared mesenchymal-like transcriptional program characteristic of fibrosis (Tsukui et al., 2020). Finally, megakaryocytes and ciliated epithelial cells showed generally relatively low correlations with other populations, though ciliated cells maintained moderate correlation with alveolar epithelial (Alv_Epi) and club cells (Figure 2D).

Impact of references atlas composition on automated annotation performance

To assess the utility of existing community resources for automated annotation of injured lung dataset, we utilized SingleR to label our Sftpc I73T dataset using two established references: TMS and LungMAP. We then used our manual annotations as the ground truth to evaluate the predictive performance of these automated annotations. SingleR identified 29 cell types using TMS and 37 using LungMAP. The matched nomenclature across these annotations is summarized in Supplementary Table S3. While a common set of 17 major cell types was shared across the manual and automated annotations, only B cells, NK cells, and lymphatics endothelial cells demonstrated high predictive sensitivity and precision across both references. Injury-specific populations arising from blood populations (MoDM and eosinophils) were uniquely identified through our manual annotation. The majority of MoDMs were misclassified as either AMs or IMs by both references; monocyte-derived dendritic cells (MoDCs) were labeled as monocytes by TMS and as Type 2 Conventional Dendritic Cells (cDC2) by LungMap. This inaccurate annotation led to high sensitivity but low precision of AM predictions and low sensitivity and precision of IM prediction in both references. Both references also misidentified eosinophils for neutrophils, likely due to shared granulocytic signatures that are not sufficiently distinguished in steady state atlases (Figures 2A, 3A,B).

The TMS reference exhibited significant limitations across immune and endothelial lineages (Figure 3C; Supplementary Table S4). Among monocyte subtypes, TMS frequently misclassified the manually annotated classical monocytes (cMono) as “intermediate,” which led to low sensitivity of the cMono prediction, but non-classical monocytes (ncMono) were accurately predicted (Figures 2A, 3A,C). Additionally, TMS failed to distinguish between Type 1 Conventional Dendritic Cells (cDC1) and cDC2 subtypes and missed multiple immune populations, such as ILCs and Mast cells (Figures 2A, 3A,C). Predictive performance was also poor for Artery, Capillary General (Cap_G), and Alveolar Capillary (Cap_A) endothelial cells, which were mislabeled as “smooth muscle cells of the pulmonary artery” and “bronchial smooth muscle cells” (Figures 2A, 3A,C). Finally, TMS annotated Mesothelial and Endothelial Progenitor Cells (EPC) as fibroblasts due to the absence of these cell types in it, which reduced predictive precision of the fibroblast category (Figures 2A, 3A,C).

The LungMAP reference demonstrated stronger predictive power among epithelial, stromal, and endothelial compartments but exhibited deficient performance in predicting immune cell types (Figure 3D; Supplementary Table S5). Beyond its inability to resolve MoDMs and eosinophils, LungMap failed to identify MoDCs, ncMono, Plasmacytoid Dendritic Cells (pDC), plasma cells, Double-Negative T cells (DN_T), and NKT cells (Figure 3D). Lastly, LungMAP mislabeled Cap_G cells as Vein, Artery, and EPC, while two subtypes of Smooth Muscle Cells (SMC) and Pericytes were erroneously annotated as Fibroblasts (Figures 2A, 3B).

Benchmarking data preprocessing and integration workflows

To mitigate technical variation across the five collection batches, we conducted a systematic benchmarking of ten integration workflows pair with two data transformation and five integration methods (Supplementary Table S1). In UMAP embeddings, only the SCTransformed unintegrated data exhibited prominent batch-related color shifts (Figures 4A–F; Supplementary Figure S1). Our quantitative evaluation of batch-effect removal (batch-mixing) and preservation of biological identity (cell-type separation) were based on Average Silhouette Width (ASW), Local Inverse Simpson’s Index (LISI), and Graph Local Inverse Simpson’s Index (GLISI) quantitative matrices. Experimental batch information generated during library preparation and sequencing was used as the batch labels for batch-mixing evaluation, whereas manually annotated cell-type labels established prior to data transformation and integration were used to calculate the cell-type separation score as an independent measure of biological signal preservation. The median was applied to represent each matrix (Table 2). Under the unintegrated circumstances, Log-normalization provided a better batch mixing and cell-type separation rank scores compared to SCTransform across most of the quantitative matrices (Figures 4G,H; Table 2). When applying the same integration method across both transformation methods, log-normalization consistently yielded higher performance rank scores, suggesting it provided a more stable foundation for our dataset (Figures 4G,H; Table 2).

Among the methods tested, RPCA achieved the highest aggregated rank scores, offering the most balanced between batch-mixing and cell-type separation quantitatively (Figures 4G,H; Table 2). When analyzed individually, Harmony demonstrated the most aggressive batch correction, whereas FastMNN proved to be the best for preserving distinct cell-type separation (Figures 4G,H; Table 2). FastMNN was utilized to display our Sftpc I73T dataset, because it provided the clearest resolution of macrophage subtypes across different conditions, with minimal overlap in UMAP space (Figures 2A,B; Supplementary Figures S2, S3). Surprisingly, scVI yielded lower overall rank scores but preserved condition-specific macrophage subtypes in UMAP (Table 2; Supplementary Figure S4). It produced narrow distributions of Batch ASW (basw) and Cell ASW (casw) and significantly reduced the skewness of cell-type GLISI (gclici) compared to other methods (Figures 4G,H; Table 2). This suggests that scVI exerts strong global regularization on the latent space.

CCC inference methods and resources benchmarking via LIANA

Given that CCC inferences are influenced by LR database curations (resource) and computational models (method), we applied the LIANA framework to systematically evaluate their impacts on the Sftpc I73T dataset (Dimitrov et al., 2022). We evaluated five re-implemented methods (connectome, logfc, natmi, sca, and cellphonedb) and three external methods (CytoTalk, CellChat, and iTalk) as defined in the LIANA framework. The re-implemented methods yielded a consistent number of total interactions, and they were shared across methods (Figure 5A). Only cellphonedb accounted for multi-subunit complexes, while other re-implemented methods were post-aggregated (Dimitrov et al., 2022). Among the external methods, CytoTalk and CellChat incorporated multi-subunit complexes (Figure 5A). Notably, iTalk predicted an extremely large amount of total and unique interactions, likely due to its reliance on expression enrichment rather than stringent matching criteria (Figure 5A) (Dimitrov et al., 2022; Wang et al., 2019). CellChat and cellphonedb provide a higher confidence signaling landscape by performing permutation tests to mitigate false-positive predictions (Dimitrov et al., 2022).

Among the 18 accessible resources through LIANA, MouseConsensus is the only natively murine; the remaining 17 are human-derived. To facilitate a broader comparison of these resources regardless of reference method, we mapped the Sftpc I73T dataset to human orthologs and converted back after the prediction generated, retaining only genes with one-to-one matches. OmniPath yielded the highest total number of interactions and second highest count of unique interactions. It shared a large core of interactions with Celllinker, MouseConsensus, and Consensus (Figure 5B). Consensus featured the largest group of unique interactions, nearly double that of the next highest resources, indicating a high degree of exclusive content (Figure 5B). CellLinker involved the greatest diversity of unique ligands and receptors within the Sftpc I73T dataset (Supplementary Figure S5, S6). MouseConsensus, Consensus, CellChatDB, CellPhoneDB, ICELLNET, Baccin 2019, and Cellinker provide complex-containing interactions (Supplementary Figure S5B, S6B). MouseConsensus demonstrated enrichment in total and unique interactions, ligands, receptors, and multi-subunit complexes, highlighting the biological value of species-specific curation (Figure 5B; Supplementary Figure S5, S6). Furthermore, while the Baccin2019 resource was smaller in size, it exhibited a high percentage of unique interactions, ligands, receptors, and multi-subunit complexes compared to larger consensus resources (Figure 5B; Supplementary Figure S5, S6).

Standardized analytical pipelines for scRNA-seq data are essential for improving computational efficiency and ensuring reproducible discovery across datasets. However, existing frameworks often lack the flexibility to support data-driven method optimization tailored to specific biological contexts, such as interspecies and organ-level variability, or disease models. To address these challenges, our study was designed with two primary objectives: first, the generation of a high-resolution and immune-focused scRNA-seq atlas detailing the progression of murine pulmonary injury; and second, the development of an interactive R Shiny platform to democratize advanced computational approaches and facilitate systematic context-aware method evaluation.

Our analysis revealed that automated cell-type prediction using healthy murine lung references (TMS and LungMAP) was inadequate to identify multiple immune sub-populations and disease-related monocyte-derived macrophages and eosinophils. Our manually annotated Sftpc I73T dataset attempts to addresses this gap by gathering information from publish scRNA-seq datasets of murine injury and therefore standardize the identification of populations engaged in sterile injury and lung remodeling (Reyfman et al., 2019; Adams et al., 2020; Habermann et al., 2020; Joshi et al., 2020). By combining cells from the Bronchoalveolar Lavage (BAL) with collagenase D digested lung tissue, our protocol intentionally enriched the dataset for immune populations, while the inclusion of multiple injury states expanded the diversity of observed cellular subtypes. Utilizing our atlas as a reference for automated cell annotation tools enhances the prediction power for future immune-focused queries in respiratory research. Despite the depth of our immune profiling, the Sftpc I73T dataset has limitations regarding mesenchymal, epithelial, and endothelial cell diversity due to the specific tissue processing methodology employed. A recent series of studies by Rodríguez et al. utilized a similar Tamoxifen-induced Sftpc I73T model with a primary emphasis on mesenchymal and epithelial populations via Dispase digestion and CD45 + cell depletion, provides a complementary reference for cross-study integration (Rodriguez et al., 2023; Rodríguez et al., 2025). This dataset also contains the steady-state and 2-week post-tamoxifen timepoint, which can serve as biological anchor for batch integration.

To expand the lineage map of the injured lung, the strategic integration of external datasets is necessary. A similar conceptual framework has been published by Sikkema et al., which focused on the re-annotation of human samples from the Human Lung Cell Atlas (Sikkema et al., 2023). However, such an endeavor brings up challenges in optimizing dataset selection, mitigating batch effects, and maintaining computational efficiency as data volume increases. To maintain a compact and performant reference, we recommend downsampling the reference to include only samples specifically matched to the query dataset. Recent reports support the notion that restricting the reference to highly variable features (HVGs) effectively reduces the feature space without sacrificing resolution (Zhao et al., 2025). More broadly, the performance of reference-based cell annotation is fundamentally determined by the biological composition of the reference atlas, not solely by the prediction algorithm itself. Since no single reference represents all biological states, our pipeline combines automated multi-atlas annotation with manual refinement using established marker genes. This strategy leverages the strengths of diverse references while minimizing the limitations of any single atlas for heterogeneous or injury-associated cell populations.

Although SCTransform is often recommended for its ability to stabilize variance and account for varying sequencing depth, our results reveal a notable departure from its widespread adoption as the preferred preprocessing methodology (Hafemeister and Satija, 2019; Germain et al., 2020). In our experiments using the Sftpc I73T dataset, Log-normalization outperformed SCTransform in both batch mixing and cell-type separation. This discrepancy likely stems from their distinct mathematical foundations: whereas preserves the relative expression structure after Log-normalization, allowing downstream highly variable gene selection and PCA to capture biologically informative variation; SCTransform’s regularized Negative Binomial model scales the Pearson residuals to a unit variance of one across all features (Hafemeister and Satija, 2019; Townes et al., 2019). This global variance stabilization may increase the influence of sampling noise in lowly expressed genes (Luecken et al., 2022), thereby distorting the biological signal in PCA and UMAP space and making it harder for integration tools to distinguish batch effects from cell-type-specific biological signals (Ahlmann-Eltze and Huber, 2023). Furthermore, because high-performing methods like FastMNN were originally optimized log-normalized data, the choice of transformation fundamentally dictates the success of downstream integration (Haghverdi et al., 2018). These findings suggest that dataset-specific evaluation of normalization and transformation methods is necessary to preserve biological resolution in downstream analyses.

The manually curated cell-type labels were established on individually processed samples before dataset merging and integration workflows. Our comparison of five widely utilized integration methods revealed a clear trade-off between technical batch correction and biological preservation. The contracting performance of Harmony and FastMNN highlights the risks of excessive batch correction on resolution of subtle cell states and rare populations. We found that quantitative metrics such as LISI, Graph LISI, and ASW should be used in conjunction with visualization inspection tools (e.g. UMAP) for optimal integration method selection. Notably, although RPCA achieved the highest aggregate quantitative rank by balancing batch-mixing and lineage preservation, FastMNN ultimately provided the most biologically meaningful UMAP embedding. This discrepancy highlights a limitation: aggregate scoring systems often weigh each metric equally, which can lead to results biased by subtle numerical differences or specific experimental conditions. In our study, the “Inflam” study group was collected as a single batch. This created a unique challenge since condition-specific populations (e.g. eosinophils) inherently penalize “batch-mixing” scores, as there are no corresponding cells in other batches available for mixing. The scVI has been reported as one of the best integration methods for large-scale atlases and cross-species context (Luecken et al., 2022; Song et al., 2023). The poor performance in our study suggests strong global regularization of latent space, likely reflecting the influence of abundant cell populations, such as macrophages. Further, our relatively small dataset contained less complex batch effect compared to large-scale atlases, which may have led to overfitting of the scVI model. The result of our comparison suggests that manual, sample-specific annotation performed prior to data merging and integration, together with the use of cell-type labels, an independent reference, are essential for evaluating integration performance, ensuring that rare or condition-specific identities are not sacrificed for global batch-homogenization.

In this work, we also comprehensively investigated multiple CCC inference methods using the LIANA framework and found that CellChat and cellphonedb methods employ permutation testing to reduce false-positive interactions among eight tested methods. However, the LIANA framework does not automatically use the most recently updated methods and resources of CellPhoneDB and CellChat. The Python-based CellPhoneDB v5 maintains its high-precision subunit architecture and expands its coverage to ∼3,000 interactions, including non-protein molecules such as steroid hormones and metabolites (Troulé et al., 2025). Its new CellSign feature is particularly powerful for linking ligand-receptor pairs to downstream transcription factor (TF) activity (Troulé et al., 2025). The R-based CellChat v2 models the global flow of information using mass-action logic and cofactors, and most importantly its paired CellChatDB includes a dedicated murine-based resource (Jin et al., 2021; Jin et al., 2025). Owing to its native murine resources, extensive visualization capabilities, and integrated network and pathway analysis tools, CellChat v2 was embedded within our R Shiny platform. Despite these benefits, our results also show that CellChat v2 generated relatively low total interaction output among the methods tested because of the stricter statistical modeling and pathway level aggregation. Our group has previously published CellChat-based cell-cell interactions in the inflamed lung. While this work predicted strong complement signaling, our expression analysis also revealed no ligand-receptor pairings for C1qa, C1qb, C1qc, resulting in a distorted network prediction for interstitial macrophages (Moos P. J. et al., 2024). These findings provide practical guidance for the CCC component of our analysis pipeline. To reduce the risks of missing potentially significant interactions for murine lung studies, we recommend pairing CellChat and cellphonedb methods with expanded databases, such as OminPath, MouseConsensus, Consensus, Cellinker, and Baccion2019, via the LIANA framework.

In this study, we established a high-resolution single-cell atlas of the Sftpc I73T murine model, providing a critical reference for the complex immune landscape characteristic of acute and chronic pulmonary injury. Our systematic benchmarking demonstrates that data normalization, transformation, and integration should be selected in a data-driven, context-dependent manner to balance between batch effect removal and biological signal preservation. Furthermore, we showed that combining statistically rigorous tools (CellChat, cellphonedb) with the broader molecular coverage of LIANA provides a practical strategy for improving qualitative and quantitative CCC inference in murine datasets. To facilitate the broader application of these results, we developed an interactive R Shiny platform that democratizes access to our optimized analytical pipelines. Ultimately, this study provides both a biological roadmap and a data-driven analytical framework for respiratory single-cell research, enabling researchers to systematically evaluate and select appropriate computational methods for their own datasets while maintaining a reproducible workflow.

This R Shiny application provides an accessible implementation of the proposed analytical framework, enabling researchers with limited programming experience to systematically evaluate alternative analytical strategies. By enforcing a traceable workflow, this platform promotes reproducibility and advances data exploration. Our future work will focus on four key directions: 1) validation of the tool with other murine models of pulmonary disease; 2) incorporation of spatial transcriptomics datasets by implementing advanced tools like LIANA+ (Dimitrov et al., 2024); 3) refinement of the machine learning models for automated annotation and biomarker discovery; and 4) extrapolation and authentication of murine cellular annotation markers and signaling pathways in the human lung. In parallel with these efforts, expansion of this application to extra-pulmonary fibrosing systems (e.g. liver, kidney) may promote scientific collaboration and further accelerate discovery.