Table of Contents

Main Description File Descriptions Linked Files Installation and Instructions

1. Main Description

This is the Zenodo repository for the manuscript titled "A TCR β chain-directed antibody-fusion molecule that activates and expands subsets of T cells and promotes antitumor activity.". The code included in the file titled marengo_code_for_paper_jan_2023.R was used to generate the figures from the single-cell RNA sequencing data. The following libraries are required for script execution:

Seurat scReportoire ggplot2 stringr dplyr ggridges ggrepel ComplexHeatmap

File Descriptions

The code can be downloaded and opened in RStudios. The "marengo_code_for_paper_jan_2023.R" contains all the code needed to reproduce the figues in the paper The "Marengo_newID_March242023.rds" file is available at the following address: https://zenodo.org/badge/DOI/10.5281/zenodo.7566113.svg (Zenodo DOI: 10.5281/zenodo.7566113). The "all_res_deg_for_heat_updated_march2023.txt" file contains the unfiltered results from DGE anlaysis, also used to create the heatmap with DGE and volcano plots. The "genes_for_heatmap_fig5F.xlsx" contains the genes included in the heatmap in figure 5F.

Linked Files

This repository contains code for the analysis of single cell RNA-seq dataset. The dataset contains raw FASTQ files, as well as, the aligned files that were deposited in GEO. The "Rdata" or "Rds" file was deposited in Zenodo. Provided below are descriptions of the linked datasets:

Gene Expression Omnibus (GEO) ID: GSE223311(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE223311)

