Breast cancer patient-derived whole-tumor cell culture model for efficient drug profiling and treatment response prediction

Significances There is an urgent demand for discovering more accurate and predictive biomarkers and tools to facilitate precision oncology. Next-generation sequencing (NGS) technologies receive much attention but will probably never decipher the biological behavior of cancer cells completely. However, our established whole-tumor cell culture model, in combination with NGS tools and functional assays, could provide us with a platform to efficiently identify drug sensitivity and resistance for individual patients. We consider it a breakthrough over other existing ex vivo models by considering the tumor-stromal interactions to represent a more unbiased snapshot of the original patient's disease.


Approval and collection of clinical tumor material
Biobank approval was obtained from the Stockholm medical biobank. Before surgery or neoadjuvant treatment, informed consent in accordance with the Declaration of Helsinki was given to patients for signature. Clinical information for each patient specimen involved in this study, as well as the individual treatment details for the validation cohort patients, can be found in List S1.

Flow cytometry
The TSCs and WTCs samples were first washed in PBS, followed by gentle digestion with Dispase at 37°C for 5-10 mins, and passed through a 70 mm nylon strainer (Miltenyi Biotec) for the generation of single-cell suspension. Cells were then re-suspended in PBS and incubated for 30 mins on ice with Human Fc Receptor Binding Inhibitor ThermoFisher Scientific). Cell viability is accessed by LIVE/DEAD Fixable Dead cell stain kit and staining is performed with the antibodies listed in the supplementary file List S2. After staining for 30 mins and subsequent washing steps, samples were analyzed on the LSR II BD flow cytometer (BD Biosciences). Data was then analyzed by FlowJo software (BD Biosciences).

Real-time Imaging
Brightfield and phase-contrast images under the 4X objective were acquired every six hours on the IncuCyte S3 system (Essen BioScience). 1x10 4 cells were seeded and cultured in an ultra-low attachment 96-well plate (174925, Thermo Fisher Scientific) with the culture medium described above. CellEvent™ Caspase-3/7 Green (C10423, Thermo Fisher Scientific) at 1 μM was also added to observe apoptosis. After four days, chemotherapeutics at 100 nM final concentration were added along with media top-up. Spheroid invasion analysis was performed using top hat segmentation with IncuCyte analysis software according to the manufacturer´s protocols.

Whole-genome sequencing analysis
The somatic SNVs and Indels were called using MuTect2 version 3.8.0(1) and Strelka2 version 2.8.2 (2). The matched normal samples were used to distinguish germline variants from somatic ones. Only the somatic SNVs and indels detected by both MuTect2 and Strelka were included in the analysis as a measure against spurious variants. The variant effect prediction was performed using SnpEff version 4.3(3), and the most deleterious variants were identified and assigned on a per gene basis (4).
The copy number analysis was performed using Control-FREEC version 11.5(5). In detecting the (somatic) copy number alterations (CNAs), the window size was set to 1 Kb, and the matched normal samples were provided as input to Control-FREEC along with the TSCs or WTCs. For each gene, the output by Control-FREEC log2 ratios (between the tumor and matched normal samples) at those bins that overlap with the gene body by at least 1 bp were identified and averaged.
The t-SNE plot was generated using the tsne R package version 0.1.3 and was based on the somatic SNVs and CNAs of the BC-specific driver genes and actionable genes (7,8).

