Single-cell barcoding, library preparation, and sequencing

Single cells were barcoded using the 10× chromium single-cell platform, and complementary DNA (cDNA) libraries were prepared according to the manufacturer’s protocol (Single Cell 3′ Reagent Kits v2, 10× Genomics, USA). Cell suspensions, reverse transcription master mix, and partitioning oil were loaded on a single-cell “A” chip with a targeted cell output of 10,000 cells per library, assuming 100% viability, and then run on the Chromium Controller. Reverse transcription was performed within the droplets at 53°C for 45 min. cDNA was amplified for 12 cycles total on a Bio-Rad C1000 Touch thermocycler. cDNA size selection was performed using SPRIselect beads (Beckman Coulter, USA) and a ratio of SPRIselect reagent volume to a sample volume of 0.6. cDNA was analyzed on an Agilent Bioanalyzer High Sensitivity DNA chip for qualitative control purposes. cDNA was fragmented using the proprietary fragmentation enzyme blend for 5 min at 32°C, followed by end-repair and A-tailing at 65°C for 30 min. cDNA were double-sided size selected using SPRIselect beads. Sequencing adaptors were ligated to the cDNA at 20°C for 15 min. cDNA was amplified using a sample-specific index oligo as a primer, followed by another round of double-sided size selection using SPRIselect beads. Final libraries were analyzed on an Agilent Bioanalyzer High Sensitivity DNA chip for qualitative control purposes. cDNA libraries were sequenced on a HiSeq 4000 Illumina platform aiming for 150 million reads per library and a sequencing configuration of 26 base pair (bp) on read1 and 98 bp on read2.

Fastq generation and read trimming

Basecalls were converted to reads with the software Cell Ranger’s (v2.2) implementation mkfastq. Multiple fastq files from the same library and strand were catenated to single files. Read2 files were subject to two passes of contaminant trimming with cutadapt (v1.17): (i) for the template switch oligo sequence (AAGCAGTGGTATCAACGCAGAGTACATGGG) anchored on the 5′ end and (ii) for poly(A) sequences on the 3′ end. Following trimming, read pairs were removed if the read2 was trimmed below 30 bp.

Cell barcode and unique molecular identifier demultiplexing, alignment, and transcript counting

Subsequent read processing was conducted with the zUMIs pipeline (v2.0). Paired reads were filtered if either the cell barcode or unique molecular identifier (UMI) sequence had more than 1 bp with a phred of <20. Reads were aligned with STAR (v2.6.0c) to the human genome reference GRCh38 release 91 from ensemble. Collapsed UMIs with reads that span both exonic and intronic sequences were retained as both separate and combined gene expression assays.

Filtering valid cell barcodes and quality control

Cell barcodes representative of quality cells were delineated from barcodes of apoptotic cells or background RNA based on a threshold of having at least 12% of transcripts arising from intron spanning, or unspliced reads, indicative of nascent mRNA (fig. S8). Cells with less than 1000 transcripts profiled or more than 20% of their transcriptome of mitochondrial origin were then removed. The average number of cells returned per library is 2925 (±1811.852). Sequencing depth comparisons across samples categorized by associated disease state were statistically assessed for disease-associated biases via Kruskal-Wallis test, where differences were observed to be statistically insignificant ( P value of 0.29; fig. S9A). A detailed summary of the average UMI and gene per cell type per library distribution per disease state is presented in fig. S9 (B and C, respectively).

Gene name conversion

Genes were originally output in ensemble gene ID format. To improve the interpretability of the variables without making compromises to sensitivity, we converted ensemble gene IDs to Hugo Gene Nomenclature Committee (HGNC) format using the R package BioMart only when an exact one-to-one translation was available.

Data normalization and cell population identification

UMIs from each cell barcode, irrespective of intron or exon coverage, were retained for all downstream analysis. Raw UMI counts were normalized with a scale factor of 10,000 UMIs per cell and subsequently natural log transformed with a pseudocount of 1. Aggregated data were subject to Louvain cluster analysis for cell type identification using the R package Seurat (version 2.3.1). Recursive clustering analysis of subpopulations of pure immune (PTPRC ⁺ ), epithelial and the remaining mesenchymal populations were conducted to improve the granularity of our cell annotations.

Multiplet cell populations were clusters, which were identified as having a transcriptomic signature that resembled the resulting combination of two or more disparate cell type signatures that already existed in the dataset. For each cell population flagged as potential multiplets, a differential expression test was performed comparing the population against their suspected combination of cell populations; a confident assessment of multiplet status was made when there were no genes uniquely expressed in the suspected population and no genes represented in the cell combination group that could not be detected in the suspected multiplet population. The average number of multiplets identified per library is 49.33% (±63.8) or 1.29% (±0.95%) of all cells identified in each library. Cell barcodes identified as multiplets were not included in downstream analyses.

