Ductal carcinoma in situ (DCIS) is a precursor to invasive breast cancer. The frequency of DCIS is increasing because of routine mammography; however, the biological features and intratumoral heterogeneity of DCIS remain obscure. To address this deficiency, we performed single-cell transcriptomic profiling of DCIS and invasive ductal carcinoma (IDC). DCIS was found to be composed of several transcriptionally distinct subpopulations of cancer cells with specific functions. Several transcripts, including long noncoding RNAs, were highly expressed in IDC compared with DCIS and might be related to the invasive phenotype. Closeness centrality analysis revealed extensive heterogeneity in DCIS, and the prediction model for cell-to-cell interactions implied that the interaction network among luminal cells and immune cells in DCIS was comparable with that in IDC. In addition, transcriptomic profiling of HER2+ luminal DCIS indicated HER2 genomic amplification at the DCIS stage. These data provide novel insight into the intratumoral heterogeneity and molecular features of DCIS, which exhibit properties similar to IDC.

Significance:

Investigation of the molecular features of ductal carcinoma in situ at single cell resolution provides new insights into breast cancer biology and identifies candidate therapeutic targets and diagnostic biomarkers.

Ductal carcinoma in situ (DCIS) is a nonobligate precursor of invasive breast cancer (1). The diagnosis of DCIS is increasing, accounting for more than 20% of all breast cancers, because of the widespread use of screening mammography, and approximately 50,000 cases of DCIS are diagnosed each year in the U.S. (2). However, overtreatment is a problem, because only 10%–50% of untreated DCIS lesions are estimated to progress to invasive ductal carcinoma (IDC) during the lifetime of the patient (3–10). Prognostic markers are required to identify patients in need of treatment, but the molecular features and the mechanism of progression from DCIS to IDC remain unclear.

To provide a molecular characterization of the transition from DCIS to IDC, a number of studies have been conducted. In an observational study, tumor size, grade, subtype and the presence of comedo necrosis were reported as risk factors for DCIS recurrence (11) but are controversial. Genomic analysis revealed the similarity of gene mutations and copy-number aberrations between DCIS and IDC. They shared TP53, PIK3CA, and MYC mutations and gain of chromosomes 1q, 8q, 11q13, 17q12, and 20q13 (12, 13). Expression analysis was conducted by bulk RNA sequencing to compare overexpressed genes in IDC with those in DCIS (14). However, key genomic and genetic factors have not been identified. On the other hand, researchers have focused on the differences in the microenvironments between DCIS and IDC. They found differences in the populations of stromal cells and immune cells between DCIS and IDC, which may be related to invasiveness in breast cancer (15, 16). The interaction of the local tumor microenvironment with epithelial carcinoma cells in DCIS is considered to be important in invasive progression.

Single-cell RNA sequencing (scRNA-seq) is a powerful approach to investigate intratumoral heterogeneity and analyze gene expression profiles of rare populations. These approaches have revealed diversity in normal breast epithelial cells (17, 18), intratumoral heterogeneity in IDC, and diversity in the tumor ecosystem (19–24). In this study, we applied scRNA-seq to examine DCIS heterogeneity and the molecular features. We obtained transcriptome profiles of 1 normal breast tissue sample, 7 DCIS samples and 6 IDC samples at single-cell resolution. DCIS consisted of several epithelial clusters as well as immune cells, endothelial cells, and stromal cells. Differential gene expression analysis between DCIS and IDC identified candidate genes that might be related to the invasive phenotype of IDC. In addition, we reported the intercellular interaction of DCIS epithelial cells with immune cells. The findings from this study suggest possible novel biomarkers for diagnosing DCIS as well as provide insight into the basic cancer biology of DCIS.

Patients and tumor specimens

This study was approved by the National Cancer Center Hospital, and all patients provided written informed consent for the collection of specimens and detailed analyses of the derived genetic materials (2019–018). We conducted this study according to the Declaration of Helsinki. Six of the 14 patients were diagnosed with IDC and underwent breast-conserving surgery or total mastectomy without prior treatment. Two patients underwent lymph node dissection because of lymph node metastasis. Seven of the 14 patients were diagnosed with DCIS and underwent surgery, but a 6-mm invasion area was detected pathologically in one patient. We obtained 1 normal breast tissue sample from a patient who was diagnosed with DCIS and underwent total mastectomy. The normal sample was used to exclude admixture of possible normal cells in DCIS and IDC cases. We performed scRNA-seq on all 14 samples. A fresh-frozen tissue from DCIS sample was subjected to spatial transcriptome analysis. A trained breast pathologist made a histological diagnosis of the area sampled for digestion.

Breast cancer tissue resection

The surgical specimens of DCIS and IDC cases were kept in tissue collection buffer (RPMI-1640 medium with 5% FBS and 1% antibiotics and antimycotics) on ice. The breast cancer tissue fragments were mechanically minced and dissociated with enzymes according to the protocol of a Tumor Dissociation Kit human (Miltenyi Biotec) within 1 hour after the sample collection. Single cells were isolated with a Miltenyi Gentle MACS dissociator and were then incubated with enzymes for 30 minutes at 4°C. The cells were then washed with PBS supplemented with 0.04% BSA (Sigma-Aldrich), passed through 70- and 35-μm filters, centrifuged at 300 × g for 10 minutes, and resuspended in PBS supplemented with 0.04% BSA. The cells were not frozen, and cDNA generation using the Chromium Next GEM Single Cell 3′ Reagent Kit was performed within the same day. Cell viability was checked by trypan blue staining, and the viability was approximately 80% in our dissociation protocol. When viability was very low; one sample showed 30% of cell viability because of a lot of necrotic tissue, we did not use this sample for single-cell RNA-seq.

Chromium 10× Genomics library preparation and sequencing

Cells were counted and resuspended in PBS supplemented with 0.04% BSA for loading for single-cell library construction on the 10× Genomics platform (10× Genomics). All experimental samples were analyzed with a Chromium Next GEM Single-Cell 3′ GEM, Library and Gel Bead Kit v3.1 and Chromium Next GEM Chip G Single-Cell Kit according to the manufacturer's instructions in the Chromium Next GEM Single-Cell 3′ Reagent Kit V3.1 User Guide. Briefly, approximately 8,000 cells were loaded into each channel and were then partitioned into gel beads in emulsion in a Chromium Controller for cell lysis and barcoded reverse transcription of RNA, followed by amplification, fragmentation, and 5′ adaptor and sample indexing. Sequencing was performed on the HiSeq platform (Illumina) with an average sequencing depth of 89,940 reads/cell.

