Patient-derived xenograft models capture genomic heterogeneity in endometrial cancer

Background Endometrial cancer (EC) is a major gynecological cancer with increasing incidence. It comprises four molecular subtypes with differing etiology, prognoses, and responses to chemotherapy. In the future, clinical trials testing new single agents or combination therapies will be targeted to the molecular subtype most likely to respond. As pre-clinical models that faithfully represent the molecular subtypes of EC are urgently needed, we sought to develop and characterize a panel of novel EC patient-derived xenograft (PDX) models. Methods Here, we report whole exome or whole genome sequencing of 11 PDX models and their matched primary tumor. Analysis of multiple PDX lineages and passages was performed to study tumor heterogeneity across lineages and/or passages. Based on recent reports of frequent defects in the homologous recombination (HR) pathway in EC, we assessed mutational signatures and HR deficiency scores and correlated these with in vivo responses to the PARP inhibitor (PARPi) talazoparib in six PDXs representing the copy number high/p53-mutant and mismatch-repair deficient molecular subtypes of EC. Results PDX models were successfully generated from grade 2/3 tumors, including three uterine carcinosarcomas. The models showed similar histomorphology to the primary tumors and represented all four molecular subtypes of EC, including five mismatch-repair deficient models. The different PDX lineages showed a wide range of inter-tumor and intra-tumor heterogeneity. However, for most PDX models, one arm recapitulated the molecular landscape of the primary tumor without major genomic drift. An in vivo response to talazoparib was detected in four copy number high models. Two models (carcinosarcomas) showed a response consistent with stable disease and two models (one copy number high serous EC and another carcinosarcoma) showed significant tumor growth inhibition, albeit one consistent with progressive disease; however, all lacked the HR deficiency genomic signature. Conclusions EC PDX models represent the four molecular subtypes of disease and can capture intra-tumor heterogeneity of the original primary tumor. PDXs of the copy number high molecular subtype showed sensitivity to PARPi; however, deeper and more durable responses will likely require combination of PARPi with other agents. Supplementary Information The online version contains supplementary material available at 10.1186/s13073-021-00990-z.

Large structural rearrangements were identified using an in-house developed tool qSV for WGS (4). Somatic copy number alterations (CNA), tumor ploidy and purity were determined with ascatNgs for WGS (6) and GAP for SNP arrays (1). Homozygous deletions (CN<1), deletions (CN<2) gains (CN≥ploidy + 1) and amplifications (CN≥ploidy +3) were considered in the analysis. Copy number alterations (CNA) and structural variant (SV) events were annotated against Ensembl. Percentage of genome with CNA was determined from SNP array and WGS CNA data for PDX tumor samples, as percentage of genome (autosomal chromosomes) with copy number ≠ 2. Tumor purity for mismatch-repair deficient (MMRd) models was determined from distribution mode of somatic variant allele frequencies multiplied by tumor ploidy of 2, since there was not sufficient copy number changes to accurately estimate tumor purity from CNA calls. Tumor and ploidy purity for WGS data was estimated by ascatNgs.

Heterogeneity analysis
Somatic mutation comparison between primary and PDX tumor samples was performed for four mismatch-repair deficient (MMRd) models with WES data and three carcinosarcoma models with WGS. A union of all pass-filter missense variants was generated for each model, and then the genomic positions of these variants were interrogated with qBasePileup to capture any variants that may have been missed during variant calling, because of low frequency or other quality filters. The overlaps of pileup variants were visualised using eulerR package. The generated pileup list of variants with allele-specific copy number information, extracted from ascatNGS (6) or GAP (1) output (for WGS and SNP arrays, respectively), was then used to identify mutational clusters using PyClone (v0.13.1) (7) with pyclone_beta_binomial emission density, random seed set to 20 and thin parameter set to 10. The top five mutational clusters with the most variants were selected for determining clonal evolution using ClonEvol package (8).
Copy number clonality analysis was performed on the WGS data of tumor-normal pairs using Battenberg (v2.2.5) with default settings.

MSI status
The level of microsatellite instability (MSI) was assessed using MSIsensor (v0.2) on tumornormal pairs of primary and PDX tumor samples using suggested parameters for WES and WGS data (9). Samples with an MSI score of >3 were classified as MSI-high.

HRD score assessment
Homologous recombination deficiency (HRD) scores were assessed on SNP array and WGS data using scarHRD (10) package on the allele-specific copy number information determined by either GAP (1) or ascatNGS (6).

Signature analysis
Mutational signature analysis was performed using two approaches with SigProfiler and deconstructSigs. De novo signatures were identified on somatic single nucleotide variants from WES data using non-negative matrix factorization (NMF) with SigProfiler, as previously described (11). The identified signatures were then compared to the 30 known COSMIC v2 signatures. The optimal number of signatures was chosen based on a number of parameters: stability, reconstruction error and cosine similarities to COSMIC signatures. The contribution of each de novo signature to a sample's mutational profile was assigned using SignatureEstimation package (12). To determine the potential relative contribution of signature 3 (HRD-associated signature), deconstructSigs package (13) was used to estimate the contribution of the sample's mutations to the full catalogue of COSMIC v2 mutational signatures, using default setting. Minimum mutation signature contribution cut-off was set to 15%.
SV signatures were identified using the same approach as used for de novo mutational signature analysis. SV events were classified into 32 previously defined categories based on event type, size and breakpoint clustering (14). SV signature analysis was performed as previously described (15).
HRdetect probability scores were calculated as previously reported (16) using the reported weights and adjustment values (based on SNV mutational signatures, SV signatures, HRD sum scores and proportion of deletions with microhomology). Deletions with microhomology were identified previously developed scripts (17). HRdetect scores >0.7 were used to categorise HRD samples.

Targetable Mutations Analysis
Cancer Genome Interpreter analysis was performed in July 2019 to identify potentially targetable mutations in the PDX models. Short somatic variants, homozygous deletions and amplifications in coding regions observed in all PDX samples were used as input. Only known pathogenic variants and predicted drivers (tier 1, 2) were considered as biomarkers. Biomarkers were matched with drugs using the Cancer Biomarker Database within Cancer Genome Interpreter. The drug prescription output from Cancer Genome Interpreter was filtered to select variants that were "complete" alterations (alterations that match the specific amino acid change in the gene which constitutes an actionable variant for a specific drug). Only drug sensitivity biomarkers with FDA approved drugs or drugs currently in clinical trials were included.
Chemotherapy and steroid compounds were not included. LG EEC -low-grade endometrioid endometrial cancer; DWD -died with disease; AWD -Alive with disease; NED -alive with no evidence of disease.  SNV and short indels UCS RPL22