Title: Gene expression profile at single cell level of CD4+ and CD8+ tumor infiltrating lymphocytes (TIL) originating from the EMT6 tumor model from mSTAR1302 treatment. Description: This submission contains the "matrix.mtx", "barcodes.tsv", and "genes.tsv" files for each replicate and condition, corresponding to the aligned files for single cell sequencing data. Submission type: Private. In order to gain access to the repository, you must use a reviewer token (https://www.ncbi.nlm.nih.gov/geo/info/reviewer.html).

Sequence read archive (SRA) repository ID: SRX19088718 and SRX19088719

Title: Gene expression profile at single cell level of CD4+ and CD8+ tumor infiltrating lymphocytes (TIL) originating from the EMT6 tumor model from mSTAR1302 treatment. Description: This submission contains the raw sequencing or .fastq.gz files, which are tab delimited text files. Submission type: Private. In order to gain access to the repository, you must use a reviewer token (https://www.ncbi.nlm.nih.gov/geo/info/reviewer.html).

Zenodo DOI: 10.5281/zenodo.7566113(https://zenodo.org/record/7566113#.ZCcmvC2cbrJ)

Title: A TCR β chain-directed antibody-fusion molecule that activates and expands subsets of T cells and promotes antitumor activity. Description: This submission contains the "Rdata" or ".Rds" file, which is an R object file. This is a necessary file to use the code. Submission type: Restricted Acess. In order to gain access to the repository, you must contact the author.

Installation and Instructions

The code included in this submission requires several essential packages, as listed above. Please follow these instructions for installation:

Ensure you have R version 4.1.2 or higher for compatibility.

Although it is not essential, you can use R-Studios (Version 2022.12.0+353 (2022.12.0+353)) for accessing and executing the code.

1. Download the *"Rdata" or ".Rds" file from Zenodo (https://zenodo.org/record/7566113#.ZCcmvC2cbrJ) (Zenodo DOI: 10.5281/zenodo.7566113).
2. Open R-Studios (https://www.rstudio.com/tags/rstudio-ide/) or a similar integrated development environment (IDE) for R.
3. Set your working directory to where the following files are located:

marengo_code_for_paper_jan_2023.R Install_Packages.R Marengo_newID_March242023.rds genes_for_heatmap_fig5F.xlsx all_res_deg_for_heat_updated_march2023.txt

You can use the following code to set the working directory in R:

setwd(directory)

1. Open the file titled "Install_Packages.R" and execute it in R IDE. This script will attempt to install all the necessary pacakges, and its dependencies in order to set up an environment where the code in "marengo_code_for_paper_jan_2023.R" can be executed.
2. Once the "Install_Packages.R" script has been successfully executed, re-start R-Studios or your IDE of choice.
3. Open the file "marengo_code_for_paper_jan_2023.R" file in R-studios or your IDE of choice.
4. Execute commands in the file titled "marengo_code_for_paper_jan_2023.R" in R-Studios or your IDE of choice to generate the plots.

Additional file 1: Figure S1. Immune cell expression heterogeneity and cell subsets distribution across patients, related to Fig. 1. (A) UMAP projection of 2068 single T cells (left), 515 B cells (middle) and 126 macrophages (right) from 19 patients. Each dot corresponds to one single cell, colored according to cell cluster. (B) Heatmap of T cell clusters (left), B cell clusters (middle) and macrophage clusters (right) with unique signature genes. Top 20 specifically expressed genes are marked alongside, if available. (C-E) Bar plots showing the number (left panel) and fraction (right panel) of cells originating from the 19 patients for each subcluster of T cells (C), B cells (D) and macrophages (E). (F) The fractions of the 15 subclusters, NK cells, CAFs and endothelial cells in each patient. Figure S2. Cell subcluster characterization of functional status. (A) Top 100 ranked (based on fold change) differentially expressed genes indicative of the functional status in each T-cell cluster (top) and z-score normalized mean expression of known functional marker sets across single T cells (bottom). The numbers in parentheses correspond to the ranks and the key markers (Table S1) are highlighted by red color. (B) Heatmap showing the log2-transformed expression of selected T cell function-associated genes in single cells. (C) Violin plots showing the expression profile of selected genes involved in T-cell cytotoxicity (top) and exhaustion (bottom), stratified by T-cell clusters. (D) Top 100 ranked (based on fold change) differentially expressed genes indicative of the functional status in each cluster (C1, C2 and C3 for B cells; C0 and C1 for macrophages). The numbers in parentheses correspond to the ranks and the key markers (Table S1) are highlighted by red color. (E) Z-score normalized mean expression of known functional marker sets across single B cells (top) and the log2-transformed expression of selected B cell function-associated genes in single cells (bottom). (F-G) Heatmaps showing the z-score normalized mean expression of known functional marker sets across single macrophages and their log2-transformed expression in single cells. Blue boxes highlight the key markers and the numbers in brackets represent the total times appeared in literature. Figure S3. MM17 reference profile and performance assessment. (A) Heatmap of MM17 reference profile depicting z-score normalized expression of each gene across 17 tumor microenvironment (TME) cell subsets. (B-C) Correlation between predicted proportions and true proportions for each individual cell state (B) and for each individual patient (C). (D) Confusion matrix of all TME cell states. Figure S4. Functional associations of tumor microenvironment (TME) cell states. (A-D) Enriched GO biological processes of T_CD8_Cytotoxic (A), B_Non-regulatory (B), T_CD8_Mixed (C) and CAF (D) based on gene set enrichment analysis (GSEA). Figure S5. Associations between cell states and clinico-pathological variables. (A-C) Associations of molecular and clinical features with cell states. (A) Boxplots showing the cell fraction distribution of each cell state stratified by tumor type (left), gender (middle) and tumor status (right). (B) Boxplots showing the cell fraction distribution of each cell state stratified by integrative age (left), tumor stage (middle), and race (right). (C) The fraction distribution of cell states stratified by TCGA subtypes. Median value difference of cell fraction among subtypes was evaluated using Mood’s test. Wilcoxon rank sum tests were used to examine the significance of the differences between two groups. For tumor stage, patients with Stage 0, Stage I, IA, IB, Stage II, IIA, IIB and IIC are grouped as “LOW” (n=154), Stage III, IIIA, IIIB, IIIC and Stage IV are grouped as “HIGH” (n=162). * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. Figure S6. Associations between cell states and immune phenotypes, related to Fig. 4. (A) Scatterplots showing relationships between T_CD8_Cytotoxic and M_M2 (top), B_Regulatory and T_CD4_Exhausted (middle), CAF and T_CD8_Mixed (Cytotoxic and Exhausted) (bottom). Pearson correlations and p values are indicated. For significant correlations, linear models are shown as blue lines. (B) Contributions of the cell states to CA-1 (top) and CA-2 (bottom). (C) Scatter chart of the Pearson correlations of CA-1 and CA-2 with cell states. Different colors indicate whether or not significant associations between CA scores and cell states were observed (p < 0.05). (D) Boxplots showing the cell fraction distribution of each cell state stratified by the median values of CA-1 (top) and CA-2 (bottom), respectively. Wilcoxon rank sum tests were used to examine the significance of the differences between two groups. (E) The distribution of cell states across the three immunophenotype groups classified by median values of CA-1 and CA-2. Median value difference of cell fraction among groups was evaluated using Mood’s test. Then the statistical significance between any two groups was evaluated by Wilcoxon rank sum test and p values are shown at the top of each panel. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. Figure S7. Assessment on association between tumor microenvironment immune phenotypes (TMIP) and response to immune checkpoint blockade (ICB) in melanoma. (A) Box plots showing differences of CA-1 (upper panel) and CA-2 (middle panel) scores between responders and non-responders in patients under immunotherapy in TCGA data. Bar charts showing numbers of responders and non-responders with different TMIPs in those patients (lower panel). (B) Projection of each patient of Riaz et al. dataset onto the first and second component of the correspondence analysis. Left panel showed pre-treatment samples and right panel denoted on-treatment patients. Non-responders were colored blue, and responders were colored orange. Points denoted Ipi-naive patients, and triangles denoted Ipi-progressed patients. (C) Box plots showing differences of CA-2 scores between responders and non-responders in anti-PD1 pre-treatment patients (upper panel) and on-treatment patients (lower panel) who progressed after a first-line anti-CTLA4 treatment (Ipi-progressed) in Riaz et al. data. (D-E) Comparison of each cell state proportion between responders and non-responders in Ipi-progressed patients based on pre-treatment (D) and on-treatment (E) transcriptomic profiles. ns: not significant; *: p < 0.05. Table S1. Gene lists used for functional analyses. Table S3. Demographics and characteristics of patients with melanoma. Table S4. Uni- and multivariate analysis for progress-free survival (316 sample). Table S5. Uni- and multivariate analysis for overall survival (316 sample).

Supplemental Dataset S3 of: Barshad et al (2018). Genome Research, 28(7): 952-967.

Additional file 3 Supplementary Table S2. Table of the Spearman correlation coefficients (SCC) shown in Fig. 2d, h filled circles between the number of cells and the correlation scores. We used the Spearman correlation coefficients (SCC) of an imputation method (values in heatmap Fig. 2d) and the number of cells in the columns as inputs to calculate the Spearman correlation between the number of cells and the scores. (CSV 1 KB)

Additional file 1: Figure S1. HE staining and QC. (A) Samples from the AM_CTRL, AM_EM and AM_EC groups were HE stained, and the black arrow shows the gland invading the muscular layer. (B) X axis represents the number of UMI in each cell, and y-axis represents the number of genes in each cell. The distribution model is fitted according to the linear relationship. Yellow dots indicate cells that deviate from the threshold and will be removed in subsequent analysis. (C) The doublet score increases gradually from light color to dark blue (left figure), the red dots represent doubles predicted by Scrublet (right figure). The proportion of mitochondrial genes (D), the number of genes expressed (E), and the number of UMIs (F) in each cell before and after QC are shown in the violin plots. (G) The mean proportion of mitochondrial genes, mean number of genes expressed, mean number of UMIs in each cell and cell number of the three sample groups before and after QC are shown. HE, Hematoxylin-Eosin; QC, Quality Control. Figure S2. Cell type identification and heatmap of gene expression in clusters. (A) Seventeen clusters were displayed in the AM_CTRL, EAM_EM and AM_EC groups. (B) Heatmap showing the expression levels of specific markers in each cluster. (C) The cell number and percentage corresponding to each cell type were counted. Complementary representative markers of different cell types (D) and corresponding violin plots (E) are shown. Figure S3 Colocalization of epithelial cell and endothelial cell markers in cluster 1. (A, B) Complementary epithelial cell markers (CDH1 and KRT7), endothelial cell markers (VWF and CDH5) and colocalization of the two cell type markers are displayed in the t-SNE map. (C) Confirmatory colocalization of EPCAM and PECAM1 was conducted in additional AM_CTRL (V) (n=3), AM_EM (V) (n=3), AM_EC (V) (n=3) samples, EPCAM (red), PECAM1 (green) and nuclei (blue) were stained. The white arrows show the colocalized cells containing EPCAM and PECAM1, scale bars = 20 μm, *P < 0.05. GO (D) and KEGG (E) analyses of upregulated genes in cluster 1 compared with all other clusters (cluster 2 to cluster 17). Representative GO (F) and KEGG (G) terms of upregulated genes in cluster 1 compared with the epithelial cell population (cluster 7 and cluster 17) are displayed. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Figure S4 Gene expression changes and functional analyses of the AM_CTRL, AM_EM and AM_EC groups. Representative GO (A) and KEGG (C) terms of upregulated DEGs in the AM_EC group compared with the AM_EM group are shown. Representative GO (B) and KEGG (D) terms of upregulated genes in the AM_EM group compared with those in the AM_CTRL group are shown. (E) Venn diagram showing the common change genes between AM_EM versus AM_CTRL and AM_EC versus AM_EM. (F) Enriched GO terms of common changed genes in (E) are shown. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Figure S5. Identification of subclusters in epithelial cell populations. (A) The t-SNE map shows 7 subclusters in the epithelial cell group. Histograms show the proportion of subcluster 1 in AM_EM and AM_EC (B), and the proportion of subcluster 4 is shown in (C). (D) The heatmap displays the expression levels of specific markers in each subcluster. Representative GO (E) and KEGG (F) terms of the top 500 genes expressed in subcluster 1 are indicated. The corresponding GO (G) and KEGG (H) analyses of subcluster 4 are also presented. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Figure S6. Change in marker distribution on the pseudotime trajectory and verification of VM formation. (A) Distribution of groups on the pseudotime trajectory. (B) Cell type, including epithelial cell, cluster 1 and endothelial cell, distributed on the pseudotime trajectory. (C, D) Epithelial cell markers (EPCAM, CDH1 and KRT7) and endothelial cell markers (PECAM1, VWF and CDH5) distributed on the pseudotime trajectory. (E) Confirmatory VM channel formations were conducted in additional AM_CTRL (V) (n=3), AM_EM (V) (n=3), AM_EC (V) (n=3) samples, VM channels were positive for PAS staining but negative for CD34 (black arrows); red arrows represent channels positive for CD34; scale bars = 50 μm. (F) Histogram indicates the number of VM channels in different groups; ns represents not significant, *P < 0.05. VM, Vasculogenic Mimicry. PAS, Periodic Acid-Schiff.

Supplementary Material 1. Supplementary Fig. 1. Quality control and baseline data of each enrolled sample. (A). Principal component analysis before and after processing in Harmony. (B). t-SNE projections among different groups. (C). Quality control of the scRNA-seq data. (D). t-SNE projections from all enrolled samples. t-SNE in the control group (left), STS Pre group (middle), and STS Post group (right). (E). Boxplot comparing the proportion of plasmacytoid dendritic cells(pDCs) across the groups. The STS Pre vs. CT and STS Post vs. CT sample comparisons show exact P values determined by the Wilcoxon rank-sum test. Pre- vs. post-STS scores were calculated via the paired two-sample Wilcoxon signed-sum test. (F). The baseline information of patients and healthy controls enrolled in the scRNA-seq cohort. Supplementary Fig. 2. Focused analysis of T cells and pDCs. (A). UMAP embedding of T lymphocytes from all profiled samples in different groups. UMAP in the control group (left), STS Pre group (middle), and STS Post group (right). (B). Boxplot comparing the proportions of CRIP + CD4 + T cells, NK T cells, and TRGC2 + CD8 + T cells across the groups. The exact P values determined by the Wilcoxon rank-sum test are shown for the STS Pre vs. CT and STS Post vs. CT comparisons. Differences between STS Pre and STS Post were evaluated by the paired two-sample Wilcoxon signed-sum test. (C). Enriched pathways from Gene Ontology Biological Process Enrichment Analysis for TRGC + CD8 + T cells. (D). Enriched pathways identified by Gene Ontology Biological Process enrichment analysis in NEAT + T cells. (E). Heatmap representing the enrichment of MSigDB Hallmark gene sets for each T lymphocyte subtype across groups. (F). Pseudotime trajectory analysis of CD4 + T lymphocyte subtypes. (G). Heatmap represents DEGs within pDCs across groups. (H). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the MSigDB of each group within pDCs. Supplementary Fig. 3. Focused analysis of B cells and myeloid cells. (A). Heatmap representing the enrichment of Hallmark gene sets in the MSigDB for each cell type within B lymphocytes across groups. (B). Boxplot comparing the proportions of myeloid cells across the groups. The exact P values determined by the Wilcoxon rank-sum test are shown for the STS Pre vs. CT and STS Post vs. CT comparisons. Differences between STS Pre and STS Post were evaluated by the paired two-sample Wilcoxon signed-sum test. (C). Enriched pathways from Gene Ontology Biological Process Enrichment Analysis for CD16 + monocytes. (D). Heatmap representing the enrichment of MSigDB Hallmark gene sets in each monocyte cell type across groups. Supplementary Fig. 4. The characteristics of IFN-related genes involved in pathogenesis. (A). Heatmap showing the differentially expressed genes (DEGs) in classical dendritic cells(cDCs) across groups. (B). Heatmap showing the genes differentially expressed in mast cells across groups. (C). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the mast cells across groups. (D). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the cDC across groups. (E). Heatmap representing the enrichment of Hallmark gene sets in the MSigDB for each cell type within Natural killer (NK) cells across groups. (F). Venn plot of the overlapping genes downregulated in the STS Pre group among B lymphocytes, T lymphocytes, monocytes, NK cells, cDCs and pDCs. (G). Receiver operating characteristic (ROC) curve and area under the curve (AUC) of overlapping genes expressed at lower levels before treatment in all cell types. (H). T-SNE analysis of CXCR4 expression in the three groups. (I). The relative expression of CXCR4 across groups was determined through qPCR. Statistical significance is denoted as P < 0.05 (), P < 0.01 (), P < 0.001 (), or P < 0.0001 (). Supplementary Fig. 5. Analysis of transcription factors (TFs) in INS. (A). Heatmap showing TFs activation across the subcluster of T lymphocytes. (B). Heatmap showing TFs activation across the subclusters of monocytes, cDCs, neutrophils and mast cells. (C). Heatmap showing TFs activation across the subclusters of B lymphocytes. (D). Heatmap showing TFs activation across the subcluster among ΝΚ cells. (E). Heatmap showing the expression level of IFNs (IFNA1, IFNA2, IFNB1, IFNG, IFNL1, IFNL2, and IFNL3). in each cell type across groups. Supplementary Fig. 6. (A). Surface expression data of TACI and BCMA on naïve B cells, unswitched memory B cells, switched memory B cells and plasma cells determined via flow cytometry. The samples were obtained from INS patients. (B). Surface expression data of TACI and BCMA on naïve B cells, unswitched memory B cells, switched memory B cells and plasma cells determined via flow cytometry. The samples were obtained from healthy donors. (C). The proportion of naïve B cells in the serum of NS patients compared to that of healthy individuals. ns, p > 0.05. (D). The proportion of unswitched memory B cells in the serum of NS patients compared to that of healthy individuals. ns, p > 0.05. (E). The proportion of switched memory B cells in the serum of NS patients compared to that in the serum of healthy individuals. ns, p > 0.05. (F). The proportion of plasma cells in the serum of NS patients compared to that in the serum of healthy individuals. ns, p > 0.05. Supplementary Fig. 7. Supplementary analysis from an extra INS cohort (GEO233277) also validates the activation of IFN. (A). UMAP dimensionality reduction embedding from GEO datasets. (B). Heatmap showing the expression levels of the markers across each cell type using scRNAseq from the GEO datasets. The color intensity indicates the marker of interest. (C). Violin plot of the ISGs across groups using scRNAseq from the GEO datasets. Significance was evaluated with the Wilcoxon rank-sum test. (D). ISG scores among cell subtypes across groups using scRNAseq from the GEO datasets. Significance was evaluated with the Wilcoxon rank-sum test. Statistical significance is denoted as P < 0.05 (*), P < 0.01 (), P < 0.001 (), or P < 0.0001 (*). (E). Incoming signaling patterns of APRIL and BAFF across groups using scRNAseq from the GEO datasets. (F). Outgoing signaling patterns of APRIL and BAFF across groups using scRNAseq from the GEO datasets. (G). Heatmap showing the expression level of all kinds of IFNs in each cell type across groups using scRNAseq from the GEO datasets. Supplementary Fig. 8. The potential mechanism underlying the pathogenesis of pediatric idiopathic nephrotic syndrome. (A). IFN-γ activation, mainly generated by T cells, stimulates the upregulation of BAFF expression among monocytes, dendritic cells, and neutrophils. Subsequently, the activation of BAFF interacts with its receptors on B cells, especially BCMA and TACI, facilitating B cell maturation and leading to autoantibody release. BAFF: B-cell activating factor; BCMA: B-cell maturation antigen; TACI: transmembrane activator and cyclophilin ligand interactor.

Supplementary Material 2: Additional file 2: Figure S2. Expression of classic markers in each T-cell subcluster and inferred TF activities of T-cell subclusters. (A) Bubble plot depicting the average expression levels and cellular fractions of selected stress-related heat shock marker genes in 3 defined T-cell subclusters. (B) Bubble plot depicting the average expression levels and percentages of exhaustion-related marker genes in 3 defined T-cell subclusters. (C) Bubble plot depicting the average expression levels and cellular fractions of effector function-related marker genes in 3 defined T-cell clusters. (D) Heatmap showing TF activity among different T-cell subclusters. (E) Heatmap showing TF activity among different T-cell subclusters in the CAG and control groups. CD4+ Tex: exhausted CD4+ T-cell subtype, CD8+ Tef: CD8+ T effector cell subtype, CD8+ Tsr: CD8+ T stress response.

This repository provides the processed data necessay to produce the figures from the Weiner, et al., Nature Communications paper entitled "Inferring replication timing and proliferation dynamics from single-cell DNA sequencing data".

The data is organized into directories based on figure numbers of the paper. For instance, the directory source_data/fig2_S3/ contains all data pertaining to main figure 2 and supplementary figure 3. Within each directory, subdirectories are organized by figure panel, meaning that the file source_data/fig2_S3/2fg_S3g/cell_metrics.tsv only pertains to data which appears in Fig 2f, Fig 2g, and Fig S3g.

Certain figure panels incorporate data from multiple samples. Sample ID subdirectories are used for these panels. For example, Fig S8a-c use data at the following paths source_data/fig5_S8/S8abc/{sample_id}/s_phase_bafs.csv.gz

All source code for upstream data preprocessing and downstream plotting can be found at the corresponding github repository for this manuscript: https://github.com/shahcompbio/scdna_replication_paper.

In addition to source data, we are also including two additional files which contain series of sample-specific heatmaps. Below are the descriptions of each file.

Additional File 1: HMMcopy and SIGNALS heatmaps of high-quality G1/2-phase cells across all samples. Matrices of somatic copy number state called by HMMcopy (Ha, et al, 2012) (left) and allelic imbalance state called by SIGNALS (Funnell, et al, 2022) (right) for all cells in a sample. Only high-quality G1/2-phase cells are included in this analysis. The rows are sorted the same in both heatmaps to preserve mapping of cell IDs. Clone IDs for all cells are shown using the colorbar to the left of both heatmaps. Each page represents a unique sample in the metacohort of breast and ovarian cell lines and PDXs (Fig 4a). The sample ID and number of high quality G1/2-phase cells (i.e. SIGNALS cells) are shown at the top of each page.

Additional File 2: PERT input and output matrices for gastric cancer cell lines at 500kb and 20kb resolution. Matrices of PERT input (left: reads per million and HMMcopy states) and PERT output (right: PERT somatic copy number and replication states) for the three gastric cancer cell lines sequenced with 10X Chromium single-cell DNA (Andor, et al, 2020). The top heatmaps on each page contains the cells predicted to be in S-phase by PERT and the bottom heatmaps contains the cells predicted to be in G1/2-phase by PERT. The rows are sorted the same in all four columns going from right to left to preserve mapping of cell IDs. Clone IDs for all cells are shown using the colorbar to the left of all four heatmaps. The first three pages represent PERT runs on each cell line at 500kb resolution. The final two pages represent PERT runs for two of the three cell lines at 20kb resolution. We did not run PERT at 20kb resolution for the SNU-668 cell line as there were too few S-phase cells to infer the RT profile at 20kb higher resolution.

For further information please reach out to Adam Weiner (weinera2@mskcc.org)

View the temporal single-cell transcriptomics data (UMAP, PCA, Heatmaps and Violin plots) using the Shiny App interface of iSEE (doi:10.12688/f1000research.14966.1) for easy visualization of the single-cell data described in "Single-cell roadmap of cardiac differentiation identifies roles for ZNF711 and retinoic acid in balanced epicardial and cardiomyocyte lineage commitment" (Snabel et al., bioRXiv).

For instructions on how to use this data, please visit https://github.com/Rebecza/scRoadmap_CardiacDiffs/.

(D) Heatmap of single-cell expression measured by qRT-PCR of major ESC and EpiSC markers in un-induced EpiSCs (black), established ESCs (red), Day 2 High/Low (dark and light blue), and Day 4 High cells (green). List of tagged entities: Multiple componant, , gene expression assay (bao:BAO_0002785),quantitative reverse transcription PCR (bao:BAO_0002090)

Additional file1: Figure S1. The information on MSC atlas taxonomy. (A) UMAP of all MSCs with cluster annotations, (B) UMAP of MSCs color-labelled by tissue, (C) Cell counts of MSCs from different tissues in each cluster, and (D) Cell counts of MSCs from different samples in each cluster. Figure S2. Differentiation scoring of MSCs on five differentiation directions. (A) Scoring of osteogenesis, chondrogenesis, adipogenesis, myogenesis and neurogenesis. (B) Scoring of representative gene expression for MSCs differentiation. Figure S3. Home page of MSCsDB. which includes website introduction, functionality overview, gene cloud, and website update news. Figure S4. Module of Dataset and link to the module of Explore. Users can view the metadata of each sample dataset, such as the original article, data repository and sequencing technology. Users can also click on the “Explore” button to view the sample’s clustering annotation, gene expression level analysis, pathway enrichment analysis, copy number variation analysis, and pseudotime analysis results. Figure S5. Functionality in the module of Atlas. (A) UMAP of MSCs with cluster annotations. Users can select specific clusters to view their distribution. The MSC atlas can also be classified by tissue or batch and shown separately. (B) Gene signature of MSCs. Users can analyze the cell percentage of all genes and click on the “View” button to view the gene expression levels in cells and clusters. The Gene Card database is also linked for users to view gene information. Users can also enter a specific gene in the search box to retrieve relevant information. Figure S6. An example of functionality in the module of Atlas. (A) Pathway enrichment analysis of MSCs from different databases. Users can switch between different databases. Users can also select specific clusters and pathways to view their enrichment status. (B) Copy number variation analysis of MSCs using copyKat and InferCNVpy packages. The copyKat software can predict whether the cells are normal cells (diploid) or tumor cells (aneuploid). The InferCNVpy package gives prediction values, so we provide chromosome heatmaps based on CNV clustering for users to distinguish between normal cells and tumor cells. (C) Pseudotime analysis of MSCs using PAGA method. We show the cell trajectory inference plot and cluster UMAP plot for a single sample. (D) Transcription factor network analysis of MSCs using pyscenic package. We provide the transcription factor network analysis result table and heatmap for a single sample’s cluster. Users can click on the “View” button in the table to view the target genes regulated by that transcription factor. Figure S7. De novo analysis for clustering, pathway enrichment, and quality evaluation. (A) UMAP plot of MSC clustering and annotation using Scanpy package for a sample dataset. (B) Pathway enrichment analysis using Clusterprofiler package for a sample dataset. (C) Copy number variation analysis using CopyKat and InferCNVpy packages for a sample dataset. Figure S8. De novo analysis for pseudotime and gene regulatory network analysis. (A) Pseudotime analysis using PAGA method for a sample dataset. (B) Gene regulatory network analysis using pyscenic package for a sample dataset. Table S1. Marker genes used for potency score analysis. Table S2. Scoring for each cluster using geneset.

Supplementary Material 4: Additional file 4: Figure S4. Gene expression of different fibroblast subpopulations and changes in the proportions of different fibroblast subpopulations between CAG and control tissues. (A) Heatmap of the top 10 genes of the cluster-defining DEGs among all fibroblasts from the control and chronic astrophic gastritis cell populations. Red on the heatmap shows the most highly upregulated genes, and blue shows the most downregulated genes. (B) Proportions of fibroblast subclusters in CAG and control samples. (C) Chemokine and cytokine molecule expression levels in different fibroblast subclusters of CAG. (D) Feature plot showing CCR5 expression in CD8+ T cells.

Table S4. List of genes used in Fig. 3a and c for heatmap. (XLSX 99 kb)

Additional file 1: Supplementary Figure 1: DDR subtyping analysis overview. A. DDR network overview. HR: Homologous recombination. MMEJ: microhomology-mediated end-joining. NHEJ: Non-homologous end-joining. BER: Base excision repair. NER: Nucleotide excision repair. MMR: Mismatch repair. DR: Direct reversal repair. TLS: Translesion synthesis. Checkpoint: Damage sensing and signaling. FA: Fanconi Anemia. Numbers in parentheses indicate number of pathway genes analyzed by our method. B. WE score formula and Essentiality Scaling Factor (ESF) criteria. Supplementary Figure 2: DDR cluster prevalence and DDR pathway single gene heatmaps. A. DDR cluster prevalence in GEMINI cohort. B. DDR cluster prevalence in IMPOWER133 cohort. C. GEMINI DDR pathway single gene expression heatmap. D. IMPOWER133 DDR pathway single gene expression heatmap. Supplementary Figure 3: IMpower133 DDR cluster differential gene expression and quantitative set analysis (QuSAGE) BG signature results. Supplementary Figure 4: SCLC CellMiner cell line DDR optimal number of k clusters elbow plot. Supplementary Figure 5: SCLC CellMiner cell line DDR cluster prevalence and marker expression. A. DDR cluster prevalence in SCLC cell line models. B. DNA damage responsive transcription factors, intra-S cell cycle checkpoint, and G2/M cell cycle checkpoint machinery expression across SCLC cell line DDR clusters. C. SLFN11 protein expression in SCLC cell line DDR clusters. D. Total RB1 and RB1-S807.S811 phospho protein expression across SCLC cell line DDR clusters. E. Cell line DDR cluster MYC family gene expression. Supplementary Figure 6: SCLC cell line DDR cluster mutation profiling. Oncoprints for TP53, RB1, and DDR gene mutations in DDR Low, Intermediate, and High clusters. Mutation data is from CCLE whole exome sequencing [17]. Supplementary Figure 7: SCLC cell line DDR cluster nonsynonymous tumor mutational burden (TMB). TMB data is from CCLE whole exome sequencing [17]. Supplementary Figure 8: SCLC/hgNEC PDX/CDX DDR gene expression and cell cycle state distribution. A. Expression of DNA damage responsive transcription factors, intra-S, and G2/M cell cycle checkpoint effectors in PDX/CDX DDR clusters. B. PDX/CDX DDR cluster scRNAseq cell cycle state distributions. Supplementary Figure 9: IMpower133 subtyping method comparison. A. Three way alluvial plot demonstrating assignment overlap between MDACC SCLC subtypes [5], Genentech (GNE) subtypes [7], and DDR status. B. IMpower133 DDR cluster and GNE subtype assignment table. Values listed represent the number of samples in each assignment category. The corresponding X2 Pearson residual dot plot comparing DDR cluster and GNE subtype assignments is shown on the right. Supplementary Figure 10: GEMINI and IMpower133 DDR cluster neuroendocrine score single gene heatmaps. A. GEMINI DDR cluster neuroendocrine score single gene expression heatmap. B. IMpower133 DDR cluster neuroendocrine score single gene expression heatmap. Supplementary Figure 11: SCLC CellMiner DDR cluster neuroendocrine features. A. SCLC cell line DDR cluster NE scores. B. SCLC cell line DDR cluster neuroendocrine status as reported by Tlemsani et al. C. SCLC PDX/CDX DDR cluster NE scores. Supplementary Figure 12: IMpower133 SCLC-A only immune checkpoint marker expression. Boxplots depict mRNA expression of CD274/PD1, CTLA4, and HAVCR2/TIM3 in SCLC-A tumors, split by their DDR status. Supplementary Figure 13: SCLC/hgNEC PDX/CDX DDR cluster MHC Class I scRNAseq expression. Supplementary Figure 14: IMpower133 DDR cluster chemotherapy response analysis. A. RECIST Best Overall Response (BOR) for IMpower133 DDR clusters following frontline EP chemotherapy. PD: Progressive disease. SD: Stable disease. PR: Partial response. CR: Complete response. B. Progression free survival Kaplan Meier plot for SCLC-A DDR clusters following frontline EP chemotherapy. C. Forest plot for progression free survival for SCLC-A DDR clusters following frontline EP chemotherapy. Supplementary Figure 15: IMpower133 DDR cluster all subtypes chemotherapy response analysis. A. Overall survival (OS) Kaplan Meier plot for all subtypes DDR clusters following frontline EP chemotherapy. B. Progression free survival (PFS) Kaplan Meier plot for all subtypes DDR clusters following frontline EP chemotherapy. C. Forrest plots for all subtype DDR cluster OS and PFS outcomes following frontline EP chemotherapy. Supplementary Figure 16: MDACC GEMINI subtype switching plots. Subtype space and SCLC-DMC subtype calls for MDACC GEMINI subtype switching patients. Supplementary Figure 17: MDACC GEMINI subtype switching patient outcomes following frontline chemoimmunotherapy. A. Overall survival outcomes following frontline chemoimmunotherapy. B. Progression free survival following frontline chemoimmunotherapy. Supplementary Figure 18: Promoter methylation changes for MDACC GEMINI subtype switching patients from baseline to progression following frontline chemoimmunotherapy. For all panels, data presented is RRBS methylation data previously published by Heeke et al [14]. Supplementary Figure 19: Circulating tumor DNA fraction for subtype switching tumors at baseline and progression following frontline chemoimmunotherapy. Supplementary Figure 20: SC53 Leiden cluster specific top drivers of plasticity. scRNAseq data for SC53 is from Gay et al. [5]. Supplementary Figure 21: WebGestalt over-representation analysis of SC53 Leiden Cluster 8 marker genes. Supplementary Figure 22: SC53 Leiden cluster specific random walk simulations. scRNAseq data for SC53 is from Gay et al. [5]. Supplementary Table 1: WE score individual DDR gene Essentiality Scaling Factor (ESF) assignments and example raw WE score generation. Supplementary Table 2: SCLC cell line whole exome sequencing mutation enrichment statistics by DDR cluster. Genes not listed in Supplementary Table 2 were not mutated in any DDR Low, DDR Intermediate, or DDR High cell line models. Supplementary Table 3: SC53 dynamic genes by Leiden cluster.

Additional file 1: Figure S1. Most insulin+ cells are ablated from 15dpDT and on. (A) Mean blood glucose levels for the mice used for the 5, 15 and 30 dpDT RNA-seq experiments. (B) Representative immunofluorescence staining of pancreas sections obtained from Control, 5dpDT and 15dpDT treated mice. (C) Percentage of cells expressing Insulin (Ins+) and YFP, expressed relative to the total number of cells in the islets. Total number of islet cells (expressed as mean ± SEM) counted per replicate: Control (784 ± 71, 3 biological replicates), 5dpDT (442 ± 111, 3 biological replicates), 15dpDT (301 ± 168, 2 biological replicates). Figure S2. Signaling pathways modulated in α-cells after acute β-cell loss. (A) Heatmap showing the Normalized Enrichment Score (NES) for selected Gene Set Enrichment Analysis (GSEA) results. T1: Ctrl vs 5dpDT α-cells, T2: 5dpDT vs 15dpDT α-cells, T3: 15dpDT vs 30dpDT α-cells. Genes with higher expression in Ctrl, relative to 5dpDT α-cells (i.e. enriched in Ctrl, with a positive NES), are associated with Smoothened signaling. Genes downregulated in 15dpDT, relative to 5dpDT α-cells (i.e. enriched in 5dpDT, with a positive NES), are associated with MET, mTOR and Rac1 signaling (among other pathways shown). Several of the pathways downregulated in T2 are upregulated in 30dpDT, relative to 15dpDT α-cells (see main text discussion for further details). All results presented are significant considering a P-value < 0.05 and FDR < 0.25. (B) Schematic summarizing the stepwise up- and downregulation of signaling pathways in α-cells as they transition from homeostasis (Ctrl α-cells) to 30dpDT. Figure S3. Epigenomic analysis reveals candidate genes and regulatory elements driving the α-cell response to β-cell loss. (A) Integrative epigenomic maps of the loci of selected cluster 1–5 marker genes showing human islet ChIP-seq signal for 5 key pancreatic islet transcription factors, histone modification enrichment profiles associated with active (H3K27ac) promoter (H3K4me3) and enhancer (H3K4me1) regions, islet enhancer hubs, islet regulome regions and arcs representing high-confidence pcHi-C interactions in human islets (data from (Pasquali et al. 2014 [18]) and (Miguel-Escalada et al. 2019 [26])). Note: OAS1 is the human homolog for Oas1g. (B) Violin plots showing the single-cell expression profiles for selected cluster marker genes in all the mouse pancreatic islet cell types identified by Baron et al. (Baron et al. 2016 [28])). (C) Western blot analysis for the anti-Ifit3 antibody using protein extract from mouse spleen. Gel image has been cropped to keep only relevant lanes. The molecular weight marker ladder and the Ifit3 western blot bands are part of the same image. Images have been taken using a LI-COR Odyssey equipment. Full-length blots/gels are presented as a supplementary file associated with Figure S3C. Figure S4. The presence of CD45+ cells is increased in islets after acute β-cell loss. (A). Immunofluorescence staining of pancreas sections obtained from Ctrl and 5dpDT treated mice. Control and 5dpDT bottom rows present zoom ins of regions indicated by white squares in the islets shown above. (B). Percentage of cells co-expressing Cd45 and YFP, expressed relative to the total number of total YFP+ cells in the islets. (C). Percentage of cells expressing Ifit3+ Cd45, and co-expressing both markers, expressed relative to the total number of Ifit3+ cells in the islets. Analysis focused on 5dpDT samples. Data for control samples is presented as reference. Values in panels (B) and (C) are expressed as the mean of 2 biological replicate experiments, with each value indicated by dots. Total number of islet cells counted per replicate: Control (replicate 1: 317, replicate 2: 76), 5dpDT (replicate 1: 368, replicate 2: 1236), 15dpDT (replicate 1: 126, replicate 2: 391). Figure S5. Ifit3 is also expressed in β-cells still present in 5dpDT islets. (A). Immunofluorescence images of pancreas sections obtained from Ctrl and 5dpDT treated mice showing Ifit3, Insulin and Glucagon co-staining. (B). Immunofluorescence images of pancreas sections obtained from 5dpDT treated mice showing YFP, Ifit3 and Glucagon co-staining. Figure S6. Il2r-g expression is strongly induced in α-cells after acute β-cell loss. (A). Representative immunofluorescence staining of pancreas sections obtained from Ctrl and 15dpDT treated mice. Note that Il2r-g is a protein mostly localized at the cell membrane, while YFP can present both nuclear and/or cytoplasmic stainings. Top right panels in 15dpDT images represent zoom ins of indicated regions. (B). Il2r-g gene expression in Control, 5dpDT, 15dpDT and 30dpDT as profiled by RNA-seq.

Additional file 1: Fig. S1 Total gene detection of all datasets compared after processing with either kallisto or Cell Ranger. The Venn diagrams show commonly detected number of genes by both pipelines and uniquely detected genes. Fig. S2 Violin-plots showing distribution of gene and UMI detection per cell of all the analyzed datasets (Table 1) run with the Cell Ranger pipeline. Fig. S3 Violin-plots showing distribution of gene and UMI detection per cell of all the analyzed datasets (Table 1) run with the kallisto pipeline. Fig. S4 Cell counts of all datasets compared after processing with either kallisto forced or Cell Ranger. The Venn diagrams show commonly detected cell barcodes by both pipelines and uniquely detected cell barcodes. Fig. S5 Alignment results of all datasets (Table 1) run with either Cell Ranger or kallisto forced against Ensembl reference. a Percent alignment rates of all reads against the reference transcriptome. b Total gene detection. c Median gene counts over all cells per dataset. d Median UMI counts over all cells per dataset. e Total cell counts of each dataset. Fig. S6 Total gene detection of all datasets compared after processing with either kallisto forced or Cell Ranger. The Venn diagrams show commonly detected number of genes by both pipelines and uniquely detected genes. Fig. S7 Violin-plots showing distribution of gene and UMI detection per cell of all the analyzed datasets (Table 1) run with the kallisto forced pipeline. Fig. S8 Violin-plots showing distribution of gene and UMI detection per cell of the dr_pineal_s2 dataset after additional filtering for downstream analysis. Run with either Cell Ranger (a), kallisto (b) or kallisto forced (c). Fig. S9 Downstream analysis of dr_pineal_s2 before cluster merging. a 2D visualization using UMAP of Cell Ranger analyzed clusters before merging, with resolution equal to 0.9. Each point represents a single cell, colored according to cell type. The cells were clustered into 21 types. b Expression profile of marker genes according to cluster [7] of (a). Clusters 0, 1, 8 and 18 are all rod-like PhRs subclusters. They expressed rod-like PhR markers (exorh, gant1, gngt1), but the expression levels differed and resulted in their separation. For simplicity, they were merged and referred as a single rod-like PhRs cluster in the main text. Similarly, cluster 7 and 12 were merged into a single Müller-glia like cluster, clusters 2, 5, 16 were merged into a single RPE-like cluster, clusters 3 and 10 were merged into a single habenula kiss1 cluster and cluster 11 and 19 were merged into a single leukocytes cluster. c. 2D visualization using UMAP of Cell Ranger analyzed clusters, with resolution equal to 2. The cells were clustered into 31 types. However, the two different cone-like PhR cell types are still not distinguished from one another. d Expression profile of marker genes according to cluster of (c). e 2D visualization using UMAP of kallisto analyzed dr_pineal_s2 clusters before merging, with resolution equal to 0.9. The cells were clustered into 24 types. f Expression profile of marker genes according to cluster of (c). Similar to the descried above, clusters 1, 2, 3, 7 and 21 were merged into a single rod-like PhRs cluster, clusters 0, 9, 17 were merged into a single RPE-like cluster, clusters 11 and 12 were merged into a single Müller-glia like cluster, clusters 4, 5 and 20 were merged into a single habenula kiss1 cluster and clusters 13 and 22 were merged into a single leukocytes cluster. g 2D visualization using UMAP of kallisto forced analyzed dr_pineal_s2 clusters, with resolution equal to 1.2. The cells were clustered into 27 types. h Expression profile of marker genes according to cluster of (g). The col14a1b gene was only detected in the kallisto and kallisto forced datasets and is the strongest DE marker within the red cone-like cluster (f, h). Fig. S10 Heatmap of genes with higher counts in kallisto pre-processed pineal data. All the UMI counts for both kallisto and Cell Ranger were summed, and the diff_ratio value was calculated ( kallisto _ counts − CellRanger _ counts kallisto _ counts + CellRanger _ counts \(\frac{\left( kallisto\_ counts- CellRanger\_ counts\right)}{\left( kallisto\_ counts+ CellRanger\_ counts\right)}\) ) for each gene (Additional file 1: Fig. 10). The top 80 diff_ratio genes, as well as the top 20 genes uniquely identified in kallisto were plotted according to the average scaled expression per cluster. Fig. S11 Heatmap of genes with higher counts in Cell Ranger pre-processed pineal data. All the UMI counts for both kallisto and Cell Ranger were summed, and the diff_ratio value was calculated ( kallisto _ counts − CellRanger _ counts kallisto _ counts + CellRanger _ counts \(\frac{\left( kallisto\_ counts- CellRanger\_ counts\right)}{\left( kallisto\_ counts+ CellRanger\_ counts\right)}\) ) for each gene (Additional file 1: Fig. S11). The top 80 diff_ratio genes, as well as the top 20 genes uniquely identified in Cell Ranger were plotted according to the average scaled expression per cluster.

Additional file 1: Figure S1. Study flow chart. Single-cell RNA-sequencing. Figure S2 a-g. Enriched functions of FBs, Eps, ECs, T cells, NK cells, M cells, ast cells and neutrophils by GO and KEGG analyses. Figure S3. Single-cell RNA-seq analysis revealed distinct characterization of cell populations in endometriosis lesions, eutopic endometrium and normal endometrium. (A) Similarity analyses of each cell population in endometriosis lesions, eutopic endometrium and normal endometrium. (B) Differential expression analysis was performed comparing whole cells from endometriosis lesions, eutopic endometrium and normal endometrium. (C) Representative significantly enriched GO and KEGG processes were performed with the significantly upregulated genes in endometriosis lesions, eutopic endometrium and normal endometrium. (D) Expression of selected pathway genes were shown for each cell of endometriosis, eutopic endometrium, or normal endometrium origin. Dot size corresponded to the percentage of cells in the cluster expressing a gene, and dot color corresponded to the average expression level for the gene in the cluster. Figure S4. DEG analysis compared cells from endometriosis lesions, eutopic endometrium, and normal endometrium within ECs and Eps by heatmap. Figure S5. Immunofluorescence in cells showed the C3, C7, StAR and S100A10 positive FBs in endometriosis lesions. Figure S6. (A) Contribution of each group to each cell state on Figure 2E. The majority of state 1 was occupied by FBs from the eutopic endometrium and normal endometrium. State 2 was primarily contained by ectopic and eutopic endometrium FBs, while state 3 contained FBs from all three groups. (B) Relative contributions of FBs from13 subclusters to each group. (C) Relative contributions of FBs from three groups to each cluster, as shown by the t-SNE plot. (D)Western blot analysis reflected si-StAR transfection efficiency. (E) Scratch wound healing determined the ability of migration between si-Ctrl and si- StAR FBs. (F) Transwell assays were used to compare the effects of si- StAR on migration and invasion. (G) Proliferative ability between si-Ctrl and si- StAR FBs. (H) Representative images of flow cytometric cell cycle analysis. (I) Apoptosis of si-Ctrl and si- StAR FBs by flow cytometric. Figure S7. (A) Bar plot demonstrated the relative ratio of each subcluster to the entire T cell population. (B) Relative contributions of T cells from three groups to each cluster. (C) Contribution of each CD4+ T cells group to each cell state on Figure 6B. The majority state 1 was occupied by CD4+ T cells from eutopic endometrium and normal endometrium. State 2 contained CD4+ T cells from all three groups while half of CD4+ T cells in state 3 were from endometriosis lesions. (D) Contribution of each CD8+ T cells group to each cell state on Figure 6F. The majority of state 1 was occupied by CD8+ T cells from endometriosis lesions. Eutopic endometrium and normal endometrium and half of CD8+ T cells in state 2 were from endometriosis lesions. State 3 and state 4 contained T cells from all three groups while normal endometrium CD8+ T cells accounted for most of state 4. State 5 was mostly occupied by CD8+ T cells from eutopic and normal endometrium. (E) Bar plot demonstrated the relative ratio of each subcluster to three NK cell groups. (F) Relative contributions of T cells from three groups to each cluster. (G) Bar plot demonstrated the relative ratio of each subcluster to three M cell groups. (H) Relative contributions of M cells from three groups to each cluster. (I) Heatmap showed average expression of specific genes of M1 and M2 including CD64, CD40, CD86, CD163 and CD206. Figure S8. Functional enrichment analysis with GO and KEGG analyses of M-1, M-3, M-4, M-5, M-6, M-7, and M-9. Figure S9. The dense network and multiple cellular connection in eutopic endometrium and normal endometrium. (A) Putative signaling between differentially expressed receptors in different cell types and their ligands. Compartments represented cell types, and their preferentially expressed receptors and ligands were labeled along the outer margin. The solid line indicated that there was a significant interaction, and the dashed line indicated that it was not significant. (B) Capacity for intercellular communication between FBs and immune cells. Each line color indicated the ligands expressed by the cell population represented in the same color (labeled). The lines connected to the cell types that expressed the cognate receptors. Figure S10. (A)Heatmap of average gene expression of ESR2, PGR, StAR, CYP19A1. (B) Relative contributions of each cluster to each sample. (C) t-SNE plots of cells from each sample profiled in this study, with each cell color coded to indicate the associated cell types.

Table S6. List of ligand-receptor pairs and cell–cell pairs used in Fig. 4c for heatmap. (XLSX 12 kb)

A. Cell type labels used based on re-analysis of IPF and healthy control data from GSE135893 (Kropski-Vanderbilt Univ single cell cohort) [24]. Clustering was performed using R package Seurat and cell types were identified using known markers. Ciliated_0 and Ciliated_1: Ciliated epithelial cell subpopulations; AT2_2, AT2_13, AT2_29, AT2_30: Alveolar epithelial cell type II subpopulations; SPP1_mac_3: SPP1+ monocytes/macrophages; C1QA_mac_4, C1QA_mac_5, C1QA_mac_9, C1QA_mac_12: C1QA+ macrophage subpopulations; Mono_7, Mono_21: Monocyte subpopulations; Tc_8: cytotoxic T cells; Th_10: helper T cells; AT1_11, MUC5Bpos_AT1_15, Basal_AT1_17: Alveolar epithelial cell type I subpopulations; ACKR1_pos_endo_14: ACKR1+ endothelial cells; ACKR1_neg_endo_16 and ACKR1_neg_endo_20: ACKR1- endothelial cell subpopulations; Diff_cil_18: Differentiating ciliated epithelial cells; Fibroblasts_19 and Fibroblasts_23: Fibroblast subpopulations; Sm_26: smooth muscle; Prolif_mac_22: Proliferating macrophages; Ly_endo_24: Lymphatic endothelium; Bcells_25: B cells; PC_28: Plasma cells; MC_27: mast cells; Mesothelial_31: mesothelial cells. B. Heatmap (left panel) and correlation matrix (right panel) in GSE47460 (Kaminski-LGRC bulk expression cohort) of genes included in the signature derived from the dataset shown in panel A. (ZIP)

Clustering was performed using R package Seurat and cell types were identified using known markers. A. Total lung cell suspension. SPP1_monocytes_0: SPP1+ monocytes; Infl_monocytes_1: Inflammatory monocytes; ACKR1pos_endo_2: ACKR1+ endothelial cells; ACKR1neg_endo_3: ACKR1- endothelial cells; Fibroblasts_4: Fibroblasts; AT2_5 and AT2_23: Alveolar epithelial cell type II subpopulations; Th_6: helper T cells; Pericytes_7 and Pericytes_22: Pericyte subpopulations; HLAhigh_mac_8 and HLAhigh_mac_10: HLA class II high macrophage subpopulations; Sm_9: smooth muscle cells; Bcells_11 and Bcells_21: B cell subpopulations; Tc_12: cytotoxic T cells; AT1_13: Alveolar epithelial cell type I; PC_14: Plasma cells; Endo_15 and Endo_24: endothelial cell subpopulations; Ciliated_16: ciliated epithelial cells; Monocytes_17 and Monocytes_18: Monocyte subpopulations. B. Lineage sorted cells. THY1high_alv_fib_0: THY1 high alveolar fibroblasts; THY1pos_sm_1: THY1+ smooth muscle; THY1neg_sm_2: THY1- smooth muscle; CTHRC1pos_3: CTHRC1+ fibroblasts; Adventitial_4: Adventitial fibroblasts; THY1neg_alv_fib_5: THY1- alveolar fibroblasts; Pericytes_6: Pericytes; Peribronchial_7: Peribronchial fibroblasts; Sm_8 and Sm_13: smooth muscle cell subpopulations; Alveolar_9 and Alveolar_10: Alveolar fibroblast subpopulations; Epi_11: Epithelial cells; Hematopoietic_12 and Hematopoietic_14: Hematopoietic cells. C. Heatmap (left panel) and correlation matrix (right panel) in GSE47460 of genes included in the signature derived from the ‘Total lung cell suspension’ (shown in panel A) dataset across each cluster shown in panel A. D. Heatmap (left panel) and correlation matrix (right panel) in GSE47460 of genes included in the signature derived from the ‘Lineage sorted’ (shown in panel B) dataset across each cluster shown in panel B. (ZIP)