scRNA-seq data alignment

Raw sequencing data were processed with the CellRanger pipeline (version 3.0.2, 10× Genomics) and mapped to the human genome (GRCh38) to generate gene count matrices by cell barcodes.

Data quality control and normalization

scRNA-seq data analyses, including normalization, scaling, clustering of cells, and identification of cluster marker genes, were performed using the R software package Seurat version 4.0.3. We extracted single cells with nFeature_RNA >100 and percent.mt <20 to remove dead cells. We applied DoubletFinder (25) for removing doublet cells in the analysis. Log normalization was performed with the function “NormalizeData” in Seurat. The expression data were normalized by log transformation and scaled by the number of genes and percentage of mitochondrial reads. For data integration, “Harmony” was applied to integrate the scRNA-seq data by sample and reduce batch effects (26, 27).

Data clustering and dimensionality reduction

We performed principal component analysis for dimensionality reduction in R with Seurat. Clustering of single cells was performed with the functions “FindNeighbors” and “FindClusters” in Seurat. The dimensionality-reduced cell clustering was shown as a UMAP plot with the function “runUMAP.” Heatmaps were generated by Morpheus software from the Broad Institute (https://software.broadinstitute.org/morpheus/). A correlation plot was generated with the “corrplot” package in R.

Pathway enrichment analysis

We performed enrichment analysis with the signature gene list from each epithelial cell cluster using the ‘ClusterProfiler’ packages in R (28). For enrichment analysis, gene symbols were converted to ENTREZ IDs using the “org.Hs.eg.db” package (Carison M, R package version 3.10.0., 2019). GO enrichment analysis using the “enrichGO” function was performed by the BH method.

Closeness centrality analysis

Correlation network analysis was performed with the “igraph” package (29). We calculated the correlation coefficients for each cell population between IDC and DCIS. Putative normal cells (cells in Cluster 5) from DCIS and IDC were excluded in the analysis. In each cluster, normalized closeness centrality was calculated in R. Where r is the absolute value of Pearson correlation coefficient and n is the number of cells in the cluster.

formula

The detailed procedures were previously described (29, 30).

Cell–cell interaction analysis

Cell–cell interactions (CCI) between luminal epithelial cells from DCIS and IDC lesions and immune cell subpopulations were performed using the algorithms of CCI score calculation. Using previously published tool (31, 32) and we originally conducted the code for CCI with the interaction database (33). Putative normal cells (cells in Cluster 5) from DCIS and IDC were excluded in the analysis. Data visualization was performed with the “circlize” package in R (34).

Spatial transcriptomics

A tissue sample was embedded in optimal cutting temperature compound (Sakura Finetek Japan) and stored at −80°C. Tissue block was cut into 10-μm sections, which were analyzed with a Visium Spatial Gene Expression Kit (10x Genomics) according to the manufacturer's instructions. First, breast tissue permeabilization conditions were optimized using the Visium Spatial Tissue Optimization Kit (10x Genomics). The optimal permeabilization time was 24 minutes. Sliced tissue pieces were stained with hematoxylin and eosin, imaged using a BZ-X710 microscope (Keyence) under a ×10 objective, and processed for spatial transcriptomics. The resulting complementary DNA library was checked for quality control and was then sequenced using a DNBSEQ system (MGI Tech). The cycling conditions were set to 28 and 90 for Read 1 and Read 2, respectively.

Normalization and feature selection for spatial transcriptomics

Raw FASTQ files and histological images were processed using Space Ranger software v1.2.2. To assign individual spots to tumor cells or cells that composed the microenvironment, we compared the clusters morphologically annotated by a pathologist. Analyses of spatially resolved RNA-seq data, including normalization, dimensional reduction, clustering, detection of spatially variable features, and interactive visualization (scaling, clustering of cells, and identification of cluster marker genes) were performed using the R software package Seurat version 4.0.3. Sctransform (35) normalization was performed with the function “SCTransform” in Seurat (https://satijalab.org/seurat/articles/spatial_vignette.html).

Ingenuity Pathway Analysis (IPA; https://www.digital-biology.co.jp/allianced/products/ipa/) was used to analyze spatial transcriptomic data for predicting enriched canonical pathways and upstream regulators.

Immunohistochemical staining

Each breast cancer tissue sample was fixed with 4% paraformaldehyde overnight and embedded in paraffin. Breast cancer sections were dewaxed with xylene and rehydrated with ethanol (100%–70%). Antigen retrieval, blocking, incubation with anti-HER2 primary antibody and visualization were performed using HercepTest II (Agilent) according to the manufacturer's instructions.

Cell culture

The MCF7 and T47D breast cancer cell lines were purchased from the ATCC and tested to verify the absence of Mycoplasma contamination when cells were prepared. All the cell lines were authenticated by genetic profiling using polymorphic short tandem repeat loci (Promega). The cell lines were maintained in RPMI-1640 medium (Nacalai Tesque Inc.) supplemented with 10% heat-inactivated FBS (Thermo Fisher Scientific) and 1 × antibiotic and antimycotic solution (Thermo Fisher Scientific) and grown at 37 °C in a 5% CO2 incubator. The medium was changed every 3 days, and the cells were passaged when they reached 70% to 80% confluence. In this study, the passage number of cell lines is fewer 15 after purchase.

siRNA transfection

One day after seeding MCF7 and T47D cells (5×103 for 96-well plates and 4×104 for 12-well plates), the medium was replaced by medium containing siRNAs (RPL17, Cat# J-013633–11–0002, Cat# J-013633–09–0002; CPB1, Cat# J-005821–07–0002, Cat# J-005821–05–0002; siControl2, Cat# D-001206–14–05; Horizon) in DharmaFECT reagent-1 (Cat# T-2001–02; GE Healthcare Dharmacon, Inc.) according to the manufacturer's instructions. The final concentration of the siRNAs was 5 nmol/L. After 3 days of transfection, cell numbers were estimated by a CellTiter-Glo2.0 cell viability assay (Promega), and luminescence was measured using a microplate reader (Molecular Devices). For qRT–PCR, total RNA was collected on day 3 after transfection.

qRT-PCR

Total RNA was extracted using the RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol. RNA was subsequently reverse transcribed to cDNA using a high-capacity cDNA reverse transcription kit (Thermo Fisher Scientific). qPCR was then performed using Power SYBR Green PCR Master Mix (Thermo Fisher Scientific). The reaction conditions were 95°C for 1 minute, 95°C for 15 seconds, 56°C for 15 seconds, and 72°C for 30 seconds. The relative expression levels of mRNAs were calculated by the ΔΔCt method. β-Actin was used as the internal control. The sequences of the primers are shown in Supplementary Table S1.

Data availability

Transcriptomic data were deposited into the NCBI GEO database and NBDC JGA database under accession number GSE195861, GSE196208, and JGAS000309.

Statistical analysis

The Welch's t test was used for comparisons of two datasets. The Mann–Whitney U test was used for calculation of statistics in closeness centrality. Correlation coefficients were calculated by Spearman correlation analysis. ANOVA was used for comparisons among more than two groups, followed by Dunnett's multiple-comparison test to identify differences. Significance was defined as a P value of <0.05. The statistical software used was Prism version 8 (GraphPad Software, Inc.).

Single-cell transcriptome atlas of DCIS and IDC

We initially performed single-cell transcriptome profiling on 7 DCIS and 6 IDC samples and 1 normal human breast tissue sample (Fig. 1A). Selected DCIS sample was also examined by spatial transcriptome analysis. Normal breast tissues were obtained from a surgical specimen collected from a primary breast cancer patient. Pathohistological information is shown in Fig. 1B; Supplementary Table S2. A summary of the sequencing results is shown in Supplementary Table S3. During quality control with the Seurat R package, low-quality sequencing data were computationally excluded, and to reduce the batch effect of each sample, the data were integrated with “Harmony” (26, 27). For removing putative doublet cells, we applied “DoubletFinder” (Supplementary Fig. S1A and S1B; ref. 25). Unsupervised clustering followed by UMAP revealed 32 distinct clusters (Supplementary Fig. S2A). According to the gene expression profiles of the cell type-specific markers, this cluster analysis broadly segregated the cells into 12 groups: Luminal epithelial cells, proliferating luminal epithelial cells, basal cells, proliferating basal cells, T cells, proliferating immune cells, B cells, plasma cells, macrophages, monocytes, erythrocytes and stromal cells (Fig. 1C, left side). The proportions of cells within the samples (Fig. 1D) and within the cell types (Supplementary Fig. S2B) showed cellular diversity in the DCIS and IDC samples. Hormone receptor genes such as ESR1 and PGR were detected selectively in epithelial cell clusters (Fig. 1E). In addition, the expression of other cell type–specific markers was confirmed by UMAP (Fig. 1E and F, Supplementary Fig. S2C and Supplementary Table S4). We identified myoepithelial cells that are ACTA2-positive cells (Fig. 1F), mixed in the basal cell cluster.

To further segregate the clusters of epithelial and basal cells in detail, “luminal epithelial cells,” “proliferating luminal epithelial cells,” “basal cells,” and “proliferating basal cells” were isolated and remapped on a separate UMAP plot (Fig. 1C, right side). Thus, we next used the clusters of luminal cells and proliferating luminal cells to analyze DCIS diversity and to search for genes responsible for invasiveness in IDC.

Variations in gene expression among DCIS and IDC samples at single-cell resolution.

To deeply investigate epithelial heterogeneity in luminal cells, a total of 30,808 cells from the clusters of “luminal cells” and “proliferating luminal cells” were re-clustered (Fig. 2A). The UMAP plot with the luminal cell lineages revealed 11 distinct clusters. We examined the distribution of the cells among the disease states (Fig. 2B and C) as well as among the patients (Supplementary Fig. S3A). In the 11 clusters, we observed a nonproportional distribution of DCIS and IDC cells. DCIS cells were found in all clusters excluding Cluster 10. Clusters 2, 7, and 10 predominantly contained IDC cells. In addition, luminal cells from the normal breast tissue were mainly found in Cluster 5. The correlation plots among the clusters also showed that Cluster 5, which mainly contained normal luminal cells, was clearly distinct from the other clusters (Fig. 2D). These data suggest luminal cell diversity among DCIS and IDC cells. The selected differentially expressed genes in each of the clusters are summarized (Fig. 2E) and exhibit the diversity of luminal cell lineages. We further investigated cluster-specific gene signatures by performing pathway enrichment analysis on the luminal cell clusters (Fig. 2F and Supplementary Table S5). We found clear differences in pathway enrichment between 11 clusters. For example, Clusters 1, 2, 3, and 6 were associated with immune-related pathways. The term “positive regulation of cell cycle” was significantly enriched in Cluster 4, reminiscent of proliferating luminal cells. In addition, module analysis of the cell-cycle confirmed increased cell proliferation in Cluster 4 (Supplementary Fig. S3B).

Identification of genes presumably related to the invasive phenotype

We sought to identify genes associated with the transition from DCIS to IDC. We compared the gene expression between all DCIS cells and all IDC cells. As shown in Fig. 3A, some clusters are DCIS or IDC specific. The majority of cells in Clusters 2 and 7 were IDC cells, with a few DCIS cells, and Cluster 10 contained only IDC cells. When differentially expressed genes between DCIS and IDC were analyzed in all cells, we found significantly upregulated genes in IDC samples (Supplementary Tables S6 and S7). Six genes on the gene list—RPL17, RNASEK, MIF, CRIP1, CEACAM6, and RPL36A—were found to be gradually upregulated from DCIS to IDC (Fig. 3B). Specifically expressed genes in IDC were shown (Fig. 3C), which were not detected in DCIS and normal breast samples. In addition, we found that CPB1 expression was the highest in Cluster 7, which predominantly contained IDC (Supplementary Table S8). As the examples, the UMAP plots for RPL17, MIF, SNHG29, and SNHG6 showed their expression patterns at the single-cell level (Fig. 3D and E). Higher expression levels of RPL17 and MIF were observed in IDC samples, although lower expression levels were detected in DCIS and normal breast samples (Fig. 3D). Specific expression of the long noncoding RNA (lncRNA) SNHG6 and SNHG29 was observed in IDC samples (Fig. 3E). Because some of the selected genes, such as SNHG6 and SNHG29, have already been reported to be associated with tumor malignancy (36, 37), we selected other genes with unknown functions in cancer progression and investigated their functions. As shown in Supplementary Fig. S4A–S4F, silencing of RPL17 or CPB1 decreased cell growth. In particular, knockdown of CPB1 with specific siRNAs resulted in morphological changes from a spindle-shaped to a cobblestone-like morphology in MCF7 cells (Supplementary Fig. S4E, bottom). Although further investigation is required, these genes might be related to the invasive phenotype of breast cancer.

We also analyzed the upregulated genes in DCIS cells (Supplementary Table S9). TFF1, ANKRD30A, TAT, and MDK were more highly expressed in DCIS than in IDC. In addition, we compared gene expression by grade in DCIS. We defined high and low-grade according to the pathological diagnosis, and then found that high-grade DCIS cells were enriched in Clusters 8, 9, and 11 (Supplementary Fig. S5A). To figure out markers for high-grade DCIS cells, we focused on these clusters (Supplementary Fig. S5B) and identified genes whose expression gradually elevated in the order of low-grade DCIS cells, high-grade DCIS cells, and IDC (Supplementary Fig. S5C). Although it is a comparison with a small number of cells, some of these genes are reported to associate with cancer progression.

HER2 amplification and heterogeneity in DCIS

As shown in Supplementary Table S2, 4 of the 7 DCIS patients were diagnosed with HER2+ luminal DCIS. IHC showed an obvious difference in the staining intensity and percentage of HER2 (Fig. 4A) To investigate the molecular features of HER2+ luminal DCIS, we compared the HER2+ luminal DCIS samples (NCCBC5, NCCBC6, NCCBC11, and NCCBC14) with the HER2 luminal DCIS samples (NCCBC2, NCCBC3, and NCCBC7) in the UMAP plot with HER2 states (Fig. 4B). The cell proportions between HER2+ and HER2 luminal DCIS samples are shown based on the clusters, as established in Fig. 2A. Because Cluster 10 contained only IDC cells, DCIS cells were not detected (Fig. 4C). The cell distribution was also examined on the basis of the clusters (Fig. 4D) and patient (Supplementary Fig. S6A). The UMAP plot with HER2 states revealed a variation between HER2+ and HER2 luminal DCIS (Fig. 4B). Because Clusters 8 and 9 were primarily comprised of HER2+ luminal DCIS cells (Fig. 4C), we investigated highly expressed genes in these clusters and identified 24 genes as selective markers. Interestingly, 11 of these 24 genes are located at the chromosome 17q12-q21 locus (Fig. 4E; Supplementary Fig. S6B), which is a well-known HER2 amplicon in breast cancer (38). These data implied that HER2 amplification may be induced at an early stage of luminal breast cancer. As the data shown, these genes were highly expressed in the cells in Clusters 8 and 9 compared with the cells in the other clusters or HER2 luminal DCIS cells (Supplementary Fig. S6C). In HER2+ luminal DCIS cells, these genes were not uniformly expressed, suggesting cellular heterogeneity among HER2+ luminal DCIS. We also compared “HER2 DCIS versus HER2 IDC” and “HER2+ DCIS versus HER2+ IDC,” respectively (Supplementary Table S10–S13).

To further investigate the cellular heterogeneity of HER2+ luminal DCIS, we performed spatial transcriptome analysis of NCCBC6. For sequencing library construction, frozen sections of DCIS tissue were applied to barcoded slides (Supplementary Fig. S7A). After sequencing, transcriptome data from 1,521 distinct spots were obtained. The transcriptome data for all the spots were analyzed by UMAP (Supplementary Fig. S7B). The spatial information provided by the barcodes was used to visualize the tissue context of gene expression (Supplementary Fig. S7C). The expression levels of CD24 and KRT18, luminal epithelial markers, indicated that Clusters 3 and 4 were composed of luminal DCIS cells (Supplementary Fig. S7D). Although ERBB2/HER2 expression was observed in Clusters 3 and 4, Cluster 4 showed a higher expression level of ERBB2 (Supplementary Fig. S7E and S7F). To see the difference between Cluster 3 (HER2-low) and Cluster 4 (HER2-high), we analyzed differentially expressed genes between these clusters (244 genes, fold change > 1.1 and P < 0.05; Supplementary Table S14). Unsupervised clustering analysis showed clear separation of the cells from the two clusters (Supplementary Fig. S7G). Analysis of canonical pathways and upstream regulators by IPA indicated functionally distinct luminal epithelial populations in one HER2+ DCIS sample (Supplementary Fig. S7H and S7I). These data also confirmed the intratumoral heterogeneity in DCIS, particularly in HER2+ luminal DCIS samples.

High heterogeneity in DCIS luminal epithelial cells

The spatial transcriptome data of HER2+ luminal DCIS prompted us to further investigate cellular heterogeneity in DCIS. We recently developed a method to quantify and visualize cellular diversity based on gene expression profiles (30). In this method, closeness centrality is calculated to show alterations in gene expression in an arbitrary population (see Materials and Methods). Higher centrality indicates lower heterogeneity, and conversely, lower centrality indicates higher heterogeneity (Fig. 5A). Unexpectedly, closeness centrality analysis of luminal cells from DCIS and IDC samples revealed lower centrality in DCIS samples, indicating higher heterogeneity in DCIS than in IDC (Fig. 5B). To exclude the possibility that the large individual differences between the DCIS samples might influence the result, we calculated closeness centrality in each sample (Supplementary Fig. S8). Although variations in closeness centrality were observed among the 7 DCIS and 6 IDC samples, many DCIS samples showed low centrality, and this result at least suggested that the heterogeneity in luminal cells in DCIS samples was similar to or greater than that in IDC samples.

Next, we examined closeness centrality in microenvironment cells in DCIS and IDC. As shown in Fig. 5C, the population of endothelial cells and fibroblasts showed lower centrality in IDC than DCIS, suggesting that endothelial cells and fibroblasts were more heterogeneous in IDC, reminiscent of the emergence of tumor endothelial cells and cancer-associated fibroblasts. Analysis of immune cell types showed lower centrality in most lymphocytes and macrophages except plasma cells and neutrophils, although almost no differences were observed in erythrocytes (Fig. 5D). These data also confirmed the diversity of cells in the tumor immune microenvironment, particularly within the lymphatic and macrophage lineages.

Differences CCIs between luminal epithelial cells and immune cells in DCIS

Although our single-cell isolation protocol could yield only a few populations of fibroblasts and endothelial cells, we collected a sufficient amount of single-cell transcriptome data from immune cells (Fig. 1C and D). The dataset comprised a number of immune cell-related clusters, representing 24.5% (11,166 cells) of the cells in the human breast tissues, including the DCIS and IDC samples. To investigate the relationship between the immune cell population and breast cancer cells, single-cell data of immune clusters based on UMAP plot (Fig. 1C) were selected, and a UMAP plot with all immune cells was generated (Fig. 6A; Supplementary Fig. S9A). The UMAP plot of immune-related cells revealed 11 distinct clusters, which were categorized by marker gene expression patterns (Supplementary Table S15). The cellular composition of immune-related cells showed differences in the cell number between the DCIS and IDC samples (Fig. 6B) as well as normal breast samples (Supplementary Fig. S9B and S9C). Although the total number was small, CD4+ T cells were only detected in IDC samples. Comparative analysis of the cell number in the clusters among the disease states revealed that the number of macrophages tended to increase in the DCIS samples (Fig. 6C, left). In contrast, the number of CD8+ T cells tended to increase in the IDC samples (Fig. 6C, right).

To understand interactions between DCIS or IDC cells (luminal epithelial cells) and immune cells, we calculated CCI scores in both the DCIS and IDC samples. For this purpose, we focused on immunological responses—we extracted chemokine and inflammatory cytokine data from a publicly available database, the CCI Database (33), and established an immune response network between luminal cells and immune cells (30). As described previously in Materials and Methods, the interaction score between luminal cells and immune cells was computed from the expression levels of the two proteins of interest and the positive interaction rate in each cluster. We identified 720 predicted interactions (479 pairs of ligands in the epithelial cluster and receptors in the immune cluster and 241 pairs of ligands in the immune cluster and receptors in the epithelial cluster) in the DCIS samples and 719 predicted interactions (480 pairs of ligands in the epithelial cluster and receptors in the immune cluster and 239 pairs of ligands in the immune cluster and receptors in the epithelial cluster) in the IDC samples. To directly compare the interactions between immune cells and luminal cells in DCIS and IDC, we plotted the scores of the CCIs (between ligands in the epithelial cluster and receptors in the immune cluster, Fig. 6D). Greater interaction was observed in the IDC samples; however, luminal cells in the DCIS samples were also predicted to interact with immune cells at levels comparable with those in the IDC samples. In the IDC samples, stronger interactions with clusters of neutrophils and CD4+ T cells were predicted. In addition, luminal cells in the DCIS samples were predicted to interact with immune cells such as NK/NKT cells and Treg. When focusing on each interaction between luminal cells and immune cells with individual interaction scores (CCI score >0.15), the data showed apparent similarities and differences in the interactomes in the DCIS and IDC samples (Fig. 6E).

In this study, we investigated intratumoral heterogeneity in DCIS and the molecular features at single-cell resolution. Accumulating evidence has demonstrated that the intratumoral heterogeneity found in IDC might already be established in advanced DCIS lesions (1, 39–41). Previous reports suggested that DCIS consists of various types of epithelial cells, including proliferating cells, immune response-related cells, and estrogen response-related cells (21, 22), which were also found in our scRNA-seq. Spatial transcriptome analysis also confirmed that these clusters were heterogeneously localized in each duct (42). Our data showed that DCIS samples exhibited higher intratumoral heterogeneity than IDC samples and contained cells with gene expression profiles similar to those of IDC cells. On the other hand, our scRNA-seq analysis identified DCIS-dominant clusters and IDC-dominant clusters. These differences in clusters between DCIS and IDC could not be detected by bulk RNA sequencing. When analyzing differentially expressed genes between the DCIS and IDC samples, we found significantly upregulated genes in the IDC samples compared with the DCIS and normal breast tissue samples. Some of these transcripts were well-known lncRNAs, such as SNHG6 and SNHG29. LncRNAs have been associated with diverse functions in cancer development (43). For example, SNHG6 is reported to be related to cell proliferation and migration in high-grade PR+ breast cancer by binding to miR-26a (36). SNHG29 is known to promote cell proliferation and invasion by regulating KLF4 and miR-380-3p in triple-negative breast cancer (37). In addition, a loss-of-function experiment with siRNA transfection revealed that CPB1 and RPL17 contributed to cancer progression. CPB1 has been reported to be related to metastasis in low-grade luminal A breast cancer (44). CPB1 is a secretory tissue protease, and it may promote invasion in relation to matrix metalloprotease activity (45). CPB1 and SNHG6 may be associated with poor prognosis in low-grade breast cancer (44, 46). These genes may have potential as biomarkers for a high recurrence risk in DCIS and as therapeutic targets.

We also identified HER2+ luminal DCIS as having cellular heterogeneity. These samples contained cells ununiformly overexpressing ERBB2 and neighboring genes. HER2 amplification is an oncogenic event, and it has already occurred at the DCIS stage (47). Our data also implied that HER2 amplification may be induced at an early-stage of HER2+ luminal breast cancer. Previous reports explained that distinct CNV events (1q gain and 8p loss) occurred before HER2 amplification in DCIS and that the diversity of CNVs was increased after HER2 amplification (47). Moreover, HER2 amplification activates the interferon signaling pathway (47). HER2+ luminal DCIS has high heterogeneity, and it has been reported that its prognosis is poor (23). HER2 positivity in DCIS has been reported to be associated with a risk for local recurrence but not for infiltration (48). Thus, HER2 positivity in luminal DCIS may be an important prognostic marker.

Intratumoral heterogeneity is currently considered to contribute to cancer aggressiveness and resistance to anticancer drugs (49). In DCIS, high intratumoral heterogeneity is associated with a high-proliferative index, poor tumor differentiation, necrosis, increased immunostaining for p53, and increased extent of disease (50, 51). Moreover, recent single-cell RNA-seq data suggest that increasing the diversity of stromal cells (22, 52) and immune cells (22, 42, 53) is associated with invasiveness. Also, the myoepithelial cell phenotype and morphology were shown to associate with invasiveness (24). Previous invasive breast cancer tissue single-cell resolution identified the heterogeneity of stromal cells such as three endothelial states, three perivascular-like states, and five CAF states. CAFs include inflammatory CAFs and myofibroblast-like CAFs. Myofibroblast-like CAFs are more enriched in invasive cancer regions than DCIS, which correlates with immunosuppressive environments that are poor in CD8+ T lymphocytes and enriched in CD4+ T lymphocytes expressing high levels of immune checkpoints such as PD-1 and CTLA4 (22, 52). It is also reported that T cells in IDC are highly heterogeneous, enriched in Treg, and express high levels of PD-L1 (53). In our study, we confirmed the heterogeneity of stromal cells and immune cells increased in IDC than in DCIS. Moreover, CD8 T cells and Treg increased in IDC, and CD4+ T cells were identified in only IDC. The immunosuppressive environment may relate to invasiveness.

This study has some limitations. First, the cohort size was limited to 14 patients (7 DCIS samples, 6 IDC samples, and 1 normal breast tissue sample), and there is a possibility that some of these tumors have a low risk of invasion, especially the low-grade DCIS tumors. One of the solutions for this issue is large-scale scRNA-seq analysis. Second, in our single-cell dissociation process, we obtained a small number of fibroblasts and endothelial cells, which might lead to sampling bias. This is more likely because of our digestion protocol. Digestion condition is reported to affect the proportion of cell type (54). With our single-cell isolation protocol, we could isolate epithelial cells and immune cells efficiently; however, endothelial and fibroblastic cells were not so effectively collected. Third, there is a possibility of admixture of normal and benign in DCIS and IDC cases. To reduce the potential normal cell contamination from DCIS and IDC cases, we performed scRNA-seq for normal breast tissue that was isolated from a breast cancer specimen but pathologically normal, and normal luminal cell cluster was excluded in the analysis. However, we assume that this attempt still cannot rule out the possibility of admixture of normal cell types in the analysis. To identify normal cells, there are some approaches to infer single-cell CNV profiles from scRNA-seq data: InferCNV (55), CaSpER (56), and CopyKAT (57). Although these approaches are still developing, combining with CNV data could realize a better performance to predict normal cells from tumor cells in the scRNA-seq analysis. Forth, differences among samples in both cell number and sequence depth (Supplementary Table S3) may influence the analysis. To solve this issue, we applied “Harmony” for fair integration of scRNA-seq data; however, further attempt could be required for better data investigation in the process of technological and analytical advances. Fifth, in this study, we collected each sample from different patients. Sampling of paired DCIS and IDC tissues from the same patients might be an optimal to investigate the differences between DCIS and IDC and to identify the invasiveness-related genes. For better understanding of molecular mechanism of the transition from DCIS to IDC, further investigation would be needed.

In summary, we generated single-cell transcriptome profiles of DCIS and identified the genes that may be related to the invasive phenotype in breast cancer. This study not only provides new insight into breast cancer biology but also contributes to the identification of candidate genes for clinically usable diagnostic and therapeutic markers in luminal DCIS patients. These biomarkers may be useful for the identification of patients who should undergo resection or more intensive therapy.

T. Kohno reports personal fees from Eli Lilly Japan, and grants from Chugai Pharmaceutical and Sysmex outside the submitted work. No disclosures were reported by the other authors.

M. Tokura: Conceptualization, data curation, software, investigation, writing–original draft. J. Nakayama: Conceptualization, data curation, software, formal analysis, supervision, validation, investigation, visualization, methodology, writing–original draft. M. Prieto-Vila: Investigation. S. Shiino: Sample collection. M. Yoshida: Sample collection. T. Yamamoto: Investigation. N. Watanabe: Investigation. S. Takayama: Sample collection. Y. Suzuki: Data curation, methodology. K. Okamoto: Software, investigation. T. Ochiya: Supervision. T. Kohno: Supervision. Y. Yatabe: Supervision, sample collection. A. Suto: Supervision, sample collection. Y. Yamamoto: Conceptualization, data curation, software, formal analysis, supervision, funding acquisition, visualization, methodology, writing–original draft, project administration.

The authors appreciate the contribution of Drs. K. Jinbo, T. Murata, and C. Watase for supporting sample collection and Dr. Yutaro Mori for technical support. The present work was supported in part by a Grant-in-Aid for Scientific Research (B), Grant Number: 21H02721 (to Y. Yamamoto), Grant-in-Aid for Early-Career Scientists, Grant Number: 21K15562, Grant-in-Aid for JSPS Fellows: 20J01794, and JSPS KAKENHI Grant Number 16H06279(PAGS) (to J. Nakayama), research grant from SGH Foundation (to Y. Yamamoto), and MSD Life Science Foundation (to Y. Yamamoto). Schematics in figures were made using an academic license of Biorender software.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Cowell
CF
,
Weigelt
B
,
Sakr
RA
,
Ng
CKY
,
Hicks
J
,
King
TA
, et al
.
Progression from ductal carcinoma in situ to invasive breast cancer: revisited
.
Mol. Oncol
2013
;
7
:
1
11
.
3.
Collins
LC
,
Tamimi
RM
,
Baer
HJ
,
Connolly
JL
,
Colditz
GA
,
Schnitt
SJ
Outcome of patients with ductal carcinoma in situ untreated after diagnostic biopsy: results from the nurses’ health study
.
Cancer
2005
;
103
:
1778
84
.
4.
Betsill
WL
Jr
,
Rosen
PP
,
Lieberman
PH
,
Robbins
GF
Intra- ductal carcinoma. Long-term follow-up after treatment by biopsy alone
.
JAMA
1978
;
239
:
1863
7
.
5.
Eusebi
V
,
Feudale
E
,
Foschini
MP
,
Micheli
A
,
Conti
A
,
Riva
C
, et al
.
Long-term follow-up of in situ carcinoma of the breast
.
Semin Diagn Pathol
1994
;
11
:
223
35
.
6.
Erbas
B
,
Provenzano
E
,
Armes
J
,
Gertig
D
The natural history of ductal carcinoma in situ of the breast: a review
.
Breast Cancer Res Treat
2006
;
97
:
135
44
.
7.
Page
DL
,
Dupont
WD
,
Rogers
LW
,
Landenberger
M
Intraductal carcinoma of the breast: follow-up after biopsy only
.
Cancer
1982
;
49
:
751
8
.
8.
Ryser
MD
,
Weaver
DL
,
Zhao
F
,
Worni
M
,
Grimm
LJ
,
Gulati
R
, et al
.
Cancer outcomes in DCIS patients without locoregional treatment
.
J Natl Cancer Inst
2019
;
111
:
952
60
.
9.
Sagara
Y
,
Mallory
MA
,
Wong
S
,
Aydogan
F
,
DeSantis
S
,
Barry
WT
, et al
.
Survival benefit of breast surgery for low-grade ductal carcinoma in situ: a population-based Cohort Study
.
JAMA Surg
2015
;
150
:
739
45
.
10.
Maxwell
AJ
,
Clements
K
,
Hilton
B
,
Dodwell
DJ
,
Evans
A
,
Kearins
O
, et al
.
Risk factors for the development of invasive cancer in unresected ductal carcinoma in situ
.
Eur J Surg Oncol
2018
;
44
:
429
35
.
11.
Parikh
U
,
Chhor
CM
,
Mercado
CL
Ductal carcinoma in situ: the whole truth
.
AJR Am J Roentgenol
2018
;
210
:
246
55
.
12.
Vargas
AC
,
Reed
AEM
,
Waddell
N
,
Lane
A
,
Reid
LE
,
Smart
CE
, et al
.
Gene expression profiling of tumour epithelial and stromal compartments during breast cancer progression
.
Breast Cancer Res Treat
2012
;
135
:
153
65
.
13.
Casasent
A
,
Edgerton
M
,
Navin
NE
Genome evolution in ductal carcinoma in situ: invasion of the clones
.
J Pathol
2017
;
241
:
208
18
.
14.
Strand
SH
,
Rivero-Gutiérrez
B
,
Houlahan
KE
,
Seoane
JA
,
King
L
,
Risom
T
, et al
.
DCIS genomic signatures define biology and clinical outcome: Human Tumor Atlas Network (HTAN) analysis of TBCRC 038 and RAHBT cohorts
.
bioRxiv
2021
.
DOI: 10.1101/2021.06.16.448585
.
15.
Trinh
A
,
Alcazar
CRGD
,
Shukla
SA
,
Chin
K
,
Chang
YH
,
Thibault
G
, et al
.
Genomic alterations during the in situ to invasive ductal breast carcinoma transition shaped by the immune system
.
Mol Cancer Res
2021
;
19
:
623
35
.
16.
Lesurf
R
,
Aure
MR
,
Mork
HH
,
Vitelli
V
,
Oslo Breast Cancer Research Consortium(OSBREAC)
,
Lundgren
S
, et al
.
Molecular features of subtype-specific progression from ductal carcinoma in situ to invasive breast cancer
.
Cell Rep
2016
;
16
:
1166
79
.
17.
Regan
JL
,
Smalley
MJ
Integrating single-cell RNA-sequencing and functional assays to decipher mammary cell states and lineage hierarchies
.
npj Breast Cancer
2020
;
6
:
32
.
18.
Nguyen
QH
,
Pervolarakis
N
,
Blake
K
,
Ma
D
,
Davis
RT
,
James
N
, et al
.
Profiling human breast epithelial cells using single-cell RNA sequencing identifies cell diversity
.
Nat Commun
2018
;
9
:
2028
.
19.
Karaayvaz
M
,
Cristea
S
,
Gillespie
SM
,
Patel
AP
,
Mylvaganam
R
,
Luo
CC
, et al
.
Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq
.
Nat Commun
2018
;
9
:
3588
.
20.
Azizi
E
,
Carr
AJ
,
Plitas
G
,
Cornish
AE
,
Konopacki
C
,
Prabhakaran
S
, et al
.
Single-cell map of diverse immune phenotypes in the breast tumor microenvironment
.
Cell
2018
;
174
:
1293
308
.
21.
Pal
B
,
Chen
Y
,
Vaillant
F
,
Capaldo
BD
,
Joyce
R
,
Song
X
, et al
.
A single-cell RNA expression atlas of normal preneoplastic and tumorigenic states in the human breast
.
EMBO J
2021
;
40
:
e107333
.
22.
Wu
SZ
,
Al-Eryani
G
,
Roden
DL
,
Junankar
S
,
Harvey
K
,
Andersson
A
, et al
.
A single-cell and spatially resolved atlas of human breast cancers
.
Nat Genet
2021
;
53
:
1334
47
.
23.
Nagasawa
S
,
Kuze
Y
,
Maeda
I
,
Kojima
Y
,
Motoyoshi
A
,
Onishi
T
, et al
.
Genomic profiling reveals heterogeneous populations of ductal carcinoma in situ of the breast
.
Commun Biol
2021
;
4
:
438
.
24.
Risom
T
,
Glass
DR
,
Averbukh
I
,
Liu
CC
,
Baranski
A
,
Kagel
A
, et al
.
Transition to invasive breast cancer is associated with progressive changes in the structure and composition of tumor stroma
.
Cell
2022
;
185
:
299
310
.
25.
McGinnis
CS
,
Murrow
LM
,
Gartner
ZJ
DoubletFinder: Doublet detection in single-cell RNA sequencing data using artificial nearest neighbors
.
Cell Syst
2019
;
8
:
329
37
.
26.
Korsunsky
I
,
Millard
N
,
Fan
J
,
Slowikowski
K
,
Zhang
F
,
Wei
K
, et al
.
Fast, sensitive and accurate integration of single-cell data with Harmony
.
Nat Methods
2019
;
16
:
1289
96
.
27.
Tran
HTN
,
Ang
KS
,
Chevrier
M
,
Zhang
X
,
Lee
NYS
,
Goh
M
, et al
.
A benchmark of batch-effect correction methods for single-cell RNA sequencing data
.
Genome Biol
2020
;
21
:
12
.
28.
Yu
G
,
Wang
LG
,
Han
Y
,
He
QY
clusterprofile: an R package for comparing biological themes among gene clusters
.
OMICS
2012
;
16
:
284
7
.
29.
Nakayama
J
,
Ito
E
,
Fujimoto
J
,
Watanabe
S
,
Semba
K
Comparative analysis of gene regulatory networks of highly metastatic breast cancer cells established by orthotopic transplantation and intra-circulation injection
.
Int J Oncol
2017
;
50
:
497
504
.
30.
Watanabe
N
,
Nakayama
J
,
Fujita
Y
,
Mori
Y
,
Kadota
T
,
Shimomura
I
, et al
.
Anomalous epithelial variations and ectopic inflammatory response in chronic obstructive pulmonary disease
.
medRxiv
2020
.
doi: https://doi.org/10.1101/2020.12.03.20242412
.
31.
Armingol
E
,
Officer
A
,
Harismendy
O
,
Lewis
NE
Deciphering cell–cell interactions and communication from gene expression
.
Nat Rev Genet
2021
;
22
:
71
88
.
32.
Kumar
MP
,
Du
J
,
Lagoudas
G
,
Jiao
Y
,
Sawyer
A
,
Drummond
DC
, et al
.
Analysis of single-cell RNA-seq identifies cell–cell communication associated with tumor characteristics
.
Cell Rep
2018
;
25
:
1458
68
.
33.
Wang
Y
,
Wang
R
,
Zhang
S
,
Song
S
,
Jiang
C
,
Han
G
, et al
.
iTALK: an R package to characterize and illustrate intercellular communication
.
bioRxiv
2019
;
507871
.
34.
Gu
P
,
Chen
H
Modern bioinformatics meets traditional chinese medicine
.
Brief Bioinform
2014
;
15
:
984
1003
.
35.
Hafemeister
C
,
Satija
R
Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression
.
Genome Biol
2019
;
20
:
296
.
36.
Wang
H
,
Zhang
W
,
Zhu
H
,
Li
Q
,
Miao
L
Long noncording RNA SNHG6 mainly functions as a competing endogenous RNA in human tumors
.
Cancer Cell Int
2020
;
20
:
219
.
37.
Li
S
,
Wu
D
,
Jia
H
,
Zhang
Z
Long non-coding RNA LRRC75A-AS1 facilitates triple negative breast cancer cell proliferation and invasion via functioning as a ceRNA to modulate BAALC
.
Cell Death Dis
2020
;
11
:
643
.
38.
Sahlberg
KK
,
Hongisto
V
,
Edgren
H
,
Mäkelä
R
,
Hellström
K
,
Due
EU
, et al
.
The HER2 amplicon includes several genes required for the growth and survival of HER2 positive breast cancer cells
.
Mol Oncol
2013
;
7
:
392
401
.
39.
Buerger
H
,
Otterbach
F
,
Simon
R
,
Poremba
C
,
Diallo
R
,
Decker
T
, et al
.
Comparative genomic hybridization of ductal carcinoma in situ of the breast-evidence of multiple genetic pathways
.
J Pathol
1999
;
187
:
396
402
.
40.
Ma
XJ
,
Salunga
R
,
Tuggle
JT
,
Gaudet
J
,
Enright
E
,
McQuary
P
, et al
.
Gene expression profiles of human breast cancer progression
.
Proc Natl Acad Sci U S A
2003
;
100
:
5974
9
.
41.
Vincent-Salomon
A
,
Lucchesi
C
,
Gruel
N
,
Raynal
V
,
Pierron
G
,
Goudefroye
R
, et al
.
Integrated genomic and transcriptomic analysis of ductal carcinoma in situ of the breast
.
Clin Cancer Res
2008
;
14
:
1956
65
.
42.
Wei
R
,
He
S
,
Bai
S
,
Sei
E
,
Hu
M
,
Thompson
A
, et al
.
Spatial charting of single-cell transcriptomes in tissues
.
Nat Biotechnol
2022
.
43.
Sanchez Calle
A
,
Kawamura
Y
,
Yamamoto
Y
,
Takeshita
F
,
Ochiya
T
Emerging roles of long non-coding RNA in cancer
.
Cancer Sci
2018
;
109
:
2093
100
.
44.
Bouchal
P
,
Dvořáková
M
,
Roumeliotis
T
,
Bortlíček
Z
,
Ihnatová
I
,
Procházková
I
, et al
.
Combined proteomics and transcriptomics identifies carboxypeptidase B1 and nuclear factorκB (NF-κB)–associated proteins as putative biomarkers of metastasis in low grade breast cancer
.
Mol Cell Proteomics
2015
;
14
:
1814
30
.
45.
Swaisgood
CM
,
Schmitt
D
,
Eaton
D
,
Plow
EF
In vivo regulation of plasminogen function by plasma carboxypeptidase B
.
J Clin Invest
2002
;
110
:
1275
82
.
46.
Ding
X
,
Zhang
Y
,
Yang
H
,
Mao
W
,
Chen
B
,
Yang
S
, et al
.
Long non-coding RNAs may serve as biomarkers in breast cancer combined with primary lung cancer
.
Oncotarget
2017
;
35
:
58210
21
.
47.
Lu
P
,
Foley
J
,
Zhu
C
,
McNamara
K
,
Sirinukunwattana
K
,
Vennam
S
, et al
.
Transcriptome and genome evolution during HER2-amplified breast neoplasia
.
Breast Cancer Res
2021
;
23
:
73
.
48.
Thorat
MA
,
Levey
PM
,
Jones
JL
,
Pinder
SE
,
Bundred
NJ
,
Fentiman
IS
, et al
.
Prognostic and predictive value of HER2 expression in ductal carcinoma in situ: results from the UK/ANZ DCIS randomized trial
.
Clin Cancer Res
2021
;
27
:
5317
24
.
49.
Prieto-Vila
M
,
Usuba
W
,
Takahashi
R
,
Shimomura
I
,
Sasaki
H
,
Ochiya
T
, et al
.
Single-cell analysis reveals a preexisting drug-resistant subpopulation in the luminal breast subtype
.
Cancer Res
2019
;
79
:
4412
25
.
50.
Allred
DC
,
Wu
Y
,
Mao
S
,
Nagtegaal
ID
,
Lee
S
,
Perou
CM
, et al
.
Ductal carcinoma in situ and the emergence of diversity during breast cancer evolution
.
Clin Cancer Res
2008
;
14
:
370
8
.
51.
Park
SY
,
Gonen
M
,
Kim
HJ
,
Michor
F
,
Polyak
K
Cellular and genetic diversity in the progression of in situ human breast carcinomas to an invasive phenotype
.
J Clin Invest
2010
;
120
:
636
44
.
52.
Kieffer
Y
,
Hocine
HR
,
Gentric
G
,
Pelon
F
,
Bernard
C
,
Bourachot
B
, et al
.
Single-cell analysis reveals fibroblast clusters linked to immunotherapy resistance in cancer
.
Cancer Discov
2020
;
10
:
1330
51
.
53.
Alcazar
CRGD
,
Huh
SJ
,
Ekram
MB
,
Trinh
A
,
Liu
LL
,
Beca
F
, et al
.
Immune escape in breast cancer during in situ to invasive carcinoma transition
.
Cancer Discov
2017
;
7
:
1098
115
.
54.
Denisenko
E
,
Guo
BB
,
Jones
M
,
Hou
R
,
Kock
L
,
Lassmann
T
, et al
.
Systematic assessment of tissue dissociation and storage biases in single-cell and single-nucleus RNA-seq workflows
.
Genome Biol
2020
;
21
:
130
.
55.
Patel
AP
,
Tirosh
I
,
Trombetta
JJ
,
Shalek
AK
,
Gillespie
SM
,
Wakimoto
H
, et al
.
Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma
.
Science
2014
;
344
:
1396
401
.
56.
Harmanci
AS
,
Harmanci
AO
,
Zhou
X
.
CaSpER identifies and visualizes CNV events by integrative analysis of single-cell or bulk RNA-sequencing data
.
Nat Commun
2020
;
11
:
89
.
57.
Gao
R
,
Bai
S
,
Henderson
YC
,
Lin
Y
,
Schalck
A
,
Yan
Y
, et al
.
Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes
.
Nat Biotechnol
2021
;
39
:
599
608
.

Supplementary data