RNA sequencing analysis
Alignment to the human genome assembly (build GRCh37) was carried out using Tophat aligner version 2.0.4. Merged bam files were sorted, and PCR duplicates were marked using Picard MarkDuplicates version 1.29 (http://broadinstitute.github.io/picard). Gene counts were calculated with HTSeq count version 0.6.1 (9) with duplicates included. The RNAseq count data were prefiltered by only keeping in the analysis of those genes that had > 10 reads per million mapped reads in at least two libraries. Data were then normalized using the TMM method(10) available in the edgeR package version 3.24.3(11). Differential expression between two conditions was analyzed with the gene-wise negative binomial generalized linear model using the R/Bioconductor edgeR package version 3.24.3(11, 12). Adjustment for differences between patients was made by an additive linear model with patients as the blocking factor. Empirical Bayes methods were used to moderate the degree of dispersions across genes. The tests can be viewed as analogous to (moderated) paired t-tests. Enrichment analysis of the Hallmarks and Canonical pathways gene-set collections in the Molecular Signatures Database (MSigDB; Broad Institute, version 7.1) was performed using the GSEA software version 4.0.3 (Broad Institute) with genes pre-ranked according to scores from the differential expression analysis (13). Multiple testing was controlled by calculating the expected false discovery rate (FDR), according to Benjamini & Hochberg(14). The absolute abundance of eight immune and two stromal cell populations was estimated with Microenvironment Cell Populations-counter (MCP-counter) method(15) using the R package MCPcounter version 1.1.0.
All gene expression analysis was done in R/Bioconductor version 3.5.3 unless otherwise specified.

Cell viability assay analysis
Trastuzumab and pertuzumab were freshly prepared in water according to the manufacturer's instructions. All the other compounds were stored as 10 mM aliquots in DMSO according to the manufacturer's recommendation. For the drug profiling assay, each compound covered five concentrations ranging from 10 μM to 1 nM (2 μM to 0.2 nM for trastuzumab and pertuzumab) in 10-fold dilutions and was dispensed using the acoustic liquid handling system Echo 550 (Labcyte Inc) to make spotted 384-well plates. For the neoadjuvant setting validation assay, we updated the cyclophosphamide into its active metabolite form 4-hydroperoxy cyclophosphamide (4-OOHcyclophosphamide). Each relevant compound covered eight concentrations ranging from 10 μM to 1 nM (2 μM to 0.2 nM for trastuzumab and pertuzumab) and was dispensed using the Tecan D300e Digital Dispenser (Tecan) to make spotted 384-well plates. In both experiment settings, a total volume of 40 nl of each compound condition was dispensed into each well, for limiting the final DMSO concentration to 0.1% during the treatment period.
After 4-5 days of cultivation, WTC spheroids were dissociated into single cells by Dispase digestion at 37°C and filtration with a 70 μm nylon cell strainer, before being suspended in the freshly prepared culture medium again. Afterward, 2x10 3 cells in 40ul/well of medium were transferred into the drug-spotted 384-well OptiPlate (6007290, PerkinElmer) using the MultiDrop Combi dispenser (ThermoFisher Scientific). The plates were cultured in an incubator with an extra humidity supply for 96 hours. The cell viability was assessed using CellTiter-Glo® 3D cell viability assay (G9683, Promega) according to the manufacturer's instruction and reading luminescence by a Tecan spark 10M microplate reader (Tecan). Positive controls (benzethonium chloride 500 μM, 53751-50G, SigmaAldrich) and negative controls (DMSO, D2650-100ML, SigmaAldrich) from the plate were used to normalize the data.
Quantitatively scoring of drug sensitivity was estimated using the drug sensitivity score (DSS), which captures and integrates the multiparametric dose-response relationships into a single metric (16). Specifically, the raw dose-response readout for each of the tested drugs was first subjected to quality assessment by two independent observers, and response values with outlying behavior were removed. Outliers were defined as response values with an irregular response pattern, or response values with a significantly higher cell count than their corresponding negative controls (maximum one point per drug-response readout was removed due to the limited number of concentrations; in the rare cases of two or more outlying points, the specific dose-response readout was excluded from downstream analysis). The raw dose-response readout was then normalized in relation to the negative and positive controls (median values used), and the relative inhibition (%) was calculated. The normalized dose-response data were fitted using a non-weighted four-parameter logistic regression function, in a similar way as in (16,17), and the bottom asymptote of the curve was fixed to zero inhibition for each drug (18). The R function "get.curve.data()" in script "calculation_of_css.Rmd" was used for curve fitting, available from GitHub at (17). The DSS (type 3), which varies between 0 (insensitive) and 100 (highly sensitive), was calculated using the R package DSS version 1.2, as previously described (16). The value of 10% for the drug response's minimum activity level was used.
Statistical analyses to assess the differences in the DSS across patients over known biomarkers were performed using either a Kruskal-Wallis test or a Wilcoxon-Mann-Whitney test. P-values were adjusted using the Benjamini & Hochberg procedure (14) (function "p.adjust" in R package stats version 3.4.0) to account for multiple comparisons. Differences were considered significant if pvalue < 0.05 unless otherwise stated. All analyses were performed in R version 3.4.0 (http://www.rproject.org/) unless otherwise specified.
The cell viability assay and DSS calculation for the validation cohort patient samples were performed by researchers, and the clinical efficacy was independently evaluated by clinicians. The researchers performing the drug profiling were blinded to the interventions and clinical outcomes of the patients during the treatment period.

NanoString nCounter® Breast Cancer 360 Panel analysis
The raw data of the assay was assessed using several quality assurance (QA) metrics to measure imaging quality, oversaturation, and overall signal-to-noise ratio. All samples satisfying QA metric checks were background corrected (background thresholding) using the negative probes and normalized with their mean minus two standard deviations. The background-corrected data were then normalized by calculating the geometric mean of five housekeeper genes, namely ACTB, MRPL19, PSMC4, RPLP0, and SF3A1 (19). Pathway scoring was then performed on the prenormalized data, where pathway-level information from a group of genes was extracted using the first principal component of their expression data (20).

cycles 4 cycles P end
Day 0

Biopsy taken for WTC EC start
Day 60