A compendium of next-generation patient-derived models for diverse cancers

Biospecimen processing and quality control

Before extracting nucleic acids, pathology quality control (QC) was performed on all frozen and FFPE tumours and normal specimens. From frozen tissue, a 30 mg or smaller piece was prepared, or an equivalent amount of scrolls were cut from a FFPE block, with sections stained with haematoxylin and eosin taken from the top and bottom of each specimen. The stained slides were scanned and their pathology reviewed to confirm that the tumour was consistent with the reported histology, and to assess the per cent tumour nuclei, per cent necrosis and other pathological features. The tumour nucleus and necrosis percentages from the top and bottom slides were averaged, and tumour specimens with an average of ≥50% tumour nuclei and ≤20% necrosis were submitted for nucleic acid extraction. Normal tissues were rejected if they had any detectable tumour cells.

DNA was extracted from normal blood and saliva specimens. DNA and RNA were co-extracted from tumours, solid normal tissues and cancer models. Before extraction, cancer models were first washed with ice-cold PBS to remove residual Matrigel. Frozen tissues and cancer models were homogenized using a Qiagen TissueLyser, and RNA and DNA were extracted using a modification of the DNA/RNA AllPrep kit (Qiagen). The homogenate was applied to a Qiagen DNA column, and the flow-through was processed using a mirVana miRNA Isolation kit (Ambion). FFPE specimens were deparaffinized and cells were lysed. The pellet underwent DNA extraction using an AllPrep FFPE kit (Qiagen), whereas the supernatant underwent RNA extraction using a Highpure miRNA kit (Roche). Blood specimens were extracted using a QiaAmp DNA Blood Midi kit (Qiagen), and saliva was extracted using a Gentra Puregene Buccal Cell kit (Qiagen).

DNA was quantified by PicoGreen assay and RNA was quantified by measuring the absorbance at 260 nm with a UV spectrophotometer. DNA quality was assessed by 1% agarose gel electrophoresis to confirm high-molecular-weight fragments. RNA was analysed using a RNA6000 Nano assay (Agilent) on an Agilent Bioanalyzer, which returns an RNA integrity number (RIN) for RNA from frozen tissues, or a DV200 for RNA extracted from FFPE specimens. All extracted DNA was subjected to a custom Sequenom single-nucleotide polymorphism (SNP) panel or an AmpFISTR Identifiler (Applied Biosystems) to verify that all specimens representing a case were derived from the same patient.

Genome sequencing

Data processing

The Genomics Data Commons (GDC) DNA-seq alignment pipeline⁵⁹, was used to map WGS reads to a customized version of GRCh38. In brief, reads were aligned using BWA-MEM (v.0.7.15)⁶⁰, followed by sorting, merging and duplicate-marking using Picard Tools (v2.26.10). Base quality scores were recalibrated using GATK (v.3.7.0)⁵⁸ and BQSR. Sample contamination levels were estimated using GATK ContEst, and samples with contamination levels exceeding 4% were excluded from further analyses.

SNV and indel calling

The New York Genome Center pipeline. The New York Genome Center (NYGC; v.6) somatic SNV–indel-calling pipeline⁶¹ was run for each tumour–normal and model–normal pair, starting from the GDC-aligned BAM files. In brief, SNVs, multi-nucleotide variants (MNVs) and indels were called using MuTect2 (GATK v.4.0.5.1)⁶², Strelka2 (v.2.9.3)⁶³ and Lancet (v.1.0.7)⁶⁴. Indels were also called using SvABA (v.0.2.12)⁶⁵. Candidate indels predicted by Manta (v.1.4.0)⁶⁶ were used as input to Strelka2, as per developer recommendations. Variants were merged across callers and annotated using Ensembl (v.93)⁶⁷, COSMIC (v.86)⁴⁷, 1000Genomes (Phase3)⁶⁸, ClinVar (201706)⁶⁹, PolyPhen (v.2.2.2)⁷⁰, SIFT (v.5.2.2)⁷¹, FATHMM (v.2.1)⁷², gnomAD (r.2.0.1)⁷³ and dbSNP (v.150)⁷⁴, using Variant Effect Predictor (v.93.2)⁷⁵. From the final callset of SNVs and indels, we filtered those that met the following criteria: occurred in two or more individuals in a panel of normal samples⁶¹; those that had a minor allele frequency (MAF) of at least 1% in 1000Genomes or in gnomAD; had a tumour VAF less than 0.0001; had a normal VAF greater than 0.2; had depth less than 2 in either the tumour or the normal sample; or a VAF in the normal sample greater than in the tumour sample. The aforementioned panel of normal samples was constructed from 242 unrelated individuals, of which 148 were sequenced using HiSeqX in the Illumina Polaris project⁷⁶, 73 were sequenced using HiSeqX at the NYGC, 11 were sequenced using NovaSeq at the NYGC, and 10 were sequenced on both HiSeqX and NovaSeq platforms at NYGC.

Broad pipeline (whole-exome sequencing).Whole-exome sequencing (WES) BAM files were processed as described in the section ‘Data processing’. Variant calling was performed as previously described⁷⁷. SNVs were called using MuTect and small indels were called using Strelka. Tumour-in-normal contamination estimation was performed using deTiN, and cross-participant contamination was estimated using ContEst. WES variant calling was performed as previously described⁷⁷. A similar pipeline was used for processing variant calls from WGS.