Generation of cell type markers

Cell type marker lists were generated with two separate approaches. In the first approach, we sought to adjust for unbalanced representation of different cell types across subjects and reduce the zero inflation in the data to generate a more reliable P value from a differential expression test. We collapsed the gene expression matrix values to one where each column represents average gene expression for a given subject and for a given cell type. To highlight gene expression among relatively similar cell types, we then grouped each cell type into one of four groups: epithelial, stromal/mesenchymal, myeloid, and lymphoid. We then applied the Seurat FindAllMarkers implementation of the Wilcoxon rank sum test for each of these four subsets of cell types to generate separate marker lists for each (data S2 to S5).

The second approach to generating cell type markers uses a binary classifier system to assess the utility of detecting a given gene, irrespective of its intensity of expression, for classifying a cell. For each cell type in the data, we identified the genes whose expression was log fold change ≥ 0.3 greater than the other cells in the data. We then calculated the diagnostics odds ratio (DOR) for each of these genes, where we binarize the expression values by treating any detection of a gene (normalized expression value > 0) as a positive value and zero expression detection as negative. We included a pseudocount of 0.5 to avoid undefined values as

where TruePositives represents the number of cells within the group detected expressing the gene (value > 0), FalsePositives represents the number of cells outside of the group detected expressing the gene, FalseNegatives represents the number of cells within the group with no detected expression, and TrueNegatives represents the number of cells outside of the group with no detected expression of the gene. The log transformed DOR marker values are contained in data S6.

Cell population composition comparisons

All figures and statistical tests for proportion makeups were performed exclusively with subjects who had nonzero proportion makeups of the given cell type. This choice was made in an effort to be conservative in our estimates, as we typically considered the lack of detection of a given cell type as a technical limitation. Test results reported in data S12 are Wilcoxon rank sum tests of nonzero proportions between IPF and control subjects. We advise caution in the interpretation of these statistics as the model does not account for unbalanced sampling rates of cells across subjects, and the proportion values of any given cell type are not independent of the proportions for any other cell types that belong to the same category.

Differential expression between disease conditions

A major limitation that persists in scRNA-seq analysis is the lack of a differential expression model that can account for cellular population differences while ultimately testing differences at the level of the subject, rather than treating each individual cell as an independent sample. While we were unable to design a model capable of handling the complex sample hierarchy scRNA-seq, we nonetheless applied a crude method of differential expression.

For each subject and each cell type, the gene expression for each gene was averaged to create a single “sample” representative of an individual. We then applied a Wilcoxon rank sum test between within each cell type and average values from subjects of each disease condition for a total of three comparisons (IPF versus control, COPD versus control, and IPF versus COPD; data S8, S9, and S10, respectively). Cell types with less than five subjects representative of a particular disease condition were not included its analysis.

UMAP visualizations

For similarity-based cell network analysis and visualization, we used tools from the Python (version 3.6.8) library Scanpy (version 1.3.7). Uniform Manifold Approximation and Projection (UMAP) figures for all subsets of cells were generated using the same sequence of implementations. Feature selection of the top genes ranked by dispersion (scaled variance/mean) across 20 bins of the expression distribution are identified. The expression values for these genes are then adjusted for differences in total UMI and the fraction of mitochondrial reads across cells during z normalization with a maximum absolute z score of 10. The scaled values are then subject to principal components analysis (PCA) for linear dimension reduction.

The top principal components are subject to exploratory analysis to identify contributions to variance at the level of library batch, subject, cell type, and disease. In addition, the residuals of each principal component were explored to ascertain functional relevance of the signatures. Feature selection of principal components was conducted on the basis of these analyses, and a shared nearest neighbor network was then created on the basis of Euclidean distances between cells in multidimensional PC space and a fixed number of neighbors per cell, which was used to generate an intermediate two-dimensional (2D) UMAP.

The resulting network was then subject to PAGA using the cell type categories as abstraction nodes. A confidence threshold of cell type interconnectivity was implemented to avoid spurious manifold connections between cell types of disparate lineages. Diffusion maps is then applied for nonlinear dimension reduction on a more limited subset of principal component dimensions, and a new neighborhood graph is then computed on the basis of representations of diffusion map distances. A final 2D UMAP is then generated using the new neighborhood graph distances, initialized by the positions from the thresholded PAGA.

