Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing

doi:10.21203/rs.3.rs-1716190/v1

Download PDF

Research Article

Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing

https://doi.org/10.21203/rs.3.rs-1716190/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

The progression of multiple primary lung cancer (MPLC) involves complex changes in cell composition and metabolic function. Here, we performed scRNA sequencing of 167,397 cells from six patients with MPLC, combined with bulk whole-exome sequencing. We revealed that both naïve and memory T cells participate in the differentiation of CD8+ T cells. The terminal states of CD8+ T cells are exhausted T cells, which respond to stimuli and positively regulate cell death, and cytotoxic T cells, which are mainly implicated in the regulation of cytokine production. Multiple functional Tregs and naïve T cells contribute to the origin of CD4+ T cells. B cells, with two main functions, also play important roles in the immune response. We also uncovered the cellular metabolic activities that occur during tumor invasion. Positive regulation of blood vessel diameter has been observed in endothelial cells, while angiogenesis has not been found in early glandular neoplasia of the lung. The development of epithelial cells involves two functional states. In one state, cells respond to stress from the immune system, and in the other, some will undergo programmed cell death or enter the cell cycle due to the selective pressure that arises over the course of tumor development in synchronous MPLC. Our study showed the complete landscape of different dynamic cellular changes, which might reveal the key cellular mechanisms of MPLC and therefore provide new clues for the pathogenesis of tumors.

multiple primary lung cancer

single-cell RNA sequencing

heterogeneity

tumor microenvironment

Lung cancer is the leading cause of cancer-related mortality worldwide ^1–3. Multiple primary lung cancer (MPLC) is defined as the occurrence of two or more primary lung cancers in the same patient. MPLC can be classified as synchronous or metachronous based on the time of occurrence; metachronous MPLC is defined by a diagnostic interval of 6 months or less between tumors^4,5. Port et al. reported that MPLC was observed in approximately 16% of patients with operable stage I, II, and III non-small-cell lung cancer (NSCLC), usually by preoperative imaging analysis⁶. Concomitant with the rising prevalence of lung cancer, the incidence of MPLC is also increasing ^4,7. Thus, an in-depth understanding of its mechanism and better treatment options for MPLC are urgently needed.

Tumor heterogeneity describes differences between tumors of the same apparent type but with different genetic backgrounds and different carcinogen exposures (intertumoral heterogeneity) or differences between cancer cells within the same tumor (intratumoral heterogeneity)⁸. MPLC, by definition, arises in different areas of the lung parenchyma within a single patient and shares an identical germline genetic background and environmental exposure profile⁹. Previous studies have demonstrated differences in genomic profiles, including somatic point mutations, copy number aberrations, chromosomal structural variations and even mutational spectra, between different MPLC lesions^9–11. The highly complex and heterogeneous tumor ecosystem contains not only malignant cells but also interacting cells from the host, such as endothelial cells, stromal fibroblasts, and a variety of immune cells that control tumor growth and metastasis¹². However, the heterogeneity of cell types in the tumor microenvironment (TME) in MPLC remains unclear. Thus, it is critically important to decode the complex interplay between tumor cells and the TME in MPLC.

In the present study, we performed whole-exome sequencing (WES) and single-cell RNA sequencing (scRNA-seq) on 17 tumor samples and 6 matched normal samples from 6 patients with MPLCs. By comparing tumors at different stages in MPLCs and adjacent normal lung tissue, we comprehensively characterized the transcriptome features of the malignant cells, immune cells, and stromal cells of MPLCs, and we decoded dynamic changes in cell percentage, cell subtype heterogeneity, and intercellular interactions, providing new knowledge regarding the biological basis of MPLC and LUAD development.

Patients and sample collection

Patients who were diagnosed with pathologically confirmed MPLC with two or three tumors in the ipsilateral or bilateral lungs and underwent surgery (surgery for contralateral lung tumors was completed within 1–2 months after the first surgery, depending on the patients’ recovery status) at the Department of Thoracic Surgery of the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China, were enrolled. After MPLC samples were resected, tumor tissues were cut into two pieces along the long axis: one was processed for pathological diagnosis, and the other was processed for scRNA-seq and WES. The clinical information and specimen information of these patients are summarized in supplementary table 1. This study was approved by the National Cancer Center/Cancer Hospital Ethics Committee (Approval number: NCC2208). Written informed consent was obtained from all participants included in this study.

Tissue dissociation.

For each sample, 100–200 mg of tissue was homogenized by mincing with scissors into smaller pieces (-0.3 mm² per piece). Before tissue enzyme digestion, blood and MACS Tissue Storage Solution were washed from the tissue pieces by the addition of 10 ml of ice-cold PBS, followed by further gentle washing and tissue collection using a 70 µm strainer. Then, the tissue pieces were transferred into 5 ml of enzyme mix prepared by combining 4.7 ml of DMEM, 200 µl of Enzyme H, 100 µl of enzyme R and 25 µl of Enzyme A. The samples were incubated at 37°C for 30–60 min in an Eppendorf thermomixer at 700 rpm. After incubation, the cell suspension was filtered through a prewetted 70 µm MACS SmartStrainer, placed in 50 ml tubes and centrifuged at 300 g for 5 min at 4°C. The cell pellet was washed with 5 ml chilled 1X PBS with 0.04% BSA and then centrifuged. One milliliter of chilled 1X red blood cell removal solution was added to the cell pellet, and the cells were resuspended by gentle pipetting with a wide-bore pipette tip and then incubated for 10 min on ice. After incubation, 10 ml chilled 1X PBS with 0.04% BSA buffer was immediately added to the sample tube. The tube was gently inverted several times to mix and centrifuged at 4°C at 300 g for 10 min. The cell pellet was resuspended in 1 ml chilled 1X PBS with 0.04% BSA buffer, and the cell concentration and viability were counted using a Countess II automated cell counter. If the percentage of viable cells was < 70%, dead cells were removed to increase the proportion of viable cells. An appropriate volume of buffer was added to the cell suspension and gently mixed to achieve a target cell concentration of 700–1200 cells/µl. Once the target cell concentration was obtained, GEM generation was immediately performed with 10X Genomics reagents. The remaining cells were frozen at -80°C for WES (no more than 5 × 10⁶ cells were used).