WashU pipeline. The somatic variant pipeline TinDaisy2 (v.2.6.2; https://github.com/ding-lab/TinDaisy) was built around the common workflow language and rooted in the methods of an earlier GenomeVIP pipeline architecture⁷⁸. For both WGS and WES data, we used TinDaisy2 for calling both SNV and indel variants. In brief, aligned BAMs for tumour and normal samples were processed using VarScan (v.2.3.8)⁷⁹, Strelka2 (v.2.9.10)⁶³, Pindel⁸⁰ and MuTect (v.mutect-1.1.7)⁸¹. Filtering was performed to retain calls with length < 100, normal VAF ≤ 0.02, tumour VAF ≥ 0.05, tumour read depth > 14 and normal read depth > 8. Only variants identified by at least two or more callers were retained. Following annotation using VEP 99 (ref. ⁸²), calls with a population allele frequency of <0.005 were retained, and those that are in dbSnP but not in COSMIC⁴⁷ or ClinVar⁸³ were excluded. Adjacent SNP variants were merged into double nucleotide polymorphisms and higher-order variants, and the resulting VCF output was normalized with bcftools (v.1.10.2)⁸⁴.

Consensus method. The workflow MergeParticipantVcfs was run for all participants in the cohort. First, single-pair VCFs and MAFs from the NYGC (v.6) somatic pipeline, the Broad somatic pipeline and the WashU somatic pipeline were prepared for merging with the MergeVcf workflow. Multi-allelic calls were split, MNVs were labelled and then split into SNVs, centre names were prepended to all annotations and indels were left-aligned and normalized. Prepared VCFs were merged using BCFTools⁸⁴. During the merge–normalization of variants, the merging of SNPs into indels was not allowed. We calculated allele depth, read depth and allele frequency from BAM pileups using a custom method as described previously desribed⁶¹. MNVs were re-established. If a pipeline reported only a subset of the SNVs that constitute an MNV, then the SNVs were also reported as a possible variant.

The merged, pair-level VCFs were merged again in the same fashion as described above using the MergeParticipantVcfs workflow. Next, we ran the MergeParticipantVcfs workflow, MuTect2 (GATK v.4.0.5.1)⁶² on the multisample VCF, in forcecalling mode. A MuTect2 filter with read orientation metrics was run on the force-called results. The filtered MuTect2 calls were used to annotate the merged centre VCF. The final genotype was taken from the MuTect2 forcecalls and all other MuTect2 annotations were included with the prefix “Mutect2Multi_.”

The non-normal aliquots that the variant was called in, and the list of callers supporting the variant, were annotated in the VCF. Variants were also annotated using Ensembl variant effect predictor (v.97)⁶⁷ as well as the databases COSMIC (v.98)⁴⁷, 1000Genomes (Phase3)⁶⁸, ClinVar (09012020)⁶⁹, Polyphen2 (v.2.2.2)⁷⁰, SIFT (v.5.2.2)⁷¹, FATHMM (v.2.1)⁷², gnomAD (Genomes (v.3.1.2) and Exomes (v.2.1.1))⁷³ and dbSNP (v.150)⁷⁴.

The annotated calls were used as input for the MakePairHighConfidenceVcfs workflow and filtered into final HighConfidence VCFs. Calls were removed if MuTect2 failed the call. Calls were added for a pair if they had read support in all non-normal WGS BAMs for that participant (even if the variant was not formally called for this pair).

Copy number calling

Ascat method. Purity and ploidy estimates were generated for each tumour–normal and model–normal pair using AscatNGS (v.4.2.1)⁸⁵ using default parameters.

ABSOLUTE method. Copy number segmentation was performed using fragcounter and ReCapSeg, and subsequently with AllelicCapSeg. AllelicCapSeg output files were used as inputs for ABSOLUTE, with WES (when available) or WGS mutation calls used to infer allelic integer copy number profiles. All ABSOLUTE calls were manually curated to identify the ABSOLUTE with integer copy number states best aligning with the observed copy number values. In samples with no CNAs, the optimal ABSOLUTE solution was inferred from SNV multiplicity.

ReMixT method. We applied ReMixT⁸⁶ to predict allele-specific and clone-specific copy numbers from WGS samples according to previously described methods⁸⁷.

PURPLE method. For each tumour–normal or model–normal sample pair, we generated a B-allele frequency of heterozygous SNP sites using AMBER (v.3.5)⁸⁸, and we determined read depth ratios using COBALT (v.1.11)⁸⁸. We integrated this information, together with somatic SNVs and structural variants (SVs) to estimate the purity, ploidy and copy number profile of the tumour using PURPLE (v.2.54)⁸⁸. Tumour and model samples that did not have a matched normal sample were processed in a similar manner, but with several version changes: AMBER (v.3.9), COBALT (v.1.13) and PURPLE (v.3.4). AMBER, COBALT and PURPLE were developed by the Hartwig Medical Foundation and are freely available on GitHub (https://github.com/hartwigmedical/hmftools).

HATCHet method. For each tumour–normal or model–normal sample pair, we ran HATCHet-2 (ref. ⁸⁹) and estimated purity and ploidy. We used default parameters, but with the following exceptions. The minimum and maximum number of states and state transition probability in the hidden Markov model were set to 20, 40 and 10⁻¹², respectively, whereas ‘diploidbaf’ and ‘maxneutralshift’ parameters were both set to 0.06. The search space for the number of clones was between two and four. A Gurobi commercial optimizer was used as the preferred engine for solving integer linear programming.

Consensus method. For each sample, the ABSOLUTE solution for which values were nearest to the consensus purity and ploidy values was selected as the consensus copy number solution. In cases when the optimal ABSOLUTE purity and ploidy solution and the consensus purity and ploidy solutions were divergent, the consensus solution was still selected. ABSOLUTE forcecalling was performed using each selected purity and ploidy solution with the output being the segmented allelic copy number state at each genomic locus.

SV calling

NYGC pipeline. The NYGC (v.6) somatic SV-calling pipeline⁶¹ was run for each tumour–normal and model–normal pair, starting from the GDC-aligned BAM files. In brief, SV breakpoints were called using SvABA (v.0.2.12)⁶⁵, Manta (v.1.4.0)⁶⁶ and Lumpy (v.0.2.13)⁹⁰. We excluded SVs below 500 bp, merged the rest across callers using bedtools (v.2.26.0)⁹¹, pair-to-pair with the slop parameter set to 300 bp, and requiring the same strand orientation and at least 50% reciprocal overlap. SVs were annotated using 1000Genomes (Phase3)⁶⁸, DGV⁹², gnomAD-SV⁹³ and the SV panel of normals (built from the same individuals as the NYGC SNV–indel pipeline) with bedtools pair-to-pair, using the same parameters as for the merge across callers. We removed from the final callset any SVs that overlapped variants in 1000Genomes, DGV, gnomAD-SV or the panel of normals.

Broad pipeline. SVs were called using Manta⁶⁶, SvABA⁶⁵ and dRanger⁹⁴. Breakpointer⁹⁴ was used to refine breakpoint locations of called variants. The same SV had to have been called across at least two of Manta, SvABA and dRanger (using a clustering window of 350 bp to match the SV to the same event) to show up in the final callset. Certain filters were applied across various steps of the pipeline based on span, mapping quality and reads of support for SV events.

WashU pipeline. We used Manta (v.1.6.0)⁶⁶ for calling SVs on matched tumour–normal WGS data. We then filtered results to retain variants that met the following criteria: (1) the sample site depth was less than 3× the median chromosome depth near one or both variant breakends; (2) the somatic score was greater than 30; and (3) for a small variant (<1,000 bases) in the normal sample, the fraction of reads with MAPQ0 around the breakend did not exceed 0.4. We then converted VCF files to BEDPE format with svtools (v.0.5.1)⁹⁵ (https://github.com/hall-lab/svtools).

MSKCC pipeline. We identified SVs using deStruct (v.0.4.18)⁹⁶ and LUMPY (v.0.2.12)⁹⁰, retaining only the breakpoints called by both methods. We applied filtering based on the following criteria: inter-breakpoint distances of ≤30 bp; deletions smaller than 1,000 bp; breakpoints with fewer than 5 supporting reads in the tumour sample; or any read support in the matched normal sample.

EMBL-EBI pipeline. For tumour–normal and model–normal sample pairs, somatic SVs were called using GRIDSS2 (v.2.12.0)⁹⁷, annotated with RepeatMasker (v.4.1.2)⁹⁸, and kraken2 (v.2.1.2)⁹⁹ and filtered with GRIPSS (v.1.9)¹⁰⁰. The final somatic SV set was further refined and annotated with the copy number profile estimated using PURPLE (v.2.54)⁸⁸. SVs and copy number profiles were visualized using the ReConPlot R package (v.1.0)¹⁰¹. Tumour and model samples without a matched normal sample were processed in a similar manner but with several version changes: GRIPSS (v.2.0.1) and PURPLE (v.3.4).

Consensus method. Consensus SVs from the NYGC Broad, WashU, MSKCC and EMBL-EBI were identified as previously described¹⁰² using bedtools pair2pair with minor modifications. The original consensus algorithm used single-caller VCFs as inputs, whereas the modified pipeline uses multicaller consensus inputs from each centre. Moreover, the original pipeline used a slop value of 400, whereas the modified pipeline uses a slop value of 50. SVs called by pipelines from two or more analysis centres were accepted into the final SV consensus.

Purity and ploidy

ESTIMATE. Tumour purity was computed according to the ESTIMATE algorithm¹⁰³, which generates a composite ESTIMATE score that represents the quantification of predefined gene signatures of non-tumour components (stromal and immune) in gene expression data. Tumour purity was then inferred using the following formula: tumour purity = cos (0.6049872018 + 0.0001467884 × ESTIMATE score).

Consensus purity and ploidy. We computed consensus tumour and model purity and ploidy from the purity and ploidy values from five DNA-based copy number callers (Ascat, ABSOLUTE, ReMixT, PURPLE and Hatchet) and purity values from a sixth RNA-based caller (ESTIMATE). Consensus purity was computed using an iterative process wherein the mean purity across all callers was first computed. If the individual purity values across all six purity callers fell within ±0.2 of the mean purity, then the mean purity was accepted as the consensus purity. If one or more callers fell outside the window, then the furthest caller from the mean was removed and the mean ploidy was recomputed. This process was repeated until all remaining callers converged on a consensus purity value or when fewer than three callers remained. In the latter case, the purity value was accepted as the value from ABSOLUTE or, in cases when there was no ABSOLUTE value, from PURPLE.

Ploidy values for each individual caller were first binned into ploidy classes: haploid (ploidy < 1.5), diploid (1.5 ≤ ploidy < 2.5), triploid (2.5 ≤ ploidy < 3.5) or tetraploid (≥3.5). If values from three or more callers fell in the same ploidy class, then the median ploidy from all callers in that ploidy class was taken as the consensus ploidy. If there was no agreement between three or more callers, then the ploidy value from ABSOLUTE was accepted as the consensus ploidy value or, in cases when there was no ABSOLUTE value, the value from PURPLE was accepted as the consensus ploidy value. All consensus purity and ploidy values went through subsequent manual curation and the ABSOLUTE solution was selected in cases when the consensus purity and ploidy values poorly reflected the underlying copy number states.

Analysis methods

Tumour–model SNV concordance. SNVs included in concordance analysis were filtered for variants with VAF > 0.15 or variants for which all non-normal aliquots contained read support.

Tumour–model LOH concordance. The LOH status for all tumour–model pairs was determined by examining whether the rounded minor copy number of each genomic bin was equal to zero. LOH status was then compared between the tumour and its corresponding model, assigning a value of 1 for concordance and 0 for discordance. The mean LOH concordance was subsequently calculated across all equal-sized bins, which provided a single LOH concordance value for each tumour–model pair.

WGD inference. We determined WGD status by assessing genome-wide CNAs. First, we calculated the fraction of the genome that was affected by copy number gains, weighting each segment by its genomic length. We evaluated the extent of duplication based on thresholds for copy number amplification across the genome. Samples in which a majority of the genome exhibited copy number gains beyond a defined threshold were classified as having undergone a single WGD event, whereas those surpassing a higher threshold were assigned multiple WGD events.

Mutation signature. SNV mutational signatures were computed using SignatureAnalyzer (GPU version)^104,105 with the COSMIC single-base substitution signature reference comprising 96 mutational signatures. During mutational signature factorization, a set of 51 non-redundant mutational signatures were selected. For each tumour and matched model, we generated per-sample mutational signature profiles represented as exposure vectors across these 51 signatures. Concordance between tumours and their matched models (n = 440 pairs) was quantified in terms of the cosine similarity between their respective signature profiles. As a negative control, we generated random tumour–model pairings of equal number, both in the same tumour type and across different tumour types, and calculated cosine similarities for the randomized pairs.

Driver oncogene selection and classification. To assess the differences in the prevalence of driver mutations between tumour and model samples, we first curated a set of driver genes based on previous TCGA studies¹⁰⁶ (Supplementary Table 10). To determine the prevalence of SNVs, small indels and copy number variations (CNVs), we integrated the results from both SNV–indel and copy number analyses. Gene annotations for SNVs and indels were obtained through the discovery pipeline. Only mutations with a predicted impact of moderate or high, as defined in the Ensembl calculated gene consequences table, were included in downstream analyses. An additional mutation category included upstream promoter missense mutations in the TERT gene. Gene-level CNVs were annotated by calculating the mean copy number across the genomic span of each gene. Only genes with high-level amplifications or homozygous deletions were included in the downstream analysis. A gene was classified as highly amplified if it met the following criteria: (1) sample ploidy ≥ 1.5; (2) either the gene copy number was ≥3 and ≥2 times the ploidy (adjusted by –0.5 if the sample was a model); or (3) the gene copy number was ≥7. Conversely, a gene was classified as having a homozygous deletion if its copy number was below 0.3.

For the final visualization in Fig. 2, mutations were classified as either gain-of-function (GOF) or loss-of-function (LOF). For each gene, only mutations that were found in both tumour and matched model samples (that is, shared in a tumour-model pair) were considered. Only mutations shared between tumour and matched model samples (that is, present in both) were considered for classification. If no shared mutations were identified for a given gene, it was labelled as unclassified. Mutation types were grouped as follows: LOF included homozygous deletions, truncating variants (frameshift indels, nonsense and splice-site mutations); GOF included missense mutations, high-level amplifications, TERT promoter mutations and in-frame insertions or deletions. A gene was labelled LOF if more than 15% of its shared mutations were LOF events, otherwise, it was labelled as GOF.

Power calculation. We calculated a measure of the power to detect model SNVs in a given tumour by considering the probability that model SNVs would be detected in the tumour, assuming that there were clonal and at one allelic copy. The expected VAF for a single-allelic clonal tumour mutation at a genomic position with tumour copy number T is given by

Here N is the copy number at the given position of the contaminating cells, which was assumed to be two for autosomes. Then, assuming a binomial distribution of read counts and one read necessary to detect, the average power to detect a clonal SNV was calculated by

where Binomial_CDF is the binomial cumulative density function. Here C_i is the read coverage in the tumour at the position of the ith model SNV. This value is calculated by summing over all N SNV positions identified as having at least one alternative read in the model and any number of reads in the tumour. Only autosomal SNVs were considered for the power calculation.

ecDNA. Raw coverage was calculated from the GDC-aligned BAM files using fragCounter (v.1.0; https://github.com/mskilab-org/fragCounter). These values were corrected using dryClean (v.1.0), a robust PCA-based method that separates the foreground from the background signal, which reduces noise and artefacts. It was run using a panel of 390 normal samples, which we built by selecting random normal samples across different datasets: ICGC DCC ESAD-UK¹⁰⁷, TCGA¹⁰⁸, the MSKCC–WCM–NYGC HRD project¹⁰⁹, The Cancer Alliance at NYGC¹¹⁰, the Hartwig Medical Foundation (https://www.hartwigmedicalfoundation.nl/en/data/), ICGC PanCancer Analysis of Whole Genomes⁷, CCLE² and other publications^{111,112,113,114,115,116}. These corrected values were used as inputs for CBS¹¹⁷ to calculate the tumour–normal coverage ratio and to segment the genome into regions of similar copy number, thereby identifying potential amplifications and deletions. This information, together with the consensus purity and ploidy values and the consensus SVs, were used to construct junction-balanced genome graphs that had high-fidelity copy number profiles using JaBbA (v.1.1)^110,118.

GDC-aligned BAM files and JaBbA-derived copy number profiles were then used as input to the AmpliconSuite-pipeline (v.0.5.2)¹¹⁹. This pipeline is a wrapper for the AmpliconArchitect (v.1.3.r5)¹²⁰ and downstream AmpliconClassifier (v.0.5.3)¹¹⁹ tools. As per developer recommendations, samples derived from the same patient (for example, paired tumour–models) were run as a group using GroupedAnalysisAmpSuite.py, using default parameters.

We then compared the resulting amplicon calls using the comparison tool feature_similarity.py, with default parameters. ecDNAs were considered concordant between tumour and model if the overlapping amplicons were classified the same way and had a Jaccard genomic interval similarity of ≧0.75. A Jaccard genomic interval similarity was calculated as the total length of the intersection of the amplicon footprints, divided by length of the union. We retained amplicon calls if the indicated filter was ‘none’. ecDNA calls, filtered by AmpliconArchitect, were ‘rescued’ back into the callset if they had a passing concordant ecDNA in the associated tumour or model.

Putative cyclic ecDNA reconstructions were generated using the AmpliconSuite-pipeline module Candidate AMplicon Path EnumeratoR. In brief, the tool searches each amplicon graph for the longest cyclic and non-cyclic paths, choosing paths that best explain the observed copy numbers, and filtering based on how well the best reconstruction fits the data. Reconstructions that passed were plotted using CycleViz (https://github.com/AmpliconSuite/CycleViz).

Clonal phylogenies. Force-called SNVs were input into pyclone and phyclone to generate clonal phylogenies. We excluded indels from the analysis. We annotated each SNV with the major and minor allele copy number of the encompassing segment from consensus copy number calling. SNV copy number and supporting read counts were input to pyclone-vi (v.0.1.6)¹²¹, which we ran with a beta-binomial observation model, 10 restarts and a maximum of 40 clusters. We removed clusters comprising less than 1% of all SNVs and clusters that were approximately 0.5 cancer cell fraction across all samples (cancer cell fraction range of 0.3–0.7). The resulting pyclone clusters were input to phyclone (v.0.5.1)¹²², which we ran with a beta-binomial observation model, 16 chains, 100,000 iterations and outlier probability set to 0.1.

Genetic ancestry estimation. Ancestry proportion was determined using ADMIXTURE (v.1.3.0)^123,124, which used a maximum likelihood-based method to estimate the proportion of reference-population ancestries in a sample. To do this, we genotyped reference markers that we generated from 1,964 unrelated 1000Genomes project samples directly on the whole-genome samples using GATK pileup (v.3.4.0). We excluded individuals from the populations MXL (Mexican ancestry from Los Angeles, United States), ACB (African Caribbean in Barbados) and ASW (African ancestry in the Southwest United States) from the reference owing to their being putatively admixed. We further filtered the reference by using only SNP markers with a minimum MAF of 0.01 overall and 0.05 in at least one 1000Genomes continental population. Variants were also pruned on the basis of linkage disequilibrium using PLINK (v.1.9) with a window size of 500 kb, a step size of 250 kb and an r2 threshold of 0.2. The analysis resulted in a proportional breakdown of each sample into five continental populations (AFR, AMR, EAS, EUR and SAS) and 23 populations. We then categorized patients by the continental population of highest proportion.

DNA methylation and epigenetic fidelity

DNA methylation data

DNA methylation was evaluated using the Illumina HumanMethylationEPIC (EPICv1) array (Illumina). We downloaded raw IDAT files produced by the Illumina iScan system from the GDC data portal (https://portal.gdc.cancer.gov). We calculated DNA methylation levels (β values) from the IDAT files using the openSesame pipeline with the default arguments implemented in the R package SeSAMe (v.1.18.4)¹²⁵.

Normal tissue methylation and probe selection

We used normal tissue DNA methylation data from external resources to investigate cancer-associated DNA methylation profiles. We had previously identified 146,385 CpGs that were unmethylated in eight normal tissue types (breast, adrenal gland, liver, lung, ovary, skin, blood and brain) on the EPICv1 array¹²⁶. We processed additional ENCODE normal tissue DNA methylation data from 23 normal tissue samples from the gastrointestinal tract, including oesophagus (muscularis mucosa n = 4, squamous epithelium n = 4), stomach (n = 3), colon (transverse n = 4, sigmoid n = 4) and pancreas (n = 4). We downloaded the EPICv1 IDAT files from the ENCODE data portal¹²⁷ and generated β values using the openSesame pipeline, as described above. We identified 159,361 CpGs that had a mean β value of <0.2 in any of the four gastrointestinal tissue types. Collectively, we selected 144,571 CpGs unmethylated in normal tissues from 12 tissue types to investigate cancer-associated DNA hypermethylation profiles (Supplementary Table 11).

TCGA–TARGET DNA methylation data

We analysed TCGA Pan-Cancer Atlas (PanCanAtlas) DNA methylation data profiled using the Infinium HumanMethylation450 (HM450) array. The raw IDAT files were obtained and processed using the R package SeSAMe. IDAT files from TARGET’s neuroblastoma and Wilms tumour projects^128,129 were downloaded from the GDC data portal and processed using the R package SeSAMe. Our analysis included 130 Wilms tumours and 91 neuroblastomas originating in the adrenal glands.

Joint HCMI and TCGA–TARGET methylation data

We merged the HCMI and TCGA–TARGET DNA methylation data profiled using EPICv1 and HM450 arrays, respectively, to generate a dataset with the probes shared between the two platforms¹²⁵. We excluded samples that had a CpG probe success rate of less than 90%. We also excluded probes with ‘NA’-masked data points that were present in more than 10% of the samples and probes on the X and Y chromosomes. To investigate the cancer-associated DNA hypermethylation profiles, we analysed the probe set that lacked tissue-specific DNA methylation (selected as described above) and then extracted 53,204 CpG sites that acquired methylation (β value of >0.3) in at least two samples in any HCMI cancer type (Supplementary Table 11).

UMAP of DNA methylation data

We performed dimension reduction using NMF on the combined HCMI and TCGA–TARGET DNA methylation data matrix described above. β Values of 0 were replaced with 1.0 × 10⁻¹², and missing values were replaced with zeros. The resulting matrix was subjected to NMF, masking zeros, as implemented in the RcppML R package (v.0.5.6)¹³⁰. We assessed the optimal rank for an NMF model by performing matrix decomposition across ranks ranging from 2 to 200, each with three random initializations, using the crossValidate function in the RcppML R package. We selected an NMF model of rank 170, as this model showed the minimum mean squared error of reconstruction consistently across the three runs. UMAP visualization of the rank-170 NMF model was generated using the umap function with the cosine distance metric implemented in the R package umap (v.0.2.10.0) (Supplementary Table 11). We produced the UMAP in Fig. 4c using the subset of the rank-170 NMF matrix, which included the cancer types represented in both the HCMI and the TCGA–TARGET projects.

Heatmap of DNA methylation profiles

For the four cancer cohorts on which we focused, and from the HCMI–TCGA merged DNA methylation data described above, we identified the top 10% of the most variably methylated CpGs separately in each cancer type (Fig. 4e). We selected the 8,614 CpGs that represented the union of the 4 variably methylated CpG sets. To minimize the influence of variable tumour purity levels on clustering results, we dichotomized the data, using a β value of ≥0.3 to define positive DNA methylation and <0.3 to define a lack of methylation. The sample distance matrix was computed using the Jaccard index, and then unsupervised hierarchical clustering was performed for each TMP cancer subtype. We generated the heatmap using the ComplexHeatmap R package (v.2.20.0)¹³¹.

Tumour–model hypermethylation similarity scores

We assessed DNA methylation-based similarity between a model and its parent tumour based on cancer-associated CpG hypermethylation profiles. We identified the 144,571 CpGs unmethylated in normal tissues as described above. For each HCMI cancer type, we further selected CpGs hypermethylated (β value of >0.3) in at least two samples. Then we compared the Pearson’s correlation coefficient (r) as a similarity score for the following samples: (1) models and their matched tumours; (2) models and unmatched tumours from the same cancer type; and (3) models and unmatched tumours from different cancer types. A P value of paired model–tumour relatedness was then determined by calculating the probability of observing an unpaired model–tumour distance equal to or smaller than the same model’s distance to all unrelated tumours that were not from the same tissue type. After FDR adjustment, model–tumour pairs with adjusted P ≥ 0.1 were considered to be models that were non-concordant with their original tumours.

Transcriptional fidelity and relatedness

RNA analyte processing and sequencing

Quality assurance and QC of RNA analytes. RNA from fresh-frozen models was sent for characterization, whereas for most tumours, samples of FFPE-derived RNA were sent. Fresh-frozen RNA analytes were assayed for RNA integrity, concentration and fragment size. Samples for total RNA-seq were quantified on a TapeStation system (Agilent) and RIN scores were calculated. Model-derived RNAs with RINs > 8.0 were considered high quality. For FFPE samples, we used DV200 and fragment size to evaluate sample quality. Although we targeted input concentrations greater than 100 ng µl^–1, some FFPE samples were lower.

Total RNA-seq library construction. To construct total RNA-seq libraries, we used Illumina Stranded Total RNA Prep with RiboZero Gold, and we barcoded samples with individual tags, following the manufacturer’s instructions (Illumina). We prepared libraries, which we then pooled using an automated liquid-handling system to minimize variance. Typically, these were pools of 38–92 samples, depending on the available capacity on a sequencer. At every step we performed QC. We used a TapeStation system to quantify library concentrations, fragment size and distribution. As needed, pool balance and library quality were assessed using miSeq Nano single-end 50-bp sequencing.

Total RNA-seq. We prepared indexed libraries and ran them on an Illumina NovaSeq 6000, using paired-end 100 bp reads and generating a minimum of 150 million reads per sample library, with a target of greater than 90% mapped reads. In all but a few cases, all data were from the same sequencing run. For the few samples that needed additional read depth, we provided this with a secondary sequencing run. We demultiplexed raw Illumina sequence data and converted these to fastq files while quantifying adapter and low-quality sequences. Samples were assessed for information quality by mapping reads to the human hg38 genome reference, estimating the total number of reads that mapped, the fraction of RNA reads that mapped to coding regions, the amount of rRNA in a sample, the number of genes expressed and the relative expression of housekeeping genes. The samples that passed this quality assurance and QC step were then clustered with other expression data from similar and distinct tumour types to confirm expected expression patterns. We SNP-typed atypical samples to confirm the source analyte. FASTQ files of all reads were then uploaded to the GDC repository and distributed to the analysis teams.

MicroRNA (miRNA)-seq library construction. miRNA-seq library construction used a v4 NEXTflex Small RNA-seq kit (PerkinElmer), then samples were barcoded with individual tags following the manufacturer’s instructions. We prepared libraries on a Sciclone liquid-handling workstation. We performed QC at every step and quantified the libraries using a TapeStation system and an Agilent Bioanalyzer using a Small RNA Analysis kit. Pooled libraries were then size-selected according to specifications from a NEXTflex kit. Post-sequencing quality assurance and QC evaluated the abundance and diversity of miRNA. Data that passed were provided to the GDC.

Transcriptional relatedness measurements

MOMA subtype identification. To compare tumour subtypes in HCMI to those previously identified in TCGA, using the network-based MOMA algorithm³³, we used the OncoMatch algorithm³⁸ (see below). Specifically, we assessed the similarity of each HCMI sample to those in the subtypes identified by MOMA in TCGA. As these methods compare samples on the basis of the conservation of their master regulator proteins, which are highly enriched in mechanistic determinants of tumour cell state^{37,132,133,134}, this approach effectively complements the Celligner and DNA-methylation analyses. Specifically, this analysis helps refine subtype classification by mitigating potential confounding effects in gene expression profiles—such as those related to tissue histology, unrelated to tumour biology—as well as effectively mitigating technical batch effects³³. Master regulators, although rarely mutated, have crucial roles in cancer progression, as they orchestrate transcriptional networks that are disrupted by upstream genomic alterations, which makes them critical therapeutic targets^38,39, including in clinical trials^135,136. As such, they are more conserved in each tumour subtype than the corresponding transcriptional profiles^33,132. In the following sections, we discuss the various algorithms used in this analysis.

Original MOMA analysis. In brief, MOMA is based on the assumption that transcriptional cell states are implemented and homeostatically maintained by small, autoregulated modules of master regulator proteins, comprising transcription factors (TFs) and co-factors (co-TFs). We assessed the activities of all TFs and co-TFs using the VIPER algorithm¹³⁷, which is based on the expression of their transcriptional targets, which we inferred using the ARACNe algorithm¹³⁸. Tumour subtypes (n = 112) were then identified by clustering samples based on TF–co-TF activity, using the ‘partitioning around medoids’ algorithm¹³⁹. Below, we describe VIPER and ARACNe.

MOMA subtype comparison. To compare HCMI samples to those in the 112 MOMA subtypes in TCGA, we first used the metaVIPER algorithm¹⁴⁰ to compute the activity of all TFs and co-TFs. MetaVIPER—a multinetwork version of the original VIPER algorithm¹³⁷ that integrates the protein activities assessed by each network—enabled the use of networks from multiple TCGA cohorts that were matched to the histology of HCMI samples (Extended Data Table 1). Enrichment of the top 50 most differentially activated genes in a HCMI sample (that is, 25 most active and 25 most inactive) in TFs–co-TFs that were differentially expressed in each TCGA MOMA subtype were used to assess their similarity (with the OncoMatch algorithm). Indeed, we have shown that across all TCGA cohorts, >80% of the functional mutations in each sample are in pathways upstream of the top 50 master regulators. Moreover, we have shown that changing the number of master regulators between 20 and 200 does not significantly affect the OncoMatch statistics³³.

MetaVIPER analysis. For the analyses in this paper, we used the most recent version of the metaVIPER algorithm¹⁴⁰. The main difference is that we replaced the original method for performing gene set analysis (aREA) by nonparametric analytical-rank-based enrichment analysis (NaRnEA)¹⁴¹, which improved the assessment of significance for differentially active proteins. For implementing NaRnEA, we used the matrix_narnea() function in the PISCES R package¹⁴². For each sample, TF–co-TF transcriptional targets (that is, regulatory networks) were inferred by ARACNe analysis of one or more lineage-matched TCGA cohorts. In brief, metaVIPER estimates the differential activity of each regulatory protein by integrating the NES statistics produced by VIPER analysis using each of the selected networks. VIPER assesses the activity of a protein by assessing the NES of its activated and repressed targets in genes that were differentially expressed in a signature of interest.

Differential expression signature generation. To remove batch effects between TCGA and HCMI samples, TPM-normalized gene expression profiles were first processed using a variational autoencoder (VAE) model called scGEN¹⁴³. For subtype stratification, for which the goal is to assess differentially active proteins in samples in the same cohort, differential expression signatures for VIPER analysis were optimally computed by comparing each sample to the centroid of the entire cohort. We accomplished this by subtracting the median expression across all samples and then dividing by the median absolute deviation (MAD). To prevent division by zero, MAD values <0.01 were set to 0.01. For both HCMI and TCGA, we downloaded the TPM-normalized gene expression for protein-coding transcripts from the GDC portal¹⁴⁴. The cancer-specific networks used in the analysis were inferred by ARACNe as described below.

Regulatory network inference. Regulatory networks for each HCMI cohort were generated by analysing the gene expression profiles in their lineage-matched TCGA cohorts using ARACNe3 (ref. ¹⁴¹). ARACNe3 is the latest incarnation of the ARACNe algorithm¹³⁸. The algorithm identifies regulatory protein–target interactions on the basis of the greatest conservation of the transferred information, as assessed by computing the mutual information on the direct path and on every indirect path traversing an intermediary TF–co-TF protein, based on the data-processing inequality¹⁴⁵.

To generate cancer-type-specific networks, we applied ARACNe3 to TPM-normalized expression profiles across 32 TCGA cohorts (ACC, BLCA, BRCA, CESC, CHOL, COAD, DLBC, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LGG, LIHC, LUAD, LUSC, MESO, OV, PAAD, PCPG, PRAD, READ, SARC, SKCM, STAD, TGCT, THCA, THYM, UCEC, UCS and UVM). Each ARACNe3 network was inferred by subsampling the gene expression profiles until ≥50 targets were identified for each TF–co-TF in the consensus network. The analysis included 1,645 TFs and 1,556 co-TFs, retrieved from ref. ¹⁴⁶ and ref. ¹⁴⁷, respectively. Networks were pruned to the 100 most significant targets for each regulatory protein, based on consensus mutual information statistics. For TCGA–ESCA, we inferred networks for adenocarcinoma and squamous cell carcinoma, separately.

OncoMatch analysis. OncoMatch³⁸ was used to assess the similarity between HCMI and TCGA samples based on a weighted enrichment analysis of their differential protein activities. First, we computed the differential protein activities for each of n = 112 MOMA subtypes, as assessed by analysis of 20 cohorts with sufficient size for the analysis, including BLCA, BRCA, COAD, GBM, HNSC, KIRC, LAML, LGG, LIHC, LUAD, LUSC, OV, PAAD, PRAD, READ, SARC, SKCM, STAD, THCA and UCEC. Specifically, for each subtype, we generated an average protein activity (consensus protein activity signature) by integrating its VIPER-inferred NES values across each sample in the subtype, using Stouffer’s z score method. Then, we generated a P value to assess the similarity of each HCMI sample to each MOMA subtype in its lineage-matched cohorts. This was accomplished by assessing the following criteria: (1) the enrichment of the top 50 most differentially active protein (top and bottom 25) in the HCMI sample in a protein differentially active in the consensus protein activity signature of each MOMA subtype using NaRnEA; (2) the enrichment of the top 50 most differentially active protein (top and bottom 25) in the consensus protein activity signature of each MOMA subtype in a protein differentially active in the HCMI sample; and (3) by integrating the two using Stouffer’s method.

HCMI samples may show high similarity to more than one subtype. As such, a final assignment was performed as follows. As the tumour and the model samples were derived from the same tumour mass, we assumed that their subtype assignment should be conserved. Thus, among all MOMA subtypes producing a significant match to a model and its parental tumour P < 0.05, one-tailed, Benjamini–Hochberg-adjusted, we selected the subtype with their best integrated Stouffer score. If no agreement was identified (that is, no MOMA subtype with significant OncoMatch score for both the model and its parental tumour), then the two were independently assigned to their best scoring MOMA subtype. For simplicity, we mapped HCMI colorectal, oesophageal–gastric and GBM samples only to COAD, STAD and GBM MOMA subtypes in TCGA, respectively, thus excluding subtypes from READ, ESCA and LGG. A null model for the OncoMatch analysis was generated by generating a probability density of OncoMatch NES scores matching each HCMI sample to all non-lineage-related TCGA subtypes. A one-tail P value was assessed, as the only relevant result would be for the HCMI OncoMatch score to be larger than the OncoMatch from the null hypothesis.

Tumour–model similarity analysis. The transcriptional state similarity between models and their parental tumours was also assessed based on OncoMatch statistics. Differential expression signatures for VIPER analysis were computed by further normalizing and scaling the TPM-normalized gene expression data by subtracting the median expression across all same type samples (that is, tumours or models) in the same HCMI cohort and dividing by the MAD. Again, to prevent division by zero, MAD values <0.01 were set to 0.01. Differential gene expression was assessed separately for tumours and models to avoid sample type related batch effects. For small cohorts (n ≤ 10 samples), we used the median and the MAD of lineage-matched samples that they belonged to. Specifically, the following cohorts were normalized together: (1) intrahepatic, extrahepatic cholangiocarcinoma and hepatocellular carcinoma; (2) all sarcoma types, including bone cancer and desmoid tumours; and (3) tubulovillous adenoma, rare gastrointestinal cancers and colorectal cancer. Similar to the above analysis, we inferred protein activity using metaVIPER with the NaRnEA enrichment analysis algorithm and cancer-type-specific networks. For instance, some HCMI cohorts (for example, COAD) could be matched to multiple lineage networks inferred from TCGA cohorts (that is, COAD and READ). Thus, we assessed protein activities for colorectal, lung, bile duct–liver, brain and gastroesophageal samples in HCMI using COAD–READ, LUAD–LUSC, LIHC–CHOL, GBM–LGG and ESCA–STAD networks, respectively. For rare cohorts that lacked sufficient samples to perform ARACNe analyses, (that is, Wilms tumour, small intestine cancer, gallbladder cancer and unknown carcinoma), we leveraged the ability of metaVIPER to automatically integrate across multiple networks. Specifically, we selected the three networks that produced the greatest differential activity for the top 50 proteins. The rationale is that incorrect networks can only decrease but not increase activity (that is, the more unrelated the network, the smaller the differential activity). To avoid diluting the results based on subpar networks, networks that produced significantly lower differential activity for each protein were excluded from the analysis. We also normalized samples in rare cohorts together with the samples identified as having the best matching networks.

To generate a conservative, nonparametric null-model, we computed the probability density function (PDF) of the NES generated by matching each model with all the non-lineage-matched tumours in HCMI. To improve the PDF estimate, which was highly non-Gaussian, we bootstrapped null-model generation 1,000 times using 60% of the non-lineage-matched samples and used the resulting PDF to convert matched tumour–model pair values to z scores. As both positive and negative NES were integrated, we used Stouffer’s method to generate integrated z scores. Tumour–model pairs with FDR > 0.1 were considered poor matches.

Celligner. To align our HCMI collection to publicly available collections of tumours (TCGA and TARGET) and models (CCLE), we used Celligner³¹, a computational framework to integrate and compare multiple gene expression datasets. Initially, we attempted to expand the tumour and model datasets by directly merging the transcriptional profile of the HCMI collection. However, this approach resulted in poor alignment of the data (Supplementary Fig. 10a). To address this, we expanded on the established Celligner approach by aligning the HCMI dataset to the existing reference datasets in a two-step process. First, we replicated the alignment of TCGA and TARGET tumours onto the CCLE cell lines to create a reference tumour–model map, as described in the original Celligner publication³¹ . Subsequently, we introduced the HCMI dataset into this integrated space.

In brief, using contrastive PCA, we identified gene expression signatures elevated in the TCGA–TARGET tumour collection compared with the CCLE models, consistent with the results reported in the original Celligner analysis (Supplementary Fig. 10b). These signatures were enriched for pathways related to immune and non-malignant features, including stromal cell enrichment (Supplementary Fig. 10b). To mitigate the influence of these non-malignant features, we removed the contrastive principal components (cPCs) associated with these signatures before proceeding to the next steps. We then applied mutual nearest neighbours (MNN) batch correction as part of the Celligner algorithm to align the datasets.

Using this aligned gene expression space as a reference, we integrated the HCMI models and tumours. Specifically, we removed the cPCs increased in the HCMI dataset, which were similarly enriched for immune-related gene signatures (Supplementary Fig. 10c), correlated to DNA-based tumour impurity assessment (Supplementary Fig. 10d) and then we applied MNN batch correction. For this integration, we removed the first three cPCs, which maximized the number of MNN pairs (Supplementary Fig. 10e) and aligned well with the expected variability in the dataset through gene-set enrichment analysis of gene expression profiles. We used a k₁ value of 20 and k₂ value of 50 for MNN, as recommended by the Celligner developers, to reflect differences in dataset size and composition. The results suggested that the HCMI tumour collection exhibited similar levels of variability and contamination as the reference TCGA–TARGET dataset.

After integrating all datasets, we performed PCA on the combined data and calculated pairwise Euclidean distances between all tumours and models in the 70-PC space. These distances formed the basis for all Celligner distance-based analyses that we described in this study. To visualize the integrated data, we created a 2D UMAP embedding based on the 70 PCs for visualization (Supplementary Table 11).

To measure transcriptional relatedness between matched HCMI model–tumour pairs, we used the Celligner distances, which represent the Euclidean distance between a model and its paired tumour in 70-PC space. A P value of paired model–tumour relatedness was then determined by calculating the probability that an unpaired model–tumour distance equal to or smaller than the same model’s distance to all unrelated tumours that were not from the same tissue type. After FDR adjustment, model–tumour pairs with adjusted P values equal to or greater than 0.1 were considered to be models that were non-concordant with their original tumours. Relatedness of tumour lineage is as indicated in Supplementary Table 12.

Expression data input for Celligner. TARGET samples (n = 784) and TCGA expression data (n = 9,806) were obtained from the Xena browser (https://xenabrowser.net). Cell line gene expression data for 1,377 samples were taken from the DepMap Public 19Q4 file. HCMI dataset expression data were downloaded from the GDC portal, and TPM unstranded data were used for only protein-coding genes. Gene expression data were then log₂ transformed after adding a pseudocount of 1. Finally, we subset gene expression data to the 18,550 protein-coding genes that were present in all datasets for all Celligner analyses.

Euclidean and latent TF (OHSU methods) preprocessing. The unique genetic and molecular distribution call for each cancer cohort (for example, pancreatic cancer) and specimen type (that is, tumour or derived model) for these groupings suggested that each be run independently through the following pipeline. We filtered RNA gene expression data for biologically relevant features (the combined set of feature-selected genes of top methods from TCGA trained algorithms)³². For instances where it was needed for analysis, we aggregated multiple TCGA cohorts to more closely match the specific cancer types included in HCMI cohorts (that is, the HCMI lung cancer cohort included TCGA LUAD and LUSC cohorts; the HCMI STAD–ESCC cohort included TCGA GEA and ESCC cohorts; the HCMI kidney cohort included KIRP, KIRC and KICH TCGA cohorts). Any remaining cohorts with low sample size (n < 8) were statistically underpowered and excluded from downstream analysis.

Euclidean distance calculation. Euclidean-based distances were measured in a group (for example, in pancreatic tumours). First, we computed the mean pairwise gene Euclidean distances. For example, we generated a list of pairwise gene distances by calculating all gene distances between pairs of samples. We reported the sample pair distance as the mean of this list. We then repeated the process for all sample pairs, reporting sample pairs for both tumour–tumour pairs for intracohort similarity and tumour–model matched pairs for derived model similarity. We report the z score of these distances and identified outliers as (>3 z score).

Latent TF distance calculation. Latent TF distances were generated using neural networks. Specifically, a variational autoencoder model (NetVae) was trained with gene expression of normal tissues from all TCGA cohorts and genotype–tissue expression (GTEx). These datasets were aligned using quantile ranking. On a gene-wise basis, we calculated the summed difference between a HCMI sample (tumour or model) and normal tissue. These deviation scores were then correlated with mutations and encoded into latent space. The same methods described in the section ‘Euclidean distance calculation’ were applied to these latent values to generate the latent TF distances. Both intracohort and intercohort distances were calculated.

Multiclass pair classification-based distance calculation. We observed that relative gene expression was effective in removing batch effects across different datasets. When we plotted TCGA data using the relative expression of the most variable 500 genes on UMAP space and mapped HCMI samples to the plot, we found that data were clustered on the basis of tissue type. This was in contrast to continuous gene expression data, which clustered data in cohorts. Given this, we used a multiclass pair classification-based approach to compute the distances between matched tumour–model pairs. The multiclass pair classification algorithm¹⁴⁸ uses the relative expression between gene pairs as input and is an extension to multiple classes of the k-top-scoring pairs algorithm¹⁴⁹ that was developed for binary classification. We regressed out the purity effects in both TMP and HCMI data and found that removing purity effects did not make subtypes indistinguishable. We applied the best model selected over cross validations in TMP data to HCMI models and tumours and computed the pairwise Euclidean distances between the assignment probabilities of each HCMI sample to TMP cohort subtypes. We labelled a matched tumour–model distance as an outlier if it was above Q3 + 1.5 × IQR of the tumour–model distances in the TMP cohort.

Canonical parallel direction. We transformed and normalized the combined gene expression matrix using the variance-stabilizing transform from DESeq2 (ref. ¹⁵⁰) then split the matrix into separate tumour and model matrices. We calculated the canonical parallel direction (CPD) by taking the difference between tumour–model pairs and then performing PCA on the resulting difference matrix with the R package irlba (https://github.com/bwlewis/irlba). We took the CPD as the first loading vector. We derived distances by projecting individual expression matrices onto the CPD to obtain a score, and we took the difference in scores between model–tumour pairs as the CPD distance. We labelled samples as outliers when their distance was greater than two standard deviations from the mean.

HCMI and CCLE model coverage comparison

A critical rationale for model generation is to provide proxies for in vitro or in vivo studies of tumour biology. As a result, it is critical to assess what fraction of the tumours in a large-scale repository (for example, TCGA) are associated with effective proxies in large model repositories such as CCLE and HCMI. For this purpose, we assessed the fraction of TCGA tumours that had at least one high-fidelity matched model in CCLE and HCMI. Specifically, we identified high-fidelity models as those producing an OncoMatch score of NES ≥ 10, as previously discussed^38,39.

For this purpose, we generated violin plots representing the probability density of the highest OncoMatch NES scores against each model in HCMI, CCLE and the joint repository of HCMI and CCLE models. As in the other analyses, batch-effect-corrected, TPM-normalized datasets were mitigated using the scGEN algorithm, without considering cancer type and tumour–model covariates. Then, TPM-normalized CCLE gene expression profiles were obtained from DepMap 21Q3 (refs. ^9,151). Gene expression profiles were then centred by subtracting the median expression across all samples and dividing by the MAD. As discussed in previous sections, we used MetaVIPER to assess the activity of TF and co-TF proteins, using regulatory networks generated by ARACNe analysis of the TCGA cohort from which each tumour sample was selected.

Tumour molecular pathology subtype predictions

We predicted the subtype of each sample using TMP models³². These models were trained on molecular data (gene expression, CNV, DNA methylation, miRNA expression and somatic mutations) from TCGA primary tumours. We included a library of models that used either one data type or multiple data types. We applied models that used only gene expression, or DNA methylation, to HCMI samples because of their top performance with TCGA data. We quantile-ranked gene expression data before application. For each sample, we report a single subtype by considering all models and their model confidence values (for example, we report the consensus subtype of the five independent models that were run for ovarian samples).

snRNA-seq analysis

snRNA-seq sample preparation and analysis

Frozen samples (n = 34, from 17 tumour–model pairs) were obtained from the BPC Nationwide Children’s Hospital, DFCI and ATCC. Nucleus isolation was performed using a Chromium Nuclei Isolation kit with RNase inhibitor (10x Genomics, 1000494). Samples were homogenized using a pestle in lysis buffer, passed through a column and centrifuged in debris-removal buffer to eliminate residual tissue and debris. The isolated nuclei were then washed, resuspended and loaded onto a 10x Chromium platform for gel bead-in-emulsion (GEM) generation and barcoding. Following GEM reverse transcription, samples underwent post-GEM RT cleanup and cDNA amplification. After quality control and quantification of cDNA, libraries were prepared using the Chromium Single Cell 3′ Gene Expression Library Construction protocol and sequenced on a NovaSeq platform.

We aligned snRNA-seq short-read data from 17 tumour–model pairs (34 samples) multiplexed across 16 runs, generated by the Columbia University Single Cell Analysis Core, including 7 GBM, 6 PAAD and 4 COAD pairs, to the prebuilt 10x Genomics human reference genome GRCh38-2024-A, using 10x Genomics Cell Ranger software (v.8.0.1)¹⁵². Illumina base call files were converted to FASTQ files with the command cellranger mkfastq. Expression data were processed with cellranger count on the pre-built human reference, which encompassed 38,606 features. Cell Ranger performsd default filtering for QC, and generates filtered feature–barcode matrix files (filtered_feature_bc_matrix.h5; barcodes.tsv, genes.tsv, and matrix.mts) containing unique molecular identifier (UMI) counts for genes for each run.

snRNA-seq demultiplexing

To recover sample-specific profiles for samples processed in the same 10x Chromium flow cell, we demultiplexed the Cell Ranger output¹⁵² using Demuxlet¹⁵³. As Demuxlet requires a list of sample-specific SNPs as an input, we called germline variants for each human sample from the matched normal WGS data. To limit the number of false-positive doublets, which are common in snRNA-seq demultiplexing owing to background RNA and potential large duplications or deletions shared by tumours in the same cohort, we used the gnomADv4 exome data to subset only the most likely germline variants. Specifically, we retained only biallelic SNPs that met the following criteria: (1) population MAF ≥ 0.1%; (2) identified in the gnomADv4 WES cohort¹⁵⁴; and (3) flagged as PASS in gnomAD. We ran Demuxlet with the following parameters: –field GT and –group-list barcodes.txt, where ‘barcodes.txt’ was the cell barcode file generated by Cell Ranger. As previously recommended¹⁵⁵ for demultiplexing snRNA-seq samples, we assigned cells to a sample and used these cells in downstream analyses only if the BEST column of the Demuxlet-generated <Cell Ranger run >.best file started with ‘SNG-’ or the BEST column started with ‘DBL-’ and the PRB.DBL column contained a value of ≤0.99. We classified droplet barcodes with ‘DBL-’ in the BEST column and >0.99 in the PRB.DBL column as doublets, whereas those labelled with AMB were unassigned; droplets labelled as doublet or unassigned were excluded from downstream analyses.

snRNA-seq gene expression analysis

UMI matrices for each demultiplexed sample were processed using Scanpy (v.1.9.3)¹⁵⁶. For QC, we only retained cells with <25% mitochondrial RNA content and UMI counts in the range (800; 50,000). To recover more cells in tumour samples with low cell counts (that is, cases HCM-BROD-0110-C25 and HCM-CSHL-0073-C25), we increased the mitochondrial RNA content threshold to 35% and removed the lower bound on UMI counts. To recover more cells from a run that included three COAD tumours, we removed the threshold on mitochondrial RNA content while maintaining an upper UMI count limit of 100,000. Owing to the low number of cells in the tumour from case HCM-CSHL-0247-C18, we excluded this tumour–model pair from downstream analyses. As a result, we analysed 16 tumour–model pairs: 7 for GBM, 6 for PAAD and 3 for COAD. For statistics on UMI counts before and after QC for individual samples, refer to Supplementary Table 8. UMI counts were normalized, log transformed and scaled, following the standard Scanpy preprocessing workflow. PCA was performed with the ‘arpack’ solver (‘scanpy.tl.pca’). The nearest-neighbour distance matrix was computed with ‘scanpy.pp.neighbors’, setting ‘n_neighbors’ = 15 (default). UMAP embeddings were computed using ‘scanpy.tl.umap’, and PCA and UMAP visualizations were generated with the ‘scanpy.pl.pca’ and ‘scanpy.pl.umap’ functions. For analyses involving dimensionality reduction and visualization of tumour–model pairs or involving multiple cases, the same single-sample workflow was applied by re-scaling the multi-sample concatenated data. Individual samples were clustered using the resolution-optimized Leiden clustering algorithm via the acdc-py wrapper for the Scanpy’s Leiden implementation (https://pypi.org/project/acdc-py/). We determined the optimal number of clusters by varying resolution values from 0.01 to 1.01 in 0.01 increments, selecting the solution that produced the highest average Silhouette score.

Copy number and putative malignant cells detection

We inferred CNAs from expression counts at the single-nucleus level using the inferCNV package¹⁵⁷. We clustered nuclei according to their unsupervised clustering labels, based on gene expression. A reference set of 2,193 nuclei (1,840 oligodendrocytes and 352 microglia) were sampled from four GBM tumours (HCM-BROD-0199-C71, HCM-BROD-0415-C71, HCM-BROD-0012-C71 and HCM-BROD-0002-C71). We used these reference nuclei as controls to infer CNAs in GBM samples. We used a reference set of 3,592 nuclei (2,112 fibroblasts and 1,480 stellate cells) sampled from 6 PAAD tumours (HCM-CSHL-0078-C25 primary and metastatic, HCM-CSHL-0089-C25 primary and metastatic, HCM-CSHL-0073-C25 and HCM-BROD-0110-C25) to infer CNAs in each PAAD sample. We labelled cell types in the references (oligodendrocytes, microglia, fibroblasts and stellate cells) via over-representation analysis (ORA) using the Python decoupler package with canonical human markers (‘brain’ and ‘immune system’ for GBMs; ‘pancreas’, ‘immune system’ and ‘connective tissue’ for PAAD) from the PanglaoDB database^158,159. We sampled a reference set of 435 nuclei (annotated as macrophages by SingleR) sampled from HCM-BROD-0001-C18, and we used the only COAD sample that had enough non-epithelial cells to infer CNAs in each COAD sample. The following parameters were used in each inferCNV run: cutoff=0.1, window_length = 101, HMM=TRUE, mode=’i6’, analysis_mode = ‘subclusters’, denoise=TRUE. To classify putative malignant cells, we compared the inferred CNAs in each snRNA-seq sample to the CNAs detected by WGS in matched bulk samples. Cells were labelled as malignant if they belonged to a cluster with those CNAs.

Cell-type calling for GBM samples

snRNA-seq profiles from GBM samples first underwent coarse-grain cell-type assignment via ORA, using the decoupler Python package with canonical human markers for the brain from the PanglaoDB database^158,159. To address the low signal-to-noise ratio inherent in snRNA-seq data, we performed soft gene imputation by summing the UMI counts of the ten nearest neighbours of each nucleus before running ORA. We labelled nuclei not labelled as ‘oligodendrocytes’ or ‘microglia’ and belonging to clusters of putative malignant cells from inferCNV analysis as ‘malignant’ and assigned these cells to specific cellular states, as described below.

Malignant state assignment for GBM samples

For each GBM nucleus labelled as malignant, a gene expression signature was computed as follows: \({z}_{i}^{(k)}=\frac{{x}_{i}^{(k)}-{\overline{x}}_{\mathrm{ref}}}{{\sigma }_{\mathrm{ref}}}\), where \({x}_{i}^{(k)}\) is the d-dimensional vector of ln(x + 1)-transformed gene expression for nucleus i in snRNA-seq sample k (with d equal to the number of genes), and \({\overline{x}}_{\mathrm{ref}}\) and \({\sigma }_{\mathrm{ref}}\) are the d-dimensional vectors representing the mean and standard deviation of ln(x + 1)-transformed gene expression across all malignant GBM snRNA-seq profiles, respectively. We performed gene set enrichment analysis of each gene expression signature across the six subtypes—NPC1-like, NPC2-like, MES1-like, MES2-like AC-like and OPC-like—as previously identified⁴⁰, using the NaRnEA algorithm¹⁴¹. For this purpose, we used the Python package ‘pyVIPER’¹⁶⁰. Each gene set comprised 39–50 genes that represented distinct cellular states. We then assigned each nucleus to the subtype that produced the maximum NaRnEA NES. Nuclei with the smallest P value (>0.15) were assigned no cellular state and were labelled as unknown. For simplicity, we combined the NPC1-like and NPC2-like subtypes into a single NPC-like subtype; similarly, we merged the MES1-like and MES-2 like subtypes into a single MES-like subtype.

Cell-type calling for PAAD and COAD samples

For each PAAD and COAD snRNA-seq profile, an initial coarse-grain cell-type assignment was performed using the Python implementation of SingleR using the Blueprint-ENCODE reference^161,162. SingleR computes the correlation between each individual nucleus and every sample in the reference. It then labels each nucleus with the cell type with the highest average correlation. To deal with the low signal-to-noise ratio characteristic of snRNA-seq data, before to running SingleR, we applied soft gene imputation using a metacell approach. That is, by adding the UMI counts of the ten nearest neighbours of each nucleus using scanpy.pp.neighbors (Euclidean distance). Nuclei with an ‘epithelial cell’ annotation, which were also classified as putative malignant cells on the basis of inferCNV analysis, were labelled as malignant and assigned to specific transcriptional states, as detailed below. In the COAD cohort, we also retained a small subset of putative malignant cells annotated as ‘neurons’ by SingleR, as these were likely to be malignant cells with neuroendocrine features.

Malignant state assignment for PAAD samples

For nuclei derived from PAAD samples, we used a subtype classification strategy similar to that for GBM. We tested three classification strategies proposed in the literature, including those proposed by ref. ¹⁶³, ref. ³⁴ ref. ⁴³; the first was from bulk profiles and the other were two from single-cell profiles. For the first two, we used the same methodology described for the GBM samples. For the latter, we performed the comparison at the protein activity level. Specifically, protein activity in malignant cells was assessed using metaVIPER¹⁴⁰ by integrating six PAAD gene regulatory networks, which we generated independently from distinct PAAD cohorts, including single-cell profiles (scNET)⁴³, laser microdissected samples (CUMC-net)¹⁶⁴, the TCGA PAAD cohort (TCGA-net)¹⁶⁵, the ICGC PAAD cohort (ICGC-net)¹⁶⁶, the UNC cohort (UNC-net)¹⁶³ and single-cell profiles from HCMI samples (HCMI-net). These networks captured regulatory interactions for TFs, co-TFs, signalling proteins and surface markers. For the first five networks, detailed information for the generation has been previously provided⁴³. The HCMI-net was generated from snRNA-seq profiles from malignant nuclei in the HCMI PAAD cohort using ARACNe3 (ref. ¹⁴¹), with metacell generation constructed using the ‘pyviper.pp.repr_metacells’ function as input (the metacell approach adaptively aims at a target median depth of 10,000 UMIs per metacell whenever feasible). Although this target was not always met, a minimum depth of 8,432 UMIs was reached, which ensured robust representation of the gene expression profiles. To avoid bias associated with different regulon size, all regulons were pruned to the 100 most significant targets based on mutual information analysis. We then classified each single cell on the basis of three main PAAD transcriptional lineages identified as previously described⁴³, including GLS, MOS and PLS. Classification was based on the significance of the NaRnEA-based NES representing the enrichment of the top 100 most differentially active proteins in each single nucleus (50 most activated and 50 most inactivated proteins) differentially active in the GLS, MOS and PLS signatures. Each nucleus was assigned to the cell state with the highest NaRnEA-inferred NES. Nuclei in which the highest-scoring subtype had a P > 0.15 were assigned to no specific transcriptional state and were labelled as ‘unknown’.

Malignant state assignment for COAD samples

For each nucleus in the COAD cohort, we computed a gene expression signature using the same methodology as applied to GBM and PAAD samples. Enrichment of each signature for two previously reported intrinsic subtypes, iCMS2 and iCMS3 (ref. ¹⁶⁷), was calculated using NaRnEA based on the 50 most overexpressed and 50 most underexpressed genes in each signature. Each nucleus was assigned to the cellular state corresponding to the highest NES, whereas nuclei with a P  >0.15 value were labelled as unknown. The vast majority of cells in HCM-CSHL-0143-C20 and HCM-CSHL-0322-C20 were classified as iCMS2, the predominant epithelial state, which is enriched in tumours with strong WNT and MYC signalling activation^167,168. We identified a subset of iCMS3 cells, typically in tumours with substantial immune activation and metabolic dysregulation^167,168, in the HCM-BROD-0001-C18 tumour–model pair (about 50% in the tumour, around 20% in the model; Extended Data Fig. 8b,c). Although the proportions of iCMS2 and iCMS3 cells differed between the tumour and model, the analysis demonstrated that the model effectively recapitulated both cell states that were observed in the parental tumour.

Tumour–model differential gene expression

Differentially expressed genes between tumour and matched model samples were called using the Scanpy function ‘scanpy.tl.rank_genes_groups’ after ln(x + 1) transformation. Statistics were assessed using the Wilcoxon rank-sum test followed by Benjamini–Hochberg correction for multiple-hypothesis testing. We used the resulting signatures to assess pathway enrichment with the ‘pyviper.tl.path_enr’ function in pyVIPER¹⁶⁰. The following gene sets were used: MSigDB Hallmark gene sets for PAAD and COAD¹⁶⁹; MSigDB Hallmark gene sets¹⁶⁹; single-cell GBM expression meta-modules⁴⁰; and bulk-based transcriptional signatures¹⁷⁰ for GBM cases.

Medium-switching assay

To assess transcriptional plasticity driven by culture media, HCM-BROD-0416-C71 GBM cells, originally maintained in NSA medium, were transitioned to formulated-conditioned medium. After initial recovery and expansion in NSA medium, cells were plated on laminin-coated culture vessels and, following 24 h of adhesion, switched to formulated-conditioned medium for either 72 h or 2 weeks. Cells were then fixed and stained for the indicated markers and imaged using the Operetta CLS High-Content Analysis system (Revvity). Full experimental details are provided in the Supplementary Methods.

In vitro drug-sensitivity testing

Patient-derived GBM cells were maintained in NSA medium supplemented with epidermal growth factor and fibroblast growth factor under standard culture conditions^28,51. For drug-sensitivity assays, cells were seeded at a density of 2,000 cells per well in ultra-low-attachment 96-well plates. Twenty-four hours after seeding, temozolomide was dispensed using a D300e digital dispenser in an 8-point titration series (0.03–300 µM). Cell viability was assessed after 5 days using a CellTiter-Glo luminescent assay (Promega).

The HCMI Explorer Suite

To facilitate the exploration of HCMI translational and clinical utility, we developed the HCMI Explorer Suite, a web-based application built using R Shiny (v.1.8.1.1)¹⁷¹. The application was developed in the R programming language (v.4.3.2) and leverages key Shiny packages such as shiny, shinydashboard and shinythemes to provide a dynamic and interactive graphical user interface. The core modules of the HCMI Explorer Suite enable users to explore treatment timelines, model-specific genomic data and transcriptional similarities between cancer models, paired tumours and reference datasets (TCGA and CCLE). All plots were generated using ggplot2 (v3.3.5), with interactive features integrated using plotly (v4.9.3). The Clinical Module uses the swimplot package¹⁷² to visualize patient treatment timelines and to track model treatment exposures. The application is deployed via Shiny Server, with all backend processing and analyses performed server-side. The source code for the app is available at GitHub (https://github.com/human-cancer-model-initiative/HCMI-Explorer-Suite), and the hosted application can be accessed online (https://appshare.cancer.gov/HCMI_Explorer_Suite/).

Ethics statement

All human tissue samples and associated clinical data used in this study were collected as part of the HCMI under protocols approved by the Institutional Review Boards (IRBs) of the respective Cancer Model Development Centres (CMDCs) and participating clinical sites. All procedures were conducted in accordance with relevant ethical guidelines and regulations. Written informed consent was obtained from all participants before sample collection, in accordance with HCMI programme requirements and institutional policies. Each participating centre fulfilled all regulatory requirements, including IRB-approved protocols, informed consent procedures and data-sharing agreements (Supplementary Note 2).

Reporting summary

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