Genomic and transcriptomic correlates of immunotherapy response within the tumor microenvironment of leptomeningeal metastases

Study design: patients

These clinical trials (Clinicaltrials.gov identifier NCT02886585 and NCT02939300) were designed by the principal investigators and the Dana-Farber Harvard Cancer Center (DF/HCC) Institutional review board approved the protocol. The study was designed by the principal investigators and conducted in accordance with the provision of the Declaration of Helsinki and Good Clinical Practice guidelines. The Dana-Farber/Harvard Cancer Center institutional review board approved the protocol. All patients provided informed consent. Eligible patients had histologically confirmed disease from any solid tumor and LMD defined by positive CSF cytology for malignant cells. Other key inclusion criteria included the following: an ECOG performance status ≤2, normal organ and marrow function, and a stable dose of dexamethasone of 2 mg or less for 7 days prior to initiation of treatment. Given the frequent occurrence of neurologic symptoms (e.g. headaches) associated with LMD, 11/20 patients included in this study were on a steroid regimen at the time of enrollment and 10/20 patients were treated with steroids while receiving ICI¹³. Written informed consent was obtained for all participants.

Further details of the subjects’ clinical courses, including cytology, steroid dosage, and Ommaya shunt status, are provided in Supplementary Data 1.

Study design: treatment and endpoints

These studies were designed as open-label, single arm Phase-II clinical trials to evaluate ICI in patients with LMD of any histology. For the first trial (NCT02886585), patients with CNS metastases were enrolled across multiple cohorts. Cohorts A, B, and D include patients with parenchymal brain metastases. The LMD cohort was Cohort C of this study. Pembrolizumab was administered intravenously at 200 mg every 3 weeks until disease progression, death, or grade 3–4 toxicity. A brain MRI and CT chest/abdomen/pelvis were obtained every 6 weeks for restaging purposes. The primary endpoint of the LMD cohort was the rate of overall survival at 3 months (OS3). 11 patients with LMD enrolled to NCT02886585 were included for single-cell analysis; 4 of whom had serum and CSF sampling at multiple time points during treatment. All 11 patients had metastatic breast cancer (Supplementary Data 1).

Treatment in the second trial (NCT02939300) consisted of ipilimumab and nivolumab. Ipilimumab and nivolumab were administered intravenously every 3 weeks for 4 doses. Afterwards, ipilimumab was given every 6 weeks and nivolumab was given every 2 weeks (for non-small cell lung cancer and head and neck cancer) or 4 weeks (for all other malignancies). Treatment was continued until disease progression, death, or grade 3–4 toxicity. A brain MRI and CT chest/abdomen/pelvis were obtained every 6 weeks for restaging purposes. The primary endpoint was OS3. 9 patients on trial were included for single-cell analysis; 5 of whom had serum and CSF sampling at multiple time points during treatment. 5 patients had metastatic breast cancer, 2 patients had high-grade glioma, 1 patient had ovarian cancer, and 1 patient had esophageal cancer.

Further details of the subjects’ clinical courses, including cytology, steroid status, and Ommaya/VP shunt status are provided in Supplementary Data 1.

CSF cell extraction

Cerebrospinal fluid (CSF) from patients was extracted via an Ommaya reservoir or ventriculoperitoneal shunt (VPS) as part of clinical care. CSF not required for clinical testing was spun at 800 G for 5 min to pellet cells, and resuspended in PBS (ThermoFisher 10010023, Ca/Mg-free). Red blood cells (RBCs) were lysed using ACK lysis buffer (ThermoFisher A1049201) for 4 min on ice to remove RBCs. Cells were then washed with sterile PBS and spun down at 800 G for 5 min, resuspended as a single-cell suspension in RPMI (Corning) with 10% FBS (ThermoFisher 10082-147) for subsequent scRNA-Seq. Cytology was performed whenever possible from available CSF; cytology results for all available samples are given in Supplementary Data 1.

Peripheral blood lymphocyte (PBL) extraction