The values for all nondefault parameters used during the generation of each UMAP are represented in data S11; if a parameter is unspecified, then the default Scanpy implementation parameter was used.

Evaluation for the potential influence of outlying subject-specific variation

To ensure that our single-cell analysis was not unduly driven by the outlying characteristics of one or few patient transcriptional profiles, we ran a parallel analysis using deep generative modeling to “correct” for any potentially aberrant subject-specific variational signatures in accordance with the single-cell variational inference (scVI) approach described by Lopez et al . ( 40 ). This was implemented through a dedicated python library downloaded from GitHub ( https://github.com/YosefLab/scVI ). We found that, even after applying this “batch effect” normalization, our cluster groups were preserved with high fidelity as evaluated by multiple metrics (fig. S1), reinforcing that outlying subject-specific effects were not driving the results of our analysis.

Evaluation for the potential influence of cell cycle state

To further demonstrate the robustness of our approach, we applied a similar method to evaluate for possible data artifacts specific to cell cycle state. In another parallel analysis, we used the approach described by Scialdone et al . ( 41 ) to extract G ₁ , G ₂ -M, and S components from each individual cell before normalization using a set of 94 known cycle-associated genes. We then regressed out for these elements using Scanpy, similar to our approach for mitochondrial RNA, and proceeded with downstream analysis. We found that the G ₁ , G ₂ -M, and S states were well distributed at both a cluster level and patient level. Furthermore, regressing out the effects of cell cycle state did not meaningfully alter our resulting cluster configurations (fig S4).

Comparison of manual cell annotations with automated methods

We compared our manual annotations to those produced through automated classification using SingleR. Strong correlation with manual annotations was found using both the Human Primary Cell Atlas (HPCA) and ENCODE human reference sets across our dataset comparable to the correlation between the HPCA and ENCODE annotations themselves for these cells (fig. S3).

PAGA connectivity analysis of fibroblast and myofibroblast

To assess the most likely trajectories of cell progression toward IPF-enriched fibrosis among fibroblast and myofibroblast, we used unsupervised Louvain clustering to generate eight subpopulations, which were then subjected to PAGA analysis to ascertain the most likely intersubcluster trajectories. The edge confidences between each subcluster node for all edges is visualized using the R package igraph (v1.2.4.1).

Scoring of regulon activity

A regulon is a group of target genes regulated by a common transcription factor. To score the activity of each regulon in each nonimmune cell, we used the package pySCENIC with default settings and the following database:

cisTarget databases (hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather, hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather), and the transcription factor motif annotation database (motifs-v9-nr.hgnc-m0.001-o0.0.tbl) were downloaded from resources.aertslab.org/cistarget /, and the list of human transcription factors (hs_hgnc_tfs.txt) was downloaded from github.com/aertslab/pySCENIC/tree/master/resources .

Archetype analysis of fibroblast and myofibroblast

We observed that many IPF-enriched features in the data were represented by a continuum of increasing phenotypic deviation from controls, rather than discrete features readily amenable to delineation (e.g., cluster analysis) from control-enriched signals. Consequently, we sought to implement archetype analysis of these continua to assess disease-enriched features rather than relying on traditional group-wise comparisons.

We first assessed the most likely trajectories toward fibrotic archetypes among fibroblast and myofibroblast using the DPT implementation from Scanpy to plot the distances along the UMAP manifold toward each archetype’s terminus. The same diffusion map component structure used to generate the 2D UMAP visualizations was used for calculating DPT distances. For fibroblast, we took the cell with the highest expression of ITGB1, which lied at the terminus of fibroblast manifold’s tendril as the root cell and calculated the relative distances of all other fibroblast cells. We used 1-DPT value for the final distance value ordering and assigned a gradient of colors to cells for the range of 0 to 0.7 to represent each cells distance on the UMAP legend and in the accompanying heat map.

For myofibroblast, we calculated DPT ordering from all three termini in the myofibroblast manifold using the highest MMP11-expressing cell to represent the IPF-enriched terminus, the furthest control cell to represent the control, and COPD-enriched terminus and the furthest cell from 225I to represent the single-subject–enriched terminus. To mitigate the impact on our analysis from the single-subject–enriched archetype, we removed the nearest 500 cells to the subject “225I-” enriched terminus that did not overlap with the nearest 399 cells to the control and COPD-enriched myofibroblast archetype. Among the remaining myofibroblast cells, we used the difference between both 1-DPT distance values from the remaining two nonsubject-specific archetypes as a singular distance vector. A color gradient for the final 1-DPT range of −1 to 1 was then applied to each unfiltered cell in the UMAP legend and accompanying heat map.

For both archetype heat maps, cells were plotted in order of the final 1-DPT values from lowest to highest. A spearman correlation test was conducted between each cell’s 1-DPT value and gene expression to ascertain which genes increased or decreased in expression along the manifold and would thus serve as candidate features in the heat map. In addition to each cell’s DPT distance value associated color, each cell’s respective disease and subject color identity were included in the annotation bar to represent the extent of disease enrichment and subject-level contribution to the feature.

Pseudotime orderings were each subject to statistical assessment to determine whether cells from IPF subjects were enriched at the far end of the DPT ordering results (the rightmost portion of Fig. 4, D and E , heat maps). The first approach was a Wilcoxon rank sum test (performed with the base R “wilcox.test” implementation) between IPF and control for the average DPT distance per subject, with the alternative hypothesis that IPF DPT distribution was right skewed. The second approach was a generalized linear mixed model (performed with the R package nlme v3.1’s implementation) applied to the DPT distance of the cells, where disease was a fixed variable and subject was treated as a random variable nested within each disease.

Archetype analysis of classical monocytes and macrophage

Among classical monocytes and macrophage, we identified one monocyte archetype connected to distinct four macrophage archetypes: an inflammatory archetype, an IPF-enriched archetype, an outlier mitochondrial transfer RNA (MT-tRNA)–enriched archetype driven by two subjects, and another outlier archetype driven by two separate subjects characterized by heat shock protein expression. 1-DPT distance values for all monocytes and macrophage were calculated from each archetype terminus as described for fibroblast and myofibroblast. We removed 1700 and 5800 most proximal cells to the MT-tRNA and heat shock protein archetypes, respectively, to avoid contributions from outlier signals.

Following outlier removal, the 1-DPT distances values from the monocyte, inflamed macrophage and IPF-enriched macrophage were each independently unity normalized to values between 0 and 1. Distances along the three normalized trajectories were then used plot the cells in a ternary plot using the R gglot2 extension package tricolore (v1.2.0) to assign a color to each cell based on its relative proximity to each of the three archetype termini, which is visualized in both the UMAP color legend and ternary plot legend itself. The 15,000, 20,000, and 25,000, closest cells to the monocyte, inflamed macrophage and IPF-enriched macrophage archetype, respectively, were then selected for correlation analysis between 1-DPT values from each respective trajectory and gene expression to assess candidate genes for plotting the heat map.

GRN construction

GRNs for both control and IPF cell populations were generated using the R package bigSCale2 ( 24 ). Control and IPF cells were split, and for each disease state, cells were randomly down-sampled using Seurat’s SubsetData implementation (seed = 7) to a maximum of 500 cells per cell type, resulting in 10,560 cells from control and 15,388 cells from IPF. The resulting matrices were then filtered to remove genes with ensembl identifiers and passed to bigSCale2, where both networks were constructed under the “normal” clustering parameter, with an edge cutoff of the top 0.993 quantile for correlation coefficient. Networks were visualized with the R package igraph, and each network’s layout is derived from 10,000 iterations of the Fruchterman-Reingold algorithm and “nogrid” parameter (seed = 7).

Assigning cell relevance scores to GRN communities

Each GRN was clustered with igraph’s “cluster_fast_greedy” function (seed = 7). Similar cell types were collapsed into the groups to avoid excessive phenotypic overlap: all dendritic cells, both macrophage varieties, both monocyte varieties, all T cells, natural killer and innate lymphocytes, both B cell varieties, both secretory, and pulmonary neuroendocrine cell and ionocytes. For each cell type group, we then took the top 500 genes ranked by log(DOR) markers as representative features. The relevance score for each cell type group in a network’s community is a z score of the cumulative log(DOR) markers that intersect with members of each GRN community, weighted by the gene’s normalized PageRank centrality within its community as follows

where PR denotes the PageRank centrality for a node in the network and DOR is the diagnostic odds ratio classifier calculated earlier in the cell type marker list. Given i genes, j communities, and k cell type groups, let PR ij * denote the normalized PageRank value for each gene i ( i = 1, …, I) within its respective network community j ( j = 1, …, J), let log(DOR) _ik be the natural log transformed DOR for gene i that belongs to cell type group k ( k = 1, …, K), let i ∈ kj represent the intersection of genes from cell type group k and GRN community j , and let CumulativeScore _kj represents the cumulative relevance score for cell type group k in GRN community j . Last, the cumulative scores are then converted to z scores by centering and scaling across cell type groups within each community, only cell type groups with a relevance score greater than or equal to 1 were retained for annotation of the community. The aberrant basaloid cell type community was excluded from the control GRN analysis because this population of cells was not detected in any control samples.

Immunohistochemistry

The FFPE (Formalin fixed paraffin embedded) blocks were cut at 5 μm and then rehydrated using standard xylene/ethanol deparaffinization. For heat-induced antigen retrieval, specimens were boiled at 95°C for 20 min in 1× tris-based Antigen Unmasking Solution (Vector Laboratories, USA). Tissue slides were incubated for 10 min in BLOXALL Endogenous Peroxidase and Alkaline Phosphatase Blocking Solution (Vector Laboratories, USA) to block endogenous peroxidase and alkaline phosphatase activity. Tissue slides were blocked using 2.5% Normal Horse Serum Blocking Solution (Vector Laboratories, USA) for 20 min, then incubated with the primary antibody (data S13), and diluted in 2.5% Normal Horse Serum Blocking Solution for 30 min at room temperature. Specimen were incubated for 30 min with secondary antibodies (anti-mouse or anti-rabbit ImmPRESS reagent, conjugated with horseradish peroxidase or alkaline phosphatase, as appropriate; all Vector Laboratories, USA). Slides were incubated for 10 min in DAB working solution or for 30 min in Vector Reds working solution, as appropriate (both Vector Laboratories, USA). For sequential double stainings, this protocol was repeated from the protein-blocking step on using the alternative enzyme and substrate reagents. Tissue slides were counterstained in Hematoxilin Solution Gill no. 1 (Sigma-Aldrich, USA) for 3 min and then washed with tap water. For permanent mounting, specimens were dehydrated in ethanol/xylene and then mounted with VectaMount permanent mounting solution (Vector Laboratories, USA). Stained slides were digitalized on a Aperio Scanner (Leica) and then analyzed using the softwares QuPath and ImageJ.

Independent validation of COL15A1 protein expression in pVE cells

Images of COL15A1 immunohistochemical stainings of lung parenchyma and bronchi specimens (fig. S7) were downloaded from the Human Protein Atlas [( 19 ); www.proteinatlas.org/ENSG00000204291-COL15A1/tissue ). All specimens had been stained using the polyclonal antibody HPA017915 (Sigma-Aldrich).

We are indebted to all patients and control subjects who participated in this study. The sequencing was conducted by M. Zhong at Yale Stem Cell Center Genomics Core facility. We thank A. Brooks and the staff of Yale Pathology Tissue Services for tissue processing. We also thank M.C. Nawijn for providing access to group’s data. Funding: This work was supported by NIH grants R01HL127349, U01HL145567, U01HL122626, and U54HG008540 to N.K., NHLBI P01 HL114501 and support from the Pulmonary Fibrosis Fund to I.O.R., an unrestricted gift from Three Lake Partners to I.O.R. and N.K., and by the German Research Foundation (SCHU 3147/1) to J.C.S. Author contributions: N.K. and I.O.R. conceptualized, acquired funding, and supervised the study. B.A.R. and G.R.W. performed sample collection and phenotyping. S.P., E.A.A, and S.G.C. procured and dissociated the lungs. T.S.A., J.C.S., S.P., F.A., and G.D. performed library construction. Data were processed, curated, and visualized by T.S.A., J.C.S., N.N., and X.Y. and analyzed by T.S.A., J.C.S., N.K., X.Y., N.N., M.J., Q.D., H.A.A., and A.S. The online tool “ IPFCellAtlas.com ” was developed by N.N. IHC was performed by J.C.S. and T.S.A and evaluated by R.H., N.K., and J.C.S. The manuscript was drafted by T.S.A, J.C.S, and N.K. and was reviewed and edited by all other authors. Competing interests: T.S.A., J.C.S., S.P., E.A.A, H.A.A., A.S., I.O.R., and N.K. are inventors on a provisional patent application (62/849,644) submitted by NuMedii Inc., Yale University, and Brigham and Women’s Hospital Inc. that covers methods related to IPF-associated cell subsets. N.K. served as a consultant to Biogen Idec, Boehringer Ingelheim, Third Rock, Pliant, Samumed, NuMedii, Indaloo, Theravance, LifeMax, Three Lake Partners, and Optikira over the past 3 years and received nonfinancial support from MiRagen. H.A.A., A.S., M.J., and Q.D. are employees of NuMedii Inc. The authors declare that they have no other competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Transcriptomic data was deposited to the GEO under accession number {"type":"entrez-geo","attrs":{"text":"GSE136831","term_id":"136831"}} GSE136831 .