Whole-exome Sequencing.

We used a QIAamp DNA Mini Kit (Catalog# 51304) to process DNA extraction. The cell pellet was resuspended in PBS to a final volume of 200 µl, and 20 µl proteinase K was added to each sample and mixed by vortexing. Then, 4 µl of RNase A stock solution (100 mg/ml) was added to the sample and mixed completely, and 200 µl of Buffer AL was added to the sample and mixed by pulse-vortexing for 15 s. The samples were incubated at 56°C for 10 min, and then 200 µl of ethanol was added. The mixture was carefully transferred to a QIAamp Mini Spin column. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. The QIAamp Mini Column was placed in a clean 2 ml collection tube, and the tube containing the filtrate was discarded. The column was carefully opened, and 500 µl Buffer AW1 was added. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. This step was repeated with buffer AW2 and eluted with 100 µl buffer AE. One microliter of DNA sample was taken for quantification on a 1% agarose gel. After passing the QC, the SureSelect XT Library Prep Kit (Catalog#: 5500 − 0132) was used for library construction and QC according to the official protocol. Finally, we used an Illumina Nova 6000 for 150 bp paired end (PE 150) sequencing.

Bulk RNA sequencing.

We used the RNeasy MinElute Cleanup Kit (Catalog: 74204) to extract RNA. The cell pellet was resuspended in TRIzol reagent to a final volume of 1000 µl. Then, the samples were vortexed, and 200 µl chloroform was added to mix the samples. The lysate was centrifuged for 15 min at 12000 x g at 4°C. The supernatant was pipetted, and one volume of RLT buffer was added and mixed well. Then, one volume of ethanol was added and mixed well. The samples were transferred to a RNeasy MinElute spin column and centrifuged for 15 s at 8000 x g. The flowthrough was discarded, and the same steps were repeated with buffer RW1. Then, 80 µl DNase I solution was added to the samples and incubated for 8 min. After incubation, we added 350 µl of RW1 to the samples and centrifuged them at 8000 x g for 15 s. The flowthrough was discarded, and the wash steps were repeated with 500 µl of RPE buffer and 500 µl of freshly prepared 80% ethanol. Finally, the RNA samples were eluted in 30 µl RNase-free water, and 1 µl of the RNA sample was taken for quantification and evaluation of RNA integrity. After passing the QC, we used the TruSeq RNA Library Prep for Enrichment Kit (Catalog#: 20394426) for library construction according to the official protocol. Finally, we used an Illumina Nova 6000 for PE 150 sequencing.

Chromium 10x Genomics library and sequencing.

For 10X Chromium Single Cell Immune Profiling, target tissue was prepared into single-cell suspensions by using commercial enzymes from Miltenyi. Then, the cell concentration and viability were determined using a Countess II Automated Cell Counter and recorded. The cells were diluted to 700–1200 cells/µl. Cell loading, GEM generation and library construction of T cells and B cells were performed based on the Chromium Next GEM Single Cell 5' Kit official protocol. After library construction, we used an Agilent 2100 Bioanalyzer to perform a library check. Finally, we used an Illumina Nova 6000 for PE 150 sequencing. Additional sequencing was performed to obtain coverage of at least 50K reads per cell.

WES analysis.

Raw reads may contain adapter sequences and low-quality bases, which may affect the subsequent analysis. Therefore, it is necessary to review the raw sequencing reads. In this project, read quality was calculated for all samples by using FastQC (version 0.11.2) ¹³ software with the default parameters. We used the BWA-MEM (Burrows–Wheeler Aligner-maximal exact match)¹⁴ algorithm to compare the sequencing data to the reference genome (GRCh38), obtain the bam file, mark the PCR reads in the bam file, rearrange the regions that may have indel mutations, and recalibrate the quality of each base pair. Gene variants refer to changes in the sequence of the genome, including point mutations caused by single base changes and deletions, duplications, and insertions of multiple bases. In this analysis, the Sentieon (201611.02) ¹⁵ workflow was carried out to detect SNVs (single nucleotide variants) and indels (small insertions and deletions). The TMB result was the sum of nonsynonymous mutations and indel numbers in the exon region divided by the length of the total exon.

Single-cell RNA sequencing data alignment.

The 10X Genomics single-cell raw sequencing paired-end reads were processed using CellRanger software version 5.0.1 and the 10x human genome GRCh38 release as the reference¹⁶. The output-filtered gene expression matrices were analyzed by R (version 4.0.3) with the Seurat package (version 4.0.0). Low-quality cells with (i) < 500 UMIs (unique molecular identifiers), (ii) > 25% UMIs derived from the mitochondrial genome, or (iii) > 5000 or < 200 genes detected (genes expressed at a proportion > 0.1%) were removed. Downstream analyses, including normalization, scaling, clustering of cells, and identifying cluster marker genes, were performed by using the function in Seurat. The gene expression matrices were normalized by the NormalizeData function, then the FindVariableFeatures function was used to select 2000 features with high cell-to-cell variation. To reduce the dimensionality of the datasets, the RunPCA function was conducted with default parameters, including the selected 2000 highly variable genes on linear transformation-scaled data generated by the ScaleData function. Finally, we clustered cells using the suitable PCs to calculate the Uniform Manifold Approximation and Projection (UMAP) for dimension reduction with the FindNeighbors and FindClusters functions. Nonlinear dimensional reduction was performed with the RunUMAP function with default settings.

Cell type annotation and cluster marker identification.

Cells could be clustered together according to common features after nonlinear dimensional reduction. Then, all cells were projected into two-dimensional space by UMAP. We used the FindAllMarkers function in Seurat to find the marker genes with significant differences in expression for each cluster. All clusters were classified and annotated on the basis of these significantly different canonical markers. Clusters that expressed no canonical cell type markers were classified as “Others”.

Subclustering of major cell types.

The subset function in Seurat was used to extract the subcluster based on the expression of canonical markers of particular cell types. Next, these major cell types were integrated for further subclustering. After integration, the genes were scaled to unit variance. Scaling, PCA, and clustering were performed as described above.

Cell number statistics.

R was used to calculate the numbers of cells of different cell types. The ggplot2 library was used to generate relevant graphical outputs. The t test was used to test the significance between different groups.

Identification of malignant epithelial cells.

The R code inferCNV (https://github.com/broadinstitute/infercnv) was used to infer the CNVs in epithelial cells from the scRNA-seq data. The epithelial cells (EPCAM+) of the normal sample were considered putative nonmalignant cells, and their CNV estimates were used to define a baseline. The epithelial cells from tumor samples were compared to the normal reference baseline. The CNV R-scores in different regions are readily apparent and are overabundant or less abundant in the regions of the tumor genome compared to that of putative nonmalignant cells.

The calculated CNV signal was defined as the mean square of the CNV estimates across all genomic locations. The calculated CNV R-scores were defined as the Pearson correlation coefficient between each cell’s CNV pattern and the average CNV pattern of the top 5% of cells from the same tumor with respect to CNV signal ¹⁷. Epithelial cells (EPCAM+) with CNV R-scores ≥ 0.4 and CNV signal ≥ 0.03 were defined as malignant cells, while those with CNV R-scores ≤ 0.4 and CNV signal ≤ 0.03 were defined as nonmalignant cells. Cells that did not meet either of the above two conditions were defined as unresolved cells.

CellPhoneDB.

CellPhoneDB was used with default parameters to infer the interactions between the different kinds of cells on a linear transformation scaled data matrix. The method infers potential receptor–ligand interactions based on the expression of a receptor from one cell type and a ligand from another¹⁸. Only the receptors and ligands that were expressed in more than 30% of the cells in the specific cluster were considered. Network visualization was performed by using Cytoscape (version 3.8.2). The network layout was set to a force-directed layout.

Trajectory analysis. Trajectory analysis was performed using Monocle version 3.0.2.3 ¹⁹. We ordered the subclusters of cells, such as T cells or epithelial cells, onto a pseudotime trajectory based on the union of highly variable genes obtained from the cells of the subcluster, which were from normal, MIA and IAC samples. All parameters of Monocle3 were set to the defaults (variable genes with the following cutoff criteria: dispersion_empirical > dispersion_fit; and mean expression > 0.001). Dimensional reduction and cell ordering were performed using the DDRTree method and the orderCells function.

Data processing of single-cell VDJ libraries.

TCR and BCR reads were aligned to the GRCh38 reference genome, and consensus VDJ annotation was performed using Cellranger vdj (version 5.0.1) with default parameters. The R library Immunarch (version 0.6.7) was used to calculate the distribution of the V gene clonotype and calculate statistics for some related indicators, including overlap and correlation.

Statistical analyses.

The statistical tools, methods, and thresholds used for each analysis are explicitly described within the Results or detailed in the figure legends or Materials and Methods.

1. Global cellular landscape in multiple primary lung cancer.

To elucidate the cellular development in multiple primary lung cancer (MPLC), we analyzed a total of 167397 cells from 23 sampling sites in 6 patients with multiple foci of early-stage lung cancer from the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital (Fig. 1a-b, Supplementary Fig. 1a, d and Supplementary Table 2). Among these cells, 62714 (23.2%) were from AIS/MIA, 65844 (37.5%) were from IAC, and 38839 (39.3%) were from Normal Lung. To increase the accuracy of cell type designation, we jointly applied principal component analysis (PCA) after correction for read depth and mitochondrial read counts. Through marker-based annotation, 8 major cell types were identified in these unbatched and comparable datasets by graph-based UMAP (Fig. 1b, Supplementary Fig. 1a, 1d): epithelial cells, alveolar cells, fibroblasts, endothelial cells, myeloid cells, T lymphocytes cells, B lymphocytes cells and mast cells. Global cell type annotations for MPLC patients were categorized according to the molecular subtype of lung tumor tissues. The marker-based annotation correlated well with the cell types obtained via Cell Marker and SingleR v.1.0 ^20,21 (Fig. 1b and Supplementary Fig. 1c). T cells and epithelial cells show great differences over the course of development of MPLC. The relative abundance of B cells increased, while that of T cells decreased, from Normal Lung to AIS/MIA and IAC (Fig. 1d, Supplementary Fig. 1b, and Supplementary 1c). The most frequently mutated genes were EGFR (53%) and RBM10 (35%) (Fig. 1e and Supplementary Fig. 1g), and the most frequent variant type was missense mutation (Supplementary Fig. 1g). C > T and C > A transversions were enriched in these patients (Supplementary Fig. 1f and Supplementary Fig. 1g). The correlation of the variant allele frequency (VAF) for each individual focus was not significant (Fig. 1f). The significant mutations and the correlation of the VAF show that the foci in MPLC are independent and do not share an origin.

2. Signatures and metabolic disturbances of malignant cells in MPLC.

A total of 13176 normal epithelial cells and 14553 malignant epithelial cells were obtained from Normal Lung, AIS/MIA and IAC samples, and large-scale copy number variations (CNVs) were inferred with Normal Lung cells as references^22–24. Malignant epithelial cells were identified by inferring multipatient copy number variations (CNVs) with normal epithelial cells (mainly immune cells, including T cells and B cells) as references. The CNV patterns showed that malignant cells formed clusters in patients with different origins, indicating a high degree of intertumoral heterogeneity (Fig. 2a, Supplementary Fig. 2c, 2e). The clusters of malignant cells from different foci were independent and rarely overlapped, while the clusters of normal cells were mixed (Fig. 2a, Supplementary Fig. 2a, 2c). Alveolar cells are epithelial cells in the lung and play a pivotal role in maintaining the integrity and function of the alveoli. Both alveolar type I cells and alveolar type II cells were positive for EPCAM gene expression ^25,26. Other EPCAM + cells, including club, ciliated, basal, goblet, mucous, serous, ionocyte, and neuroendocrine cells, were categorized as epithelial cells. We identified them according to the expression of canonical gene markers (Fig. 2b). We clustered the normal epithelial cells as alveolar type I cells (AT1; AGER+), alveolar type II cells (AT2; SFTPA1+), ciliated airway epithelial cells (ciliated; TPPP3+), secretory club cells (club; SCGB1A1+) and other cells (no canonical markers to classify) based on the previously described markers (Fig. 2a, Supplementary Fig. 2a, and Supplementary Fig. 2d) ^17,27, independent of patient origin. The relative abundance of malignant epithelial cells, i.e., Ciliated, Club and AT1 cells, was increased in malignant cells compared with Normal Lung cells, and the number of AT2 cells was decreased from Normal Lung to AIS/MIA and IAC. We used lung cancer cells to recapitulate the multilineage differentiation processes of normal epithelia, and subclustering analysis of all epithelial cell datasets revealed divergent differentiation states in MPLC (Fig. 2d). Intriguingly, the projection of tissue cells and pathological feature classification onto the trajectory assigned the upper right branch and the lower right branch as malignant cells, mainly from IAC tissue. The trajectory indicated that there are two states in the evolution of malignant cells from AIS/MIA to IACs (Fig. 2c). MSI showed slight changes along the total pseudotime trajectory in MPLC, but TMB did not. In the pseudotime model, this trajectory began from nonmalignant AIS/MIA cells. Then, the trail branched in two directions: some malignant cells became state 1, and some malignant cells became state 2 in IAC. Branch 1 exclusively expressed immune response markers (e.g., CEACAM6, CXCL14, HLA_DPB1, CD74, S100A14, SPINK5, etc.), while branch 2 expressed cell death markers (e.g., PDCD4, AKR1C1, IER3, RGCC, APOH, etc.) (Supplementary Fig. 2g). Functional gene set enrichment analysis showed that the transition from State 1 was involved in the immune response, immune system process and response to transforming growth factor (Fig. 2e). In addition, the expression of genes implicated in cell death, response to stress and cell cycle increased in State 2 (Fig. 2f), which shows that some malignant cells will respond to stress from the immune response, some will undergo programmed cell death, and some will enter the cell cycle due to selective pressure during the development of tumors in MPLC.

3. Clonal dynamics and phenotypic transitions of tumor-infiltrating T cells in MPLC.

A total of 81,971 T cells were obtained from all samples and clustered into 10 subtypes. Three subtypes of CD4 + T cells were further clustered as Treg cells (TIGIT+, ENTPD1+, FOXP3+), Trm cells (CXCR6+) and CD4 + naïve T cells (CCR7+, IL7R+). Three subtypes of CD8 + T cells were further clustered as activated T cells/exhausted T cells (TNF+, IFNG+, GZMA+, GZMB+, KLRD1+), activated/effector T cells (TNF+, GZMA+, EOMES+), and effector T/Trm cells (ITGAE+, CXCR6+, GZMA+, GZMB+). One subtype of both CD4- and CD8-positive T cells was further clustered as naïve/memory T cells (IL7R+), and 2 subtypes with no CD4 or CD8 expression were further clustered as NK cells/exhausted T cells (HAVCR2+, IFNG+, GZMB+) and cytotoxic T cells/NK cells (XCL1+, XCL2+) (Fig. 3a, Supplementary Fig. 3a) based on canonical markers (Fig. 3b,3c, Supplementary Fig. 3b), as described previously ^28–32. The relative abundances of classified naïve T cells, Treg cells and Trm cells were increased, while those of activated/exhausted T cells were decreased, from Normal Lung to AIS/MIA and IAC (Fig. 3d). We next analyzed the lineage relationships between phenotypes and clonotypes in T cells. Globally, we found that cells in normal and IAC samples shared more common clonotypes in TCRs than AIS/MIA-group cells (Fig. 3e, Supplementary Fig. 4a). To investigate whether the development of tumors within clones could inform lineage transitions between T-cell phenotypes, we aggregated all clonotypes and measured the fraction shared in different clusters. There was minimal significant clonotype sharing between T cells in different cell groups. Trm and naïve T cells shared more significant clonotypes with other T cells, and the shared clonotypes increased with the development of infiltration in the tumor tissue, suggesting an increasing transition from Trm and naïve T cells to other activated T-cell states (Supplementary Fig. 4a). The T-cell phenotypes in AIS/MIA samples include more clonotype 1, for which TRB-TRA is “CASSGLAAKPGELF-CAVRRGQNFV”, than other normal and IAC samples (Supplementary Fig. 4b). This sequence may play an important role in antigen presentation during tumor invasion. In normal tissue, there is not enough antigen from tumors to stimulate T cells, and the expression of this antigen presentation factor is low. In IAC tissue, the tumor has experienced stress from T cells in AIS/MIA, and new subclone evolution may not cause the same antigen stimulation. Therefore, clonotype 1, in which TRB-TRA is “CASSGLAAKPGELF-CAVRRGQNFV”, was expressed significantly only in AIS/MIA samples. The TCR clonotype shared more overlap in focus grouped according to the source of the patient compared to those grouped by anatomical pathology info and still had more correlation (Fig. 3f).

4. Trajectory profiling reveals the branched progression of T cells in MPLC.

We next used Monocle3 to analyze the CD8 + and CD4 + T-cell trajectories to further investigate the potential transitions in each of these cell types. The pseudotime trajectory axis derived from Monocle3 indicates that CD8 + naïve and memory T cells transdifferentiate into effector T cells first. Then, some of these cells become exhausted T cells, and others transform into cytotoxic T cells (Fig. 4a). We compared the expression of genes in the next-stage cells with those in the previous stage (Supplementary Fig. 3c, 3d). Cells on the memory T to effector T branch expressed immune response markers (e.g., TNFRSF14, NFE2L2, TNFRSF1B, IL2RG, etc.) (Supplementary Fig. 3c,3d), while naïve T cells transitioning to effector T cells expressed lymphocyte activation and T-cell differentiation markers (e.g., FOXP1, TGFBR2, TNFRSF4, PRDM1, etc.) (Supplementary Fig. 3c, 3d). Cells on the effector T to exhausted T branch expressed responses to stress and cell death markers (e.g., GNLY, CCL4, and FCGR3A) (Supplementary Fig. 3c, 3d), while effector T cells expressed cytokine production markers (e.g., CXCR4, CRTAM, etc.) (Supplementary Fig. 3c, 3d). Gene Ontology enrichment analysis confirmed that the transition from memory T cells to effector T cells was mainly implicated in the immune response, immune system process and immune response regulation. The transformation from memory to effector in T cells is mainly a process of immune response and regulation (Fig. 4b (left-upper)). The function of the geneset in the transformation from naïve T to effector T cells was mainly related to lymphocyte activation, T-cell and leukocyte differentiation, T-cell activation and immune system development, which showed that the transformation from naïve to effector T cells was mainly a process of leukocyte and T-cell differentiation (Fig. 4b (upper right)). Gene Ontology enrichment analysis confirmed that signature genes were enriched for response to stress, positive regulation of cell death and response to stimulus pathway during effector T-cell transformation to exhausted T cells (Fig. 4b (lower right)). The functional gene set enrichment analysis from effector to cytotoxic T cells showed that the transition was mainly implicated in the regulation of cytokine production, which included the positive regulation of myeloid leukocyte cytokine production involved in the immune response, immune effector processes, macrophage cytokine production and the positive regulation of cytokine production involved in the immune response (Fig. 4b (lower left)). The pseudotime trajectory axis indicates that CD4 + naïve and Treg cells have four final differentiation states in themselves, which are labeled as Naïve_C1 to Naïve_C4 and Treg_C1 to Treg_C4 (Fig. 4c, Supplementary Fig. 5a). Pathway analyses showed that Naive_C2 cells were enriched in some effector pathways [Hallmark pi3k_akt_mtor signaling, chemokine signaling pathway, and natural killer cell-mediated cytotoxicity], while Naive_C4 cells did not express some pathways [T-cell receptor signaling pathway, regulation of lymphocyte differentiation, regulation of T-cell differentiation, T-cell differentiation, T-cell proliferation, and T-cell activation]. The function of the T cells in the Naïve_C2 cluster is closer to that of effector or exhausted T cells, while the function of the Naïve_C4 cluster is more closely related to that of other naïve clusters (supplemental Fig. 5b (Left)). The GSVA results showed that Treg_C3 cells were enriched in the regulation of myeloid cell differentiation [myeloid leukocyte differentiation, myeloid cell differentiation, and regulation of myeloid cell differentiation], while Treg_C2 cells mainly participated in the regulation of cell differentiation, activation of the T-cell receptor signaling pathway [regulation of cell differentiation, T-cell receptor signaling pathway, T-cell receptor signaling pathway and T-cell activation]. Treg_C2 cells also had high expression levels of genes related to lymphocyte differentiation, including lymphocyte differentiation, regulation of T-cell differentiation, and regulation of lymphocyte differentiation. Treg_C1 cells mainly participates in the regulation of stem cell differentiation [regulation of stem cell differentiation, stem cell differentiation, and hematopoietic progenitor cell differentiation] (Supplementary Fig. 5b (Right)).

5. B cells have two functional distribution subtypes in MPLC.

A total of 4279 B cells were analyzed and classified into five subtypes (Fig. 5a, Supplementary Fig. 6a, b, supplementary Table 3): Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) (49.7% in IAC, 10.3% in AIS/MIA, 6.4% in Normal); Bcells_C2 (NEIL1+, RGS13+) (0.9% in IAC, 0.4 in AIS/MIA, 0.0% in Normal); Bcells_C3 (CYTOR+, NUCB2+, ANKRD28+, TRIB1+) (7.8% in IAC, 2.9% in AIS/MIA, 1.9% in Normal); Bcells_C4 (6.2% in IAC, 3.8% in AIS/MIA, 0.9% in Normal); and Bcells_C5 (IL7R+, IL32+, CST7+, FYB1+) (5.9% in IAC, 1.5% in AIS/MIA, 1.4% in Normal) (Supplementary Fig. 6e, Supplementary Fig. 7a). The relative abundance of the B-cell subtype showed no uniform pattern in different patient foci and or when grouped by anatomical pathology (Fig. 5c, Supplementary Fig. 6c, d). The pathway analysis of B cells in MPLC further confirmed that there are two kinds of B cells with different functions. The regulation and activation of the innate immune response was highly activated in the Bcell-C4 and Bcell-C5 clusters (Fig. 5d, red pathway). However, the Bcell-C1, Bcell-C2, and Bcell-C3 clusters were enriched in pathways related to antigen binding, processing and presentation, including Reactome MHC class II Antigen Presentation (Fig. 5d, green pathway). Some B cells are involved in immune activation and regulation of the immune response, while others are involved in antigen presentation. We next used Monocle3 to analyze the trajectory performance of B cells to further investigate the potential transition in each of their cell types. The pseudotime trajectory axis derived from Monocle3 indicates that Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) have two differentiation states; one is transdifferentiating into Bcell_C3 (NEIL1+, RGS13+) and Bcell_C5 (IL7R+, IL32+, CST7+, FYB1+), and the other is transitioning to Bcell_C4 (Fig. 5b, Supplementary Fig. 7b, 7c, 7d). Combined with the pathway enrichment results, these findings show that B cells mainly differentiated into two different types of cells with different functions in MPLC; one is involved in immune activation and regulation of the immune response, while the other is involved in antigen presentation. Different from the TCR clonotypes, which shared more overlap and correlation in foci grouped according to the source within the patient, the BCR clonotypes shared little overlap or correlation in MPLC (Supplementary Fig. 6f).

6. Myeloid cells have two functional distribution subtypes in MPLC.

Subclustering of 46836 myeloid cells revealed 10 subsets (Fig. 5e, Supplementary Fig. 9a.): M1 macrophages (MSR1+, MARCO+, etc.), M2 macrophages (MERTK+, MMP9+, etc.), MDSCs (CEACAM8+, ITGAM+, etc.), monocyte cells (TNFSF10+, FCER2+, CLEC10A+, etc.), dendritic cells (FCER1A+, CD1C+, CLEC4A+, etc.), myeloid dendritic cells (CCR7+, IL15+, CD1C+, etc.), plasma cells (IRF4+, CD36+, CXCR4+, etc.), secretory cells (CXCL10+) (Fig. 5g, Supplementary Fig. 9d). The percentage of M2 macrophages increased from normal to AIS/MIA to IAC, while the other subtypes were not significantly different in MIA/AIS and IAC (Fig. 5f, Supplementary Fig. 9b). Macrophages play an important role in maintaining organismal integrity by either directly participating in pathogen elimination or repairing tissue under sterile inflammatory conditions^33, and Microphage_2 cells support tissue integrity by providing growth factors and healing capacity³⁴. From a functional perspective, myeloid cells are divided into two main types: one exhibits high expression of IL1A, FCGR3A and IFN-γ (M1 macrophages), which can positively regulate the innate immune response. The other negatively regulates the adaptive immune response and regulation of stem cell differentiation (M2 macrophages). The increase in M2 macrophages will lead to malignant cells adapting to the immune response and escaping the response from the immune system (Supplementary Fig. 9c). The pseudotime trajectory axis derived from Monocle3 indicates that myeloid cells have four main branch trajectories (Fig. 5h). M2 macrophages were mainly concentrated in the left branch trajectory, M1 macrophages were mainly in the opposite branch trajectory, monocyte cells were in the middle to left stage, secretory cells were mainly concentrated in the lower left branch trajectory, and other Myeloid2_cells were mainly in the lower right branch trajectory (Fig. 5h). Next, we evaluated the expression of DEGs in the M1 macrophage (MSR1+; MARCO+) and M2 macrophage (MERTK+; MMP9+) populations in the whole trajectory. The expression of macrophage colony-stimulating factor genes, such as CSF1, MCEMP1, APOE, CD52 and FABP4, increased gradually during the process of transdifferentiating into M1 macrophages, and the expression of genes related to cytomembrane transport, namely, GPR183, CD48, CD93, and CCL5, increased during the process of transdifferentiating into M2 macrophages (Supplementary Fig. 9f).

7. Positive regulation of endothelial cell development is extensive in IAC.

We found four distinct subtypes by reclustering 2081 endothelial cells (Fig. 6a, Supplementary Fig. 8a): Endothelial_C1 (VIPR1+, IL1RL1+, TMEM100+, TNFSF10+, etc.), Endothelial_C2 (EMP1+, CD93+, PLVAP+, etc.), Endothelial_C3 (TNFAIP3+, CXCR4+, CREM+, etc.), and Endothelial_C4 (TFF3+, TIMP2+, STMN1+, etc.) (Fig. 6c, Supplementary Fig. 8d). The percentage of Endothelial_C1 and Endothelial_C2 in MIA/AIS and IAC were more different than those in normal samples, and the distribution of endothelial subtypes was not significantly different in MIA/AIS and IAC (Fig. 6b, Supplemental Fig. 8c). Then, we focused on the transcriptomic features and functions of Endothelial_C2 cells. The genes EMP1, TSPAN7, VWF, PLVAP and CD93 were expressed at high levels in these cells, and pathway analyses showed that Endothelial_C2 cells were enriched in positive regulation of endothelial cell development pathways [positive regulation of blood vessel diameter, regulation of endothelial cell differentiation pathway, and abnormal morphology of the great vessels] (supplemental Fig. 8d, Fig. 6d). Thus, Endothelial_C2 cells are related to the regulation of endothelial development, and their prevalence increases with the development of tumor invasion.

8. Fibroblast subtype distributions are similar in IAC, AIS/MIA and normal samples.

Fibroblast cells maintain a very stable state in MPLC (Fig. 6c, Supplementary Fig. 8c). A total of 2523 fibroblast cells were analyzed, and five subtypes were identified (Fig. 6c). In MPLC, there was no significant difference in the number of different fibroblast clusters, and there were also no significant changes in the pathway and function of genes expressed in different fibroblast clusters (Fig. 6f, 6g, Supplementary Fig. 8f). Fibroblast_C1 cells expressed high levels of genes associated with the Reactome elastic fiber formation annotation, such as COL3A1, COL1A1, COL6A3, COL10A1, and COL8A1, whereas Fibroblast_C5 cells highly expressed genes related to positive regulation of fibroblast proliferation, like RASD1, SFTPB and FOS (Fig. 6g, Supplementary Fig. 8g).

9. Characterization of cell–cell interactions involved in MPLC.

There are many specific signaling network changes in MPLC, which depend on complex interactions between structural and inflammatory cells³⁵. One cell expresses a receptor or ligand, and this “ligand–receptor” pair will be defined as an incoming or outgoing interaction for this cell. We used CellPhoneDB to identify potential cell–cell interactions and define their changes in primary lung cancer. Compared with the MIA samples, the cells in tumors exhibited significantly increased interactions (Fig. 7a). Alveolar and myeloid cells had the most outgoing interactions of all cell types in both MIA and IAC samples. In comparison with MIA, the interaction between T cells and alveolar cells was significantly increased (Fig. 7b). Myeloid, cytotoxic T, natural killer and B cells are involved in the action of recognizing and presenting the antigens of epithelial cells among the interactions from immune cells to epithelial cells. Analysis of the biological functions revealed that interactions related to antigen recognition and presentation (CD74_COPA, HLA-C_FAM3C, and HLA-DRB1_OGN) were less abundant in IAC than in MIA and normal samples, which means that some escape of immune recognition occurred with deeper invasion in MPLC (Fig. 7c upper). The interactions related to independent cellular adhesion (NECTIN1_CADM3, FGF1_FGFR4, and CD44_FGFR2) were still less abundant in epithelial cells than in fibroblasts, which means that epithelial cells promoted structural changes in fibroblasts and reduced cell adhesion. The interactions regulating the growth and differentiation of numerous types of cells (ACVR1_BMPR2_BMP6 and ACVR1_BMPR2_BMP5) were more abundant (Fig. 7c, middle). On the one hand, epithelial cells reduce fibroblast adhesion; on the other hand, they strengthen the differentiation of other mesenchymal cells with deeper invasion in MPLC. The interactions related to endothelial and hematopoietic stem cell differentiation and development (SELL-CD34, SELP-CD34, NPR1-NPPC, and NPR2-NPPC) were more abundant in IAC than in normal and MIA samples (Fig. 7c, down).

Our knowledge about the dynamics of gene expression and metabolic pathways during cellular changes and evolution in MPLC remains incomplete. In this study, we comprehensively provide important new insights into these changes during the progression from normal AIS to IAC in MPLC.

Recent studies have performed similar development trajectory analyses for CD8 + T cells, with naïve T cells as the root and exhausted T cells and cytotoxic T cells as the two end states ^17,36,37. However, naïve T cells are not the only root for CD8 + T-cell development and evolution in MPLC. We found that memory T cells are also a main source of these cells and play an important role in CD8 + T-cell differentiation. The immune response and regulation of the immune response were significantly enriched during this differentiation process. The transition from CD8 + memory T cells to effector T cells were mainly a process of positive regulation of immune system activities; however, the progression of CD8 + naïve to effector T cells is mainly a process of leukocyte activation and T-cell differentiation. The progress of effector T-cell transformation to exhausted T cells involves responses to stimuli and positive regulation of cell death. The other pathway, of effector T-cell transformation to cytotoxic T cells, mainly involves the regulation of cytokine production, including the positive regulation of myeloid leukocyte cytokine production involved in the immune response and the positive regulation of cytokine production involved in the immune response. CD4 + T cells have two distinct functional statuses: Treg and naïve. We also found that T cells shared more consistent TCR clonotypes than AIS/MIA-group cells in normal and IAC samples. Unfortunately, further in-depth study of this interesting phenomenon is beyond the scope of this study.

In contrast to previous studies on the metabolism of malignant cells ^17,38−40, we used trajectory analysis to study the transition of epithelial cells from normal to malignant. We found that although different nodules in the patients were heterogeneous, they all exhibit two states during the evolution of malignant cells from AIS/MIA to IAC. Specifically, some cells respond to stress from immunity, while others undergo programmed cell death or enter the cell cycle due to selective pressure from the development of tumors in MPLC. The proportion of AT2 cells gradually decreases with the development of LUAD from normal to AIS/MIA to IAC, which also provides evidence that AT2 is the most likely origin of malignant cells in MPLC ³².

The functional distribution subtypes of B cells can be distinguished into two different stages. One is the regulation and activation of innate immune response antigen binding, and the other is processing and presentation, including MHC class II antigen presentation. Myeloid cells are divided into two major types, one with high expression of positive regulators of the innate immune response and the other performing negative regulation of the adaptive immune response and regulation of stem cell differentiation presented by M2 macrophages. With the progression of LUAD from normal to AIS/MIA to IAC, the increasing proportion of M2 macrophages causes malignant cells adapt to the immune response and escape the immune response ^41–43. We also found that some endothelial cell clusters also participate in the positive regulation of blood vessel development in early glandular neoplasia of the lung. Moreover, the proportion of these endothelial cells gradually increases with the progression of LUAD from normal to AIS/MIA to IAC ^44–46.

Comprehensive analysis of the cell–cell interactions in MPLC identified that interactions related to antigen recognition and presentation were less abundant in IAC than in AIS/MIA and normal samples, and the interactions related to independent cellular adhesion were still less abundant in epithelial cells than in fibroblasts. The interactions regulating the growth and differentiation of numerous types of cells were more abundant, and the interactions related to endothelial and hematopoietic stem cell differentiation and development were more abundant in the IAC than in the normal and AIS/MIA samples.

Despite the striking findings identified in this study, there are still several limitations, including sample collection, technical challenges and software limitations. We were unable to evaluate the continuous development of different nodules in the same patient, and the cells used for single-cell sequencing were from just a part of the tumor tissue, so there might be differences in biological complexity between our samples and the tissues used for pathological diagnosis and real cells in developing tumors. The 10x Genomics single-cell technology cannot cover all the transcripts in cells. There are also still limitations in batch effect removal, clustering and trajectory inference analysis, as these analyses are still maturing, and the complexity of the underlying topology could be underestimated ⁴⁷.

Data availability

The raw data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

We appreciate the support and participation of the physicians and patients in this study.

Competing interests

Hao Xia is employee of Mingma Technologies Co., Ltd. All remaining authors have declared no conflicts of interest.

Author contributions

W.G. and S.G. conceived the study. J.H., J.Y., F.T., Q.X., and S.G. provided project management and supervision. W.G. and L.G. provided or facilitated the accrual of patient samples, pathology, and/or clinical data. W.G. and H.X. performed bioinformatics and genomic analyses. W.G. and H.X. performed statistical analyses. W.G. and Z.B. wrote the original draft, with input from all authors.

Funding

This study was supported by the National Key R&D Program of China (2021YFC2500900), the National Natural Science Foundation of China (82002451), the CAMS Initiative for Innovative Medicine (2021-1-I2M-015), the Special Research Fund for Central Universities, Peking Union Medical College (3332020024), the Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences (2018PT32033), and the Beijing Hope Run Special Fund of Cancer Foundation of China (LC2019B15).

Chen, W. et al. Cancer statistics in China, 2015. CA: a cancer journal for clinicians 66, 115–132 (2016).
Siegel, R. L., Miller, K. D., Fuchs, H. E. & Jemal, A. Cancer statistics, 2021. CA: a cancer journal for clinicians 71, 7–33 (2021).
Sung, H. et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71, 209–249 (2021).
Shintani, Y. et al. Clinical features and outcomes of patients with stage I multiple primary lung cancers. Cancer Sci 112, 1924–1935, doi:10.1111/cas.14748 (2021).
Liu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nat Commun 7, 13200, doi:10.1038/ncomms13200 (2016).
Port, J. L. et al. Surgical resection for multifocal (T4) non-small cell lung cancer: is the T4 designation valid? The Annals of thoracic surgery 83, 397–400 (2007).
Asamura, H. Multiple primary cancers or multiple metastases, that is the question. J Thorac Oncol 5, 930–931, doi:10.1097/JTO.0b013e3181e15d8c (2010).
Meacham, C. E. & Morrison, S. J. Tumour heterogeneity and cancer cell plasticity. Nature 501, 328–337 (2013).
Ma, P. et al. Simultaneous evolutionary expansion and constraint of genomic heterogeneity in multifocal lung cancer. Nature communications 8, 1–11 (2017).
Liu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nature communications 7, 1–8 (2016).
Murphy, S. J. et al. Identification of independent primary tumors and intrapulmonary metastases using DNA rearrangements in non–small-cell lung cancer. Journal of Clinical Oncology 32, 4050 (2014).
Giraldo, N. A. et al. The clinical role of the TME in solid cancer. British journal of cancer 120, 45–53 (2019).
S. A., S. B. a. FastQC: a quality control tool for high throughput sequence data. (2010).
Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589–595, doi:10.1093/bioinformatics/btp698 (2010).
Weber JA, A. R., Gallagher BD, Edwards JS. Sentieon DNA pipeline for variant detection - Software-only solution, over 20× faster than GATK 3.3 with identical results. PeerJ PrePrints 4:e1672v2, doi:https://doi.org/10.7287/peerj.preprints.1672v2 (2016).
https://github.com/10XGenomics/cellranger. Cellranger.
Xing, X. et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv 7, doi:10.1126/sciadv.abd9738 (2021).
Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc 15, 1484–1506, doi:10.1038/s41596-020-0292-x (2020).
Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–502, doi:10.1038/s41586-019-0969-x (2019).
Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20, 163–172, doi:10.1038/s41590-018-0276-y (2019).
Zhang, X. et al. CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res 47, D721-D728, doi:10.1093/nar/gky900 (2019).
Patel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396–1401, doi:10.1126/science.1254257 (2014).
Puram, S. V. et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 171, 1611–1624 e1624, doi:10.1016/j.cell.2017.10.044 (2017).
Timothy Tickle, I. T., Christophe Georgescu,Maxwell Brown,Brian Haas. Visualizing Large-scale Copy Number Variation in Single-Cell RNA-Seq Expression Data. (2021).
Khan, P. et al. Culture of human alveolar epithelial type II cells by sprouting. Respir Res 19, 204, doi:10.1186/s12931-018-0906-9 (2018).
Travaglini, K. J. et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619–625, doi:10.1038/s41586-020-2922-4 (2020).
Deprez, M. et al. A Single-Cell Atlas of the Human Healthy Airways. Am J Respir Crit Care Med 202, 1636–1645, doi:10.1164/rccm.201911-2199OC (2020).
Morgan, A. J. et al. Expression of CXCR6 and its ligand CXCL16 in the lung in health and disease. Clin Exp Allergy 35, 1572–1580, doi:10.1111/j.1365-2222.2005.02383.x (2005).
Wein, A. N. et al. CXCR6 regulates localization of tissue-resident memory CD8 T cells to the airways. J Exp Med 216, 2748–2762, doi:10.1084/jem.20181308 (2019).
Ghorani, E. et al. The T cell differentiation landscape is shaped by tumour mutations in lung cancer. Nat Cancer 1, 546–561, doi:10.1038/s43018-020-0066-y (2020).
Delacher, M. et al. Single-cell chromatin accessibility landscape identifies tissue repair program in human regulatory T cells. Immunity 54, 702–720 e717, doi:10.1016/j.immuni.2021.03.007 (2021).
Wang, Z. et al. Deciphering cell lineage specification of human lung adenocarcinoma with single-cell RNA sequencing. Nat Commun 12, 6500, doi:10.1038/s41467-021-26770-2 (2021).
Chavez-Galan, L., Olleros, M. L., Vesin, D. & Garcia, I. Much More than M1 and M2 Macrophages, There are also CD169(+) and TCR(+) Macrophages. Front Immunol 6, 263, doi:10.3389/fimmu.2015.00263 (2015).
Ley, K. M1 Means Kill; M2 Means Heal. J Immunol 199, 2191–2193, doi:10.4049/jimmunol.1701135 (2017).
Holgate, S. T. et al. Asthma. Nat Rev Dis Primers 1, 15025, doi:10.1038/nrdp.2015.25 (2015).
Guo, X. et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat Med 24, 978–985, doi:10.1038/s41591-018-0045-3 (2018).
Qian, J. et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell Res 30, 745–762, doi:10.1038/s41422-020-0355-0 (2020).
Goveia, J. et al. An Integrated Gene Expression Landscape Profiling Approach to Identify Lung Tumor Endothelial Cell Heterogeneity and Angiogenic Candidates. Cancer Cell 37, 421, doi:10.1016/j.ccell.2020.03.002 (2020).
Ireland, A. S. et al. MYC Drives Temporal Evolution of Small Cell Lung Cancer Subtypes by Reprogramming Neuroendocrine Fate. Cancer Cell 38, 60–78 e12, doi:10.1016/j.ccell.2020.05.001 (2020).
He, Y. et al. Mechanisms of Progression and Heterogeneity in Multiple Nodules of Lung Adenocarcinoma. Small Methods 5, e2100082, doi:10.1002/smtd.202100082 (2021).
Heusinkveld, M. & van der Burg, S. H. Identification and manipulation of tumor associated macrophages in human cancers. J Transl Med 9, 216, doi:10.1186/1479-5876-9-216 (2011).
Conway, E. M. et al. Macrophages, Inflammation, and Lung Cancer. Am J Respir Crit Care Med 193, 116–130, doi:10.1164/rccm.201508-1545CI (2016).
Liu, M. et al. Transcription factor c-Maf is a checkpoint that programs macrophages in lung cancer. J Clin Invest 130, 2081–2096, doi:10.1172/JCI131335 (2020).
Lamalice, L., Le Boeuf, F. & Huot, J. Endothelial cell migration during angiogenesis. Circ Res 100, 782–794, doi:10.1161/01.RES.0000259593.07661.1e (2007).
Betsholtz, C. Cell-cell signaling in blood vessel development and function. EMBO Mol Med 10, doi:10.15252/emmm.201708610 (2018).
Sweeney, M. & Foldes, G. It Takes Two: Endothelial-Perivascular Cell Cross-Talk in Vascular Development and Disease. Front Cardiovasc Med 5, 154, doi:10.3389/fcvm.2018.00154 (2018).
Lahnemann, D. et al. Eleven grand challenges in single-cell data science. Genome Biol 21, 31, doi:10.1186/s13059-020-1926-6 (2020).

Competing interest reported. Hao Xia is employee of Mingma Technologies Co., Ltd. All remaining authors have declared no conflicts of interest.

Download PDF

Version 1

posted

You are reading this latest preprint version

Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing

Status:

Version 1

Abstract

Figures

Introduction

Methods

Patients and sample collection

Results

Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1