A 10 cc tube of blood was collected and processed within 3 hours of blood draw. 15 mL of Lymphoprep (STEMCELL Technologies, Catalog #07801) and 10 mL of phosphate-buffered saline was added to the blood. This mixture was then centrifuged at 1200 G for 12 minutes. The supernatant was then poured out and 10 mL of phosphate-buffered saline was added. This mixture was centrifuged a second time at 500 G for 5 min. This supernatant was poured out and 1 mL of CryoStor CS10, Cryopreservation Freeze Media (STEMCELL Technologies, Catalog MSPP-07930) was added to the pellet. This mixture was frozen at −80 °C, then later thawed on ice, then to room temperature, then processed using Seq-Well as described below for CSF-derived cells.

Extraction and sequencing of cell-free DNA

For blood, a 10 cc tube was first centrifuged at 500 g for 10 min. Afterwards, the supernatant was extracted and then centrifuged again at 1000 g for 10 min. The second supernatant was then used for serum cell-free DNA extraction and sequencing. For CSF samples, a 3 cc tube of sample was centrifuged at 400 G for 5 min.

Extraction of cell-free DNA from banked plasma and centrifuged CSF was done using an automated liquid handler at the Broad Institute’s Blood Biopsy Lab. Sequencing was then conducted by the Broad Institute’s core facility.

scRNA-Seq with Seq-Well

Resuspended CSF cells were profiled using the Seq-Well platform for massively parallel, high-throughput scRNA-seq for low-input clinical samples. A complete description of methods is available online⁴³. A complete list of primers described in Gierahn et al.⁴³ is additionally provided in Supplementary Data 10. Briefly, cells from each CSF sample were manually counted (Bal Supply 808CI) and loaded onto one array preloaded with barcoded mRNA capture beads (ChemGenes). All samples retained fewer than 10,000 cells with the exception of two (CSF029-4 & DFCI010-4; ~100,000 cells). Thus, all available cells were loaded onto a single array, except CSF029-4 and DFCI010-4 where ~10,000 cells were loaded. The loaded arrays containing cells and uniquely barcoded oligo-dT beads were then sealed using a polycarbonate membrane with a pore size of 0.01 μm, which allows for the exchange of buffers but retains biological molecules confined within each nanowell. Subsequent buffer exchanges facilitated cell lysis, transcript hybridization, and bead recovery before performing reverse transcription en masse. Following reverse transcription using Maxima H Minus Reverse Transcriptase (ThermoFisher EP0753) and an Exonuclease I treatment (NewEngland Biolabs M0293L) to remove excess primers, PCR amplification was carried out using KAPA HiFi PCR Mastermix (Kapa Biosystems KK2602) with approximately 2,000 beads per 50 μl reaction volume. Libraries were then pooled into one tube (with the exception of CSF014-4, CSF029-4, CSF046-2, CSF104-1, CSF104-3, and CSF123-4, which were pooled to two tubes) and purified using Agencourt AMPure XP beads (Beckman Coulter, A63881) by a 0.6X SPRI followed by a 0.8X SPRI and quantified using Qubit hsDNA Assay (Thermo Fisher Q32854). The quality of each WTA product was assessed using the Agilent hsD5000 Screen Tape System (Agilent Genomics) with an expected peak ranging between 800 and 1500 bp tailing off to beyond 5000 bp, and a small or non-existent primer peak (~100–200 bp).

3′ digital gene expression (DGE) libraries were constructed using the Nextera XT DNA tagmentation method (Illumina FC-131-1096) using index primers as described in Gierahn et al.²⁸. Loaded samples ranged from 600 to 2,000 pg of pooled cDNA, depending on the peak distribution of the WTA product for the sample. Tagmented and amplified sequences were purified at a 0.6× SPRI ratio followed by a 0.9X SPRI yielding library sizes with an average distribution of 400–750 base pairs in length as determined using the Agilent hsD1000 Screen Tape System (Agilent Genomics). Samples DFCI010-4, CSF011-1, CSF011-7, CSF014-1, CSF014-2, and CSF014-4, CSF029-2, CSF029-5, DFCI037-1, CSF046-2, CSF050-4, CSF050-7, and CSF073-4 were sequenced using an Illumina 75 Cycle NextSeq500/550v2 kit (Illumina 20024906) at a final concentration of 2.2–2.8 pM. Samples CSF029-4, CSF043-2, CSF043-3, CSF043-4, CSF050-3, CSF050-12, CSF050-19, DFCI056-3, DFCI058-1, DFCI058-5, DFCI058-7, DFCI061-1, DFCI061-2, DFCI062-2, DFCI062-3, DFCI062-4, DFCI062-5, DFCI062-6, CSF091-3, CSF104-1, CSF104-3, CSF119-1, CSF123-4, CSF127-3, and CSF129-1 were sequenced using an Illumina 100 Cycle NovaSeq6000S kit (Illumina 20027464). The read structure in both cases was paired end with read 1 starting from a custom read 1 primer containing 20 bases with a 12-bp cell barcode and 8-bp unique molecular identifier (UMI) and read 2 containing 50 bases of transcript information.

Alignment & Pre-processing of scRNA-Seq data

Read alignment was performed as in Macosko et al.²¹. In brief, for each Illumina sequencing run, raw sequencing data were converted to demultiplexed FASTQ files using bcl2fastq2 based on Nextera N700 & N500 indices corresponding to individual samples/arrays. Reads were then aligned to hg19 genome using the dropseq_tools v2.1.0 pipeline maintained by the Broad Institute using standard settings. Individual reads were tagged according to the 12-bp barcode sequenced and the 8-bp UMI contained in Read 1 of each fragment. Following alignment, reads were binned onto 12-bp cell barcodes and collapsed by their 8-bp UMI with a hamming distance correction of 1. DGE matrices (genes-by-barcode) for each sample were obtained from quality filtered and mapped reads, with an automatically determined threshold for barcode count.

DGEs from each sample were individually culled and inspected by unsupervised analysis before inclusion into the full analysis by a combination previously described methods^29,30. Each barcode was initially scored on: (1) average expression of a list of curated housekeeping genes (Supplementary Data 5) and (2) complexity, estimated by the total number of genes detected. All sequenced samples were cut to exclude barcodes with low complexity/housekeeping gene expression (400 gene complexity cutoff, housekeeping gene expression cutoff of 1.6 log₂(tp10k)). Each sample was then inspected using unsupervised analysis to further identify and curate potential analyzable single cells. Individual arrays were analyzed to determine the extent of cross-cell type gene expression contamination. Minimal cross-cell type gene expression contamination existed between immune subsets, and select barcodes were filtered out according to cross expression of marker genes from other immune subsets. Further restrictive analyses incorporating lowered complexity thresholds and count-based downsampling was performed to control for technical confounders wherever relevant. Following curation, all samples were combined and genes expressed in at least 1% of the remaining barcodes were retained (in the case of WME calculations used in Fig. 4c we retained genes expressed in at least 0.0875% of remaining barcodes; in the case of WME calculations used in Supplementary Fig. 9, we retained all genes). Consecutive samples from the same patient were combined by assigning zeros to all undetected genes per sample and concatenating columns. miRNA and T cell receptor chain genes were subset and saved before cutting genes to ensure information was not lost. This curated, UMI-collapsed data was then normalized to 10,000 UMIs per barcode (tp10k) and log₂-normalized before being input into Seurat¹⁴ v2.3.4 (https://github.com/satijalab/seurat) for further analysis. This yielded a Seurat object of 34,742 single cells and 8,156 genes, with different genes being used for more specific analysis (such as T cell analysis). The 37 individually sampled time points averaged 890.8 cells per sample with a range between 103 cells and 1,946 cells (Supplementary Data 1).

Alignment & Pre-processing of PBL-derived scRNA-Seq data

Read alignment was performed identically as with CSF-derived scRNA-Seq data. Barcodes were selected from DGEs with a 200 gene complexity cutoff. Unsupervised analysis was performed jointly with CSF-derived cells, with the same complexity cutoff, from patients with PBL-derived data (P010, P043, P046, P073). All genes were retained. In total, 810 PBL-derived immune cells were detected.

For PBL vs CSF comparison, to account for differences in cell qualities between these samples, all cells had their total UMI count adjusted to 500 (the approximate mean UMI count of PBL-derived cells with complexity greater than 200). This was done by randomly selecting 500 UMIs for each cell, with sampling probabilities given by the pre-adjusted UMI count for that gene, in that cell. This data was then normalized to 10,000 UMIs per barcode (tp10k) and log₂-normalized. Differential expression and module score calculations were performed as with CSF-derived scRNA-Seq data.

Unsupervised transcriptomic analysis

Before performing dimensionality reduction, a list of the 2,359 most variable and highly expressed genes was generated by including genes with an average normalized and scaled expression value greater than 0.1 and with a dispersion (variance/mean) greater than 0.1. We then performed principal component analysis (PCA) over the list of variable genes. For both uniform manifold approximation and projection (UMAP) and SNN (shared nearest neighbor) clustering, we used the first 30 principal components. We used FindClusters within Seurat (which utilizes a SNN modularity optimization-based clustering algorithm) with a resolution of 0.7 and UMAP with minimum distance of 0.2 and number of neighbors of 50 to identify 27 clusters across the 37 input samples. 11 of these clusters were collapsed due to gene expression similarity evaluated by prior biological knowledge (7 extraneous divisions in cluster 0, 1 extraneous division in cluster 1, 2 extraneous divisions in cluster 3, and 1 extraneous division in cluster 4) to arrive at 17 total biological clusters across all samples.

Dimensional reduction on data from the CD8 + T cells and myeloid cells alone was similarly performed using PCA followed by UMAP and SNN clustering, all implemented in Seurat. For CD8 + T cells, principal components 1-6 were used with UMAP parameters of minimum distance 0.3 and number of neighbors 20; a resolution of 0.4 was used to identify clusters.

Cell type identification and within cell type analysis

To identify genes that defined each cluster, we performed differential expression using the “bimod” test implemented with the FindMarkers function in Seurat based on a likelihood ratio test designed for single-cell differential expression incorporating both a discrete and continuous component. Thresholds were set at an average log-fold difference 0.2 and adjusted p-value (Bonferroni) less than 0.05. Top marker genes with high specificity were used to classify cell clusters into cell types (Supplementary Data 2,4) based on literature precedent. Two closely related clusters (T/NK clusters) were merged based on biological curation and analysis of hierarchical cluster trees yielding the twelve unique clusters. For T cells, we subclustered first on treatment condition, as we found that the original clusters were dependent on this metadata. The process used for clustering and subset identification was adapted for each iteration of clustering to optimize the parameters of variable genes, principal components, and resolution of clusters desired. Following identification of canonical subsets – CD4 + T cells, CD8 + T cells, and NK cells – these identities were assigned to the main T/NK cluster cells.

One cluster, cluster 15, containing 50 cells was not classified as immune or malignant. All cells in this cluster came from sample CSF011-1 and upregulated genes associated with neuronal expression, and this cluster was classified as “other.”

NK cell clusters within the pre-treatment and post-treatment T/NK groups were annotated based on expression of CD2 and FCRG3A (CD16), lack of expression of CD3 genes (CD3D, CD3E, CD3G).

The plasmacytoid DC (pDC) and conventional DC (cDC) clusters were distinguished from the other innate cells as dendritic cells, and then the differentially expressed genes between the two clusters were enriched using GSEA. The top GSEA hits on both gene lists distinguished cDCs and pDCs (Supplementary Data 9).

Comparisons of abundance of T cells were made across time points with at least 20 T cells detected (34 of 37 time points). Comparisons of proliferation of CD8 + T cells were made across time points with at least 15 CD8 + T cells detected (31 of 37 time points).

Differential expression and scoring over gene sets

To identify differentially expressed genes within cell types and subtypes across treated and untreated conditions, we again used the ‘bimod’ setting in FindMarkers implemented in Seurat. To determine the scores of gene sets and pathways, such as IFN-γ response and antigen processing, we used the ‘AddModuleScore’ function in Seurat to construct a mean score of supplied genes subtracting a background score constructed from a random selection of genes in bins of average expression across all cells. When comparing scores within a specific subset of cells, AddModuleScore was constructed only over that subset, and recalculated if the subset was further partitioned. Tumor cell scores were calculated both across all patients (to compare pre-treatment and post-treatment time points across patients) and within individual patients (to compare across time points within patients). For specific comparisons of AddModuleScore-derived signatures with large differences in complexity between groups of cells, an upper complexity threshold and count-based downsampling were used to examine the possibility of complexity-confounded effects. No such effects were observed in comparing between tumor cells across patient and within patient.

IFN-γ Response, Antigen Presentation, and Exhaustion Signatures

IFN-γ response signature and exhaustion signatures were obtained from GSEA (HALLMARK_INTERFERON_GAMMA_RESPONSE, various signatures from Wherry et al. 2007), provided in Supplementary Data 5. Antigen presentation signature was compiled following a search of the literature and is provided in Supplementary Data 5.

Inferred CNV analysis of Patient 043

To more accurately infer CNV patterns in high-complexity (complexity > 1000) tumor cells with sub-chromosomal resolution, we group genes into fixed length windows of 200 genes consecutive along the genome, removing from consideration those genes in the uppermost decile of dropout rate, as well as all immunoglobulin genes. All possible windows were considered where all included genes reside on the same chromosome. We converted the log (TP10k + 1) gene expression profiles to TP10k ones. We then took the mean TP10k expression over genes in a window, neglecting the highest 5 gene expressions in that window. This vector of values is hereafter referred to as the unnormalized Windowed Mean Expression (uWME).

Having identified the malignant cells for each patient, we additionally computed a normalized version of the uWME as follows: the uWME from all patients’ non-malignant cells were averaged for each window across patients. HLA-* and associated genes on the 6p arm exhibited particularly strong hematopoietic expression; therefore, the means of these windows were imputed with the mean (windows) of the mean (patients) WME for all other windows. These values we refer to as the mean non-malignant uWME.

We normalize uWME for malignant cells by dividing the window uWME by the mean non-malignant uWME for each window, hereafter referred to as Windowed Mean Expression (WME).

To reduce possible confounding factors due to experimental or batch effects during subsequent clustering analysis, we converted the WME values in each single cell to ranks, hereafter referred to as the ranked, normalized WME (rWME). PCA and UMAP were performed on the rWME using the first 50 principal components of all tumor cells. In this P043, the UMAP-obtained clustering was concordant with that achieved via agglomerative clustering. To perform this clustering, we used as a distance metric 1-τ_K, where τ_K is the Kendall’s т coefficient between the WME of all pairs of cells. Agglomerative clustering was performed with a weighted linkage to obtain four clusters; two clusters contained single cells, and two other clusters contained 128 and 62 cells and were denoted descendant and ascendant respectively.

To support the subclonal hypothesis without relying on either inferred CNV profiles or unsupervised clustering thereof, we performed a supervised comparison of single-cell expression profiles at each time point to both the early and late cfDNA-derived WES copy ratios (Fig. 4d). We calculate the Kendall’s т correlation for all genes’ total copy ratio and single-cell expression, for all single cells. Then we calculate the difference in correlation for all single cells when using total copy ratio from time point 4 (late) vs. time point 2 (early). We observe that CSF043-2 single cells exhibit correlations more similar to WES from time point 2, and that CSF043-4 single cells exhibit correlations more similar to WES from time point 4. At CSF043-3, we observe bimodality in the distribution of the difference of Kendall’s т correlations. Additionally, we observe that single cells derived from post-treatment time points (CSF043-3 and CSF043-4) exhibit anti-correlation between their IFN-γ response score (Fig. 4d) and the difference in Kendall’s т correlations between total copy ratios derived from WES at time point 4 vs. time point 2.

We note that the relative populations of the two identified clusters in P043 varied significantly across time (Fig. 4e). We plotted, for each gene, the mean purity corrected tCR vs change in the WME between all possible pairs of time points. The purity corrected tCR has the following form:

where p is purity of sample calculated by ABSOLUTE⁴⁴ and tCR_germline = 1.

This relationship is demonstrated in Supplementary Fig. 11, showing that the windowed expressional change between these clusters is concordant with the change in WES-derived tCR between any two time points. This concordance is robust to considering only the cells obtained at time point 3 (i.e., the correlated changes in single-cell expression and cfDNA-derived CNV profile cannot be attributed to batch effects confounding the observed scRNA-seq profiles).

Windowed-mean expression results were compared to the InferCNV R package from the Broad Institute, and broad amplifications and deletions were concordant between the two approaches (Supplementary Fig. 11a,b).

Statistical analyses

Statistical analyses of differential expression were performed using Seurat v2.3.4 implemented in RStudio. All statistical tests of distributions, cluster diversity, and change in representation, etc. were performed using base R packages implemented in RStudio. All statistical tests of gene set enrichment were performed using piano v1.22.0 and implemented in RStudio for all except enrichments of cluster markers for the full dataset, which was implemented in R. All violin plots and boxplots were generated using ggplot2 without modifications to smoothing or density. Boxplot rectangles encompass the 25^th to 75^th percentile with outliers as individual points above and below the rectangle. Overlapping genes between IFN response and antigen processing signatures were removed from both before utilization.

As scores followed non-normal distributions as tested for using a Lilliefors normality test, we used a Wilcoxon rank-sum test where indicated for determining statistical significance. For scores in single-cell data, we report effect sizes in addition to statistical significance as an additional metric to capture the magnitude of the effect observed. The calculation was performed as Cohen’s d where: effect size d = (Mean₁-Mean₂)/(s.d. pooled). In Supplementary Fig. 5, the calculation of Cohen’s d was modified to d_pair= (Mean₁-Mean₂)/(s.d.₂), where the difference in means is normalized by the standard deviation of the pre-treatment group. All p-values subject to the multiple comparisons problem (such as marker identification by differential expression) were adjusted by Bonferroni correction. Wilcoxon rank-sum tests were calculated via the R command wilcox.test. Related Student’s t-test p-values were computed in python 3.7.7 using the function scipy.stats.ttest_rel from scipy v1.5.4. Theil-Sen slopes were computed in python 3.7.7 using the function scipy.stats.theilslopes from scipy v1.5.4. Kendall-tau correlations and associated p-values were computed using the function scipy.stats.kendalltau from scipy v1.5.4.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.