Data preparation and filtering

We obtained breast-cancer data from TCGA for 1101 patients in total. To determine germline-mutation status, we downloaded raw sequencing data from CGHub [37] for normal (blood) samples. We limited our analysis to whole-exome sequencing samples that had been sequenced using Illumina Genome Analyzer or HiSeq equipment. Because the sequencing data were stored in BAM format, we used Picard Tools (SamToFastq module, version 1.131, http://broadinstitute.github.io/picard) to convert the files to FASTQ format. We used the Burrows-Wheeler Alignment (BWA) tool (version 0.7.12) [38] to align the sequencing reads to version 19 of the GENCODE reference genome (hg19 compatible) [39]. We used sambamba (version 0.5.4) [40] to sort, index, mark duplicates, and flag statistics for the aligned BAM files. In cases where multiple BAM files were available for a single patient, we used bamUtil (version 1.0.13, https://github.com/statgen/bamUtil) to merge the BAM files.

When searching for relevant germline variants, we examined 61 genes from the BROCA Cancer Risk Panel (http://tests.labmed.washington.edu/BROCA) [41, 42]. We extracted genomic data for these genes using bedtools (intersectBed module, version 2) [43].

We used Picard Tools (CalculateHsMetrics module) to calculate alignment metrics. For exome-capture regions across all germline samples, the average sequencing coverage was 44.4. The average percentage of target bases that achieved at least 30X coverage was 33.7%. The average percentage of target bases that achieved at least 100X coverage was 12.3%.

To call DNA variants, we used freebayes (version v0.9.21-18-gc15a283) [44] and Pindel (https://github.com/genome/pindel). We used freebayes to identify single-nucleotide variants (SNVs) and small insertions or deletions (indels); we used Pindel to identify medium-sized insertions and deletions. Having called these variants, we used snpEff (version 4.1) [45] to annotate the variants and GEMINI (version 0.16.3) [46] to query the variant data.

To expedite execution of the above steps, we used the GNU Parallel software [47]. The scripts and code that we used to process the germline data can be found in an open-access repository: https://bitbucket.org/srp33/tcga_germline.

Geneticists experienced in variant interpretation (BHS, TW, SG, MCK) evaluated each of the germline variants for pathogenicity. Following accepted guidelines for variant classification [42], variants were reviewed individually, followed by a group discussion. In evaluating each variant, we first considered the likelihood of pathogenicity as reported in the University of Washington, Department of Laboratory Medicine clinical-variant database. In addition, we used functional annotations from SIFT [48], Polyphen2 [49], and GERP [50]. When evaluating splice-site variants, we assessed pathogenicity based on whether the variants had been shown experimentally to cause truncations. Additionally, we used NNsplice [51] and Rescue ESE [52] to evaluate splicing variants. We also used maximum entropy modeling [53] and Human Splicing Finder [54], which aggregate data from other splice prediction tools. We classified as benign any in-frame deletions that had been observed as naturally occurring transcripts. When evaluating potential effects of a variant on protein function, we evaluated the extent to which a given protein domain varies within the general population. This process for evaluating candidate cancer-predisposing variants is used in clinical practice and has been demonstrated to maximize actionable results and minimize the frequency of variants of unknown significance [42]. Personal and family histories were unavailable for TCGA patients, so this information could not be included in the evaluation process. The germline calls that we made for these patients were independent of variant-classification calls used in prior studies of TCGA data [25, 55].

To assess loss of heterozygosity (LOH), we used data from Riaz et al. [56]. They had made LOH calls for a large proportion of TCGA breast-cancer patients. Their process included an evaluation of data from Affymetrix SNP 6.0 arrays, genotyping via Affymetrix’s birdseed algorithm, and calculation of the log ratios and B-allele frequencies using PennCNV [57]. The ASCAT algorithm [58] was used to determine allele-specific copy number and loss of heterozygosity for each mutation.

We identified somatic SNVs and indels for each patient by examining variant calls that had been made using Mutect [59]; these variants had been made available via the Genomic Data Commons [60]. Somatic variants that 1) were synonymous 2) snpEff classified as having a “LOW” or “MODIFIER” effect on protein sequence, 3) SIFT [48] and Polyphen2 [49] both suggested to be benign [61], and 4) were observed at greater than 1% frequency across all populations in ExAC [62] were excluded. For BRCA1 and BRCA2, we examined candidate variants based on a priori observations of pathogenicity as reported in the University of Washington, Department of Laboratory Medicine clinical-variant database [63]. Based on these criteria, we categorized each variant as pathogenic, likely pathogenic, variant of uncertain significance (VUS), likely benign, or benign. Then we examined ClinVar [64] for evidence that VUS or likely benign variants had been classified by others as pathogenic; however, none met this criterion. To err on the side of sensitivity, we considered any BRCA1 and BRCA2 mutation to be “aberrant” if it fell into our pathogenic, likely pathogenic, or VUS categories. The final classifications are shown in S1 Table.

Using the somatic-mutation data for each patient, we derived mutation-signature profiles using the deconstructSigs (version 1.8.0) R package [65]. As input to this process, we used somatic-variant calls that had not been filtered for pathogenicity, as a way to ensure adequate representation of each signature. The output of this process was a vector for each tumor that indicated a weight for each signature [19]. S1 and S2 Figs illustrate these weights for two tumors that we analyzed.

We downloaded DNA methylation data via the Xena Functional Genomics Explorer [66]. These data were generated using the Illumina HumanMethylation27 and HumanMethylation450 BeadChip platforms. For the HumanMethylation27 arrays, we mapped probes to genes using a file provided by the manufacturer (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL8490). For the HumanMethylation450 arrays, we mapped probes to genes using an annotation file created by Price, et al. [67] (see http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL16304). Typically, multiple probes mapped to a given gene. Using probe-level data from BRCA1, BRCA2, PTEN, and RAD51C, we performed a preliminary analysis to determine criteria for selecting and summarizing these probe-level values. We started with the assumption that in most cases, the genes would be methylated at low levels. We also assumed that probes nearest the transcription start site would be most informative. Upon plotting the data (S3 Fig), we decided to limit our analysis to probes that mapped to the genome within 300 nucleotides of each gene’s transcription start site. In some cases, probes appeared to be faulty because they showed considerably different methylation levels (“beta” values) than other probes in the region (S3 Fig). To mitigate the effects of these outliers, we calculated gene-level methylation values as the median beta value across any remaining probes for that gene.

To identify tumors that exhibited relatively high beta values—and thus could be considered to be hypermethylated—we used a univariate, outlier-detection algorithm, implemented in the extremevalues R packages (version 2.3.2) [68]. This enabled us to look for extreme values using one side of a specified distribution. This package supports five options for the distribution: normal (default), lognormal, exponential, pareto, and weibull. None of these distributions was a consistently good fit for the methylation (beta) values, in part because the shape of the data differed considerably across the genes (S4–S7 Figs). We used the exponential distribution because it identified hypermethylated genes at rates that were largely consistent with prior work [21]. We used the “getOutliersII” function with default parameter values.

We downloaded copy-number-variation data from the Xena Functional Genomics Explorer [66]. These data had been generated using Affymetrix SNP 6.0 arrays; CNV calls had been made using the GISTIC2 method [69]. The CNV calls had also been summarized to gene-level values using integer-based discretization. We focused on tumors with a gene count of “-2”, which indicates a homozygous deletion.

We used RNA-Sequencing data that had been aligned and summarized to gene-level values using the original TCGA pipeline [25]. To facilitate biological and clinical interpretation, we examined relationships between germline and somatic events and the Prosigna™ Breast Cancer Prognostic Gene Signature (PAM50) subtypes [70]. Netanely, et al. had previously published PAM50 subtypes for TCGA breast cancer samples; we reused this information in our study [71]. We also sought to identify tumors with unusually low expression levels. To do this, we used the getOutliersI function in the extremevalues package. We used the following non-default parameter values: alpha = c(0.000001, 0.000001), distribution = "lognormal", FLim = c(0.1, 0.9). RNA-Sequencing data have been shown to fit the log-normal distribution [72].

We parsed demographic, histopathological, and surgical variables for TCGA samples from the repository prepared by Rahman, et al. [73]. We obtained drug-response data from the TCGA legacy archive (https://portal.gdc.cancer.gov/legacy-archive) and standardized drug names using synonyms from the National Cancer Institute Thesaurus [74].

Quantitative analysis and visualization

To prepare, analyze, and visualize the data, we wrote computer scripts in the R programming language [75]. In writing these scripts, we used the following packages: readr [76], dplyr [77], ggplot2 [78], tidyr [79], reshape2 [80], ggrepel [81], cowplot [82], data.table [83], UpSetR [84], BSgenome.Hsapiens.UCSC.hg38 [85, 86], and Rtsne [87].

To reduce data dimensionality, we applied multidimensional scaling (MDS) [88] to the somatic-mutation signatures and gene-expression profiles. This reduced the data to two dimensions. To quantify homogeneity within a group of tumors that harbored a particular aberration, we calculated the pairwise Euclidean distance between each patient pair in the group and then calculated the median pairwise distance, based on the two-dimensional values [89]. As an additional measure of homogeneity, we used logistic regression to predict BRCA mutation status using the dimensionally reduced coordinates. In evaluating this approach, we used two-fold cross validation and configured the model to weight the minority class in inverse proportion to the frequency of the minority class.

When comparing two given groups, we used a similar resampling approach but instead calculated the median distance between each pair of individuals in either group. To determine whether the similarity within or between groups was statistically significant, we used a permutation approach. We randomized the patient identifiers, calculated the median pairwise distance within (or between) groups, and repeated those steps 100,000 times. This process resulted in an empirical null distribution against which we compared the actual median distance. We then derived empirical p-values by calculating the proportion of randomized median distances that were larger than the actual median distance. We adjusted the empirical p-values for multiple testing using Holm’s method [90]; we applied this correction to the p-values across all aberration types.

For visualization, we plotted the MDS values as Cartesian coordinates; we also used Barnes-Hut t-distributed Stochastic Neighbor Embedding (t-SNE) [91, 92] as an alternative method for visualizing the data.

We created a series of R scripts that execute all steps of our analysis and that generate the figures in this paper; these scripts are available at https://osf.io/9jhr2.

Ethics statement

Brigham Young University’s Institutional Review Board approved this study under exemption status. This study uses data collected from public repositories, other than the germline-variant data, which are restricted by TCGA data-access policies. We played no part in patient recruiting or in obtaining patient consent. We have adhered to guidelines from TCGA on handling data.

Somatic-mutation signatures and gene-expression profiles across all breast-cancer patients

First, we identified the primary somatic-mutation signature associated with each TCGA breast-cancer patient. Approximately 91% of the patients were associated with somatic-mutation signature “1A” (n = 670), “1B” (n = 48), “2” (n = 98), or “3” (n = 130). The remaining patients were assigned to 12 other signatures (S8 Fig). Second, we identified the primary PAM50 gene-expression subtype associated with each patient. The Luminal A and Luminal B subtypes were most common, but each subtype was represented by at least 37 tumors (S9 Fig).

Although it is useful to evaluate the primary somatic-mutation signature or PAM50 subtype associated with each tumor, tumors are aggregates of multiple signatures and subtypes. To account for this diversity, we characterized the tumors based on 1) all 27 somatic-mutation signatures or 2) expression levels for all available genes. To enable easier interpretation of these profiles, we reduced dimensionality of the data using the MDS and t-SNE techniques (see Methods). Generally, tumors with the same primary somatic-mutation signature or PAM50 subtype clustered together in these visualizations (Figs 1 and 2; S10 and S11 Figs); however, in some cases, this did not happen. For example, the dimensionally reduced gene-expression profiles for Basal-like tumors formed a cluster that was mostly separate from the other tumors; but some Basal-like tumors were modestly distant from this cluster, and some Normal-like tumors clustered closely with the Basal-like tumors (Fig 2; S11 Fig). Tumors assigned to somatic-mutation “Signature 3” formed a cohesive cluster (Fig 1; S10 Fig), but some “Signature 3” tumors were modestly distant from this cluster. These observations highlight the importance of evaluating molecular profiles as a whole, not just using the primary category for each tumor.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Two-dimensional representation of somatic-mutation signatures using multidimensional scaling.

We summarized each tumor based on their somatic-mutation signatures, which represent overall mutational patterns in a trinucleotide context. We used multidimensional scaling (MDS) to reduce the data to two dimensions. Each point represents a single tumor, overlaid with colors that represent the tumor’s primary somatic-mutation signature. Mutational Signature 1A (A) was the most prevalent; these tumors were widely dispersed across the signature landscape. Signatures 1B (B), 2 (C), and 3 (D) were relatively small and formed cohesive clusters. The remaining 23 clusters were rare individually and were dispersed broadly.

https://doi.org/10.1371/journal.pone.0239197.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Two-dimensional representation of gene-expression levels using multidimensional scaling.

We used multidimensional scaling (MDS) to reduce the gene-expression profiles to two dimensions. Each point represents a single tumor, overlaid with colors that represent the tumor’s primary PAM50 subtype. Generally, the PAM50 subtypes clustered cohesively, but there were exceptions. For example, some Basal-like tumors (A) exhbited expression patterns that differed considerably from the remaining Basal-like tumors. The normal-like tumors (E) showed the most variability in expression. This graph represents patients for whom we could identify a PAM50 subtype.

https://doi.org/10.1371/journal.pone.0239197.g002

Aberrations in BRCA1 and BRCA2

Of 993 breast-cancer patients with available germline data, 22 harbored a pathogenic SNV or indel in BRCA1; 22 harbored a BRCA2 variant (Fig 3A). We were able to identify loss of heterozygosity (LOH) for all but 3 BRCA1 carriers and 7 BRCA2 carriers (S12 and S13 Figs). A somatic mutation, homozygous deletion, or hypermethylation event occurred in BRCA1 or BRCA2 for 80 patients (Fig 3B–3D). Most of these events were mutually exclusive with each other and with germline variants (S14 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Molecular aberrations in BRCA1 and BRCA2 across all breast-cancer patients.

A) Germline mutations, B) Somatic mutations, C) Copy-number variations, D) DNA methylation levels. SNV = single nucleotide variation.

https://doi.org/10.1371/journal.pone.0239197.g003

BRCA1 carriers fell into the Basal (n = 17); Her2-enriched (n = 1), Luminal A (n = 2), and Luminal B (n = 1) gene-expression subtypes (S15 Fig) [22, 93, 94]. We were unable to assign a gene-expression subtype to one BRCA1 carrier due to missing data. Most BRCA2 carriers fell into the Luminal A subtype (n = 13); the remaining individuals were dispersed across the other subtypes. As demonstrated previously [19], the primary somatic-mutation signature for most BRCA1 and BRCA2 carriers was “Signature 3”; however, other signatures (especially “1A”) were also common (S16 Fig). S17 Fig shows the overlap between these two types of molecular profile.

Homogeneity of somatic mutation signature and expression profiles of germline BRCA1/2 carriers

The somatic-mutation signatures of BRCA1 germline carriers were more homogeneous than expected by chance (p = 0.00056; S18A, S19A and S20A Figs), as were those from BRCA2 carriers (p = 0.0003; S18B, S19B and S20B Figs). As an additional measure of homogeneity, we used logistic regression to predict BRCA aberration status based on the dimensionally reduced data. Using somatic-mutation signatures, we could predict the presence of BRCA1 germline mutations with a sensitivity of 0.72 and a specificity of 0.82. Additional classification results are shown in S2 Table.

None of the three BRCA1 carriers who lacked LOH events clustered closely with the remaining BRCA1 tumors (S18A and S19A Figs). Of the 7 BRCA2 tumors without detected LOH events, 3 clustered closely with the remaining BRCA2 tumors, while 4 did not (S18B and S19B Figs). It has been shown previously that germline BRCA1/BRCA2 mutations leave a recognizable imprint on a tumor’s mutational landscape [19]. This effect may be more consistent when a LOH event has occurred as a second “hit” within the same gene [97], but this was difficult to confirm with the available sample sizes.

Under the assumption that BRCA1/BRCA2 germline variants induce recognizable effects on tumor transcription, we assessed whether tumors from BRCA1/BRCA2 carriers have homogeneous gene-expression profiles. As expected based on the tumors’ primary PAM50 classification, 17 of 21 BRCA1 carriers (for whom we had gene-expression data) overlapped closely with the Basal-like subtype (S15 Fig). As a whole, the expression profiles for this group were more homogeneous than expected by chance (p = 0.0318; S21A, S22A and S23A Figs). However, expression values for BRCA2 carriers were not significantly homogeneous (p = 1.0; S23B Fig). Tumors from these individuals were dispersed across the gene-expression landscape (S21B and S22B Figs).

Similarities among individuals with BRCA1/BRCA2 aberrations

Next we evaluated similarities between BRCA1 and BRCA2 germline carriers. Somatic-mutation signatures for BRCA1 and BRCA2 carriers were highly similar to each other (p = 0.00014; S18A, S18B, S19A, S19B and S24A Figs). However, this pattern did not hold for gene-expression profiles. Although some BRCA2 carriers fell into the Basal-like gene-expression subtype, overall profiles for these patients were dissimilar to those from BRCA1 carriers (p = 1.0; S21A, S21B, S22A, S22B and S25A Figs).

Whether for somatic-mutation signatures or gene-expression profiles, tumors with BRCA1 hypermethylation were relatively homogeneous and highly similar to tumors from BRCA1 germline carriers (Table 1; S18G, S19G, S20G, S21G, S22G and S23G Figs). For gene-expression data, no other aberration type showed significant similarity to BRCA1 germline mutations. Somatic-mutation signatures from tumors with BRCA1 somatic mutations were significantly similar to those from BRCA1 germline mutations (Table 1). Only 2 tumors had BRCA2 hypermethylation, but the mutational signatures for these samples were significantly similar to tumors from BRCA2 germline carriers (p = 0.00054; S19H Fig). Likewise, tumors with a BRCA2 somatic mutation or homozygous deletion had mutational signatures that were similar to germline BRCA2 carriers (Table 1; S18D, S18F, S19D and S19F Figs). Aberrations in BRCA1 and BRCA2 appear to induce similar effects on somatic-mutation signatures—but not necessarily gene expression—whether those disruptions originate in the germline or via somatic events.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Results of similarity comparisons among BRCA aberration groups.

We compared somatic-mutation signatures or gene-expression profiles between groups of patients who harbored aberrations in BRCA1 or BRCA2. We evaluated whether patients in one group (e.g., those who harbored a BRCA1 germline mutation) were more similar to patients in a second group (e.g., those with BRCA2 germline mutation) than random patient subsets of the same sizes. The numbers in this table represent empirical p-values from our resampling approach. In cases where an individual harbored an aberration in both comparison groups, we excluded that patient from the comparison. We used Holm’s method to correct for testing multiple hypotheses.

https://doi.org/10.1371/journal.pone.0239197.t001

Evaluation of aberrations in other cancer-predisposing genes and clinical factors

Next we aggregated all patients who had any type of BRCA1 or BRCA2 aberration into a “BRCA-like reference group.” As a whole, mutational signatures for this group were more homogeneous than expected by chance (p = 0.00001; S26 Fig). We used this reference group to evaluate 59 other cancer-predisposition genes that might be associated with being BRCA-like. For the remaining evaluations, we used somatic-mutation signatures only.

We evaluated whether molecular aberrations in the cancer-predisposition genes resulted in mutational signatures that were similar to our BRCA-like reference group. We found pathogenic and likely pathogenic germline mutations in 13 genes (ATM, BARD1, BRCA1, BRCA2, BRIP1, CDH1, CHEK2, NBN, PALB2, PTEN, RAD51B, RAD51C, SLX4, TP53, XRCC2). The most frequently mutated were CHEK2, ATM, and NBN (Fig 4; S27 and S28 Figs). We found potentially pathogenic somatic mutations in 55 genes, most frequently in TP53, PIK3CA, CDH1, and PTEN (Fig 5; S29 and S30 Figs). Homozygous deletions occurred most frequently in PTEN, CDKN2A, RB1, and CDH1 (Fig 6; S31 and S32 Figs). Hypermethylation occurred in 22 genes, most commonly GALNT12, PTCH1, CDKN2A, and RAD51C (Fig 7, S33 and S34 Figs). Using our resampling approach, we compared each aberration type in each gene against the BRCA-like reference group. In cases where an aberration overlapped between the reference and comparison groups, we excluded individuals who harbored both aberrations. For 11 genes (ATR, BARD1, CDKN2A, CTNNA1, PALB2, PALLD, PRSS1, RAD51B, SDHC, SMARCA4, VHL), at least one type of aberration attained statistical significance after multiple-testing correction (Table 2). A total of 8 aberrations occurred in BARD1: a germline mutation, 2 somatic mutations, and 5 homozygous deletions; as a group, these tumors were statistically similar to the BRCA-like reference group (p = 0.0018). ATR was mutated in 15 tumors and hypermethylated in 1 tumor; together, these tumors were statistically similar to the BRCA-like reference group (p = 0.0035). Tumors with an aberration in PRSS1 were also statistically similar to the BRCA-like reference group (p = 0.0069), but we only observed 2 aberrations in this gene. Other genes showed high similarity to the BRCA-like reference group for one type of aberration only; these included homozygous deletions in CDKN2A (n = 47; p = 0.00001), homozygous deletions in CTNNA1 (n = 6; p = 0.00004), germline mutations in PALB2 (n = 3; p = 0.00001), and homozygous deletions in PALLD (n = 9; p = 0.00001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Non-BRCA germline mutations on the somatic-mutation signature landscape using multidimensional scaling.

Using the same two-dimensional representation of mutational signatures shown in Fig 1, this plot indicates which patients had germline mutations in non-BRCA cancer-predisposition genes. Diamond shapes indicate patients for whom no loss-of-heterozygosity was observed.

https://doi.org/10.1371/journal.pone.0239197.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Non-BRCA somatic mutations on the somatic-mutation signature landscape using multidimensional scaling.

Using the same two-dimensional representation of mutational signatures shown in Fig 1, this plot indicates which patients had somatic mutations in non-BRCA cancer-predisposition genes. Diamond shapes indicate patients for whom no loss-of-heterozygosity was observed.

https://doi.org/10.1371/journal.pone.0239197.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Non-BRCA homozygous deletions on the somatic-mutation signature landscape using multidimensional scaling.

Using the same two-dimensional representation of mutational signatures shown in Fig 1, this plot indicates which patients had homozygous deletions in non-BRCA cancer-predisposition genes. Diamond shapes indicate patients for whom no loss-of-heterozygosity was observed.

https://doi.org/10.1371/journal.pone.0239197.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Non-BRCA hypermethylation events on the somatic-mutation signature landscape using multidimensional scaling.

Using the same two-dimensional representation of mutational signatures shown in Fig 1, this plot indicates which patients had hypermethylation events in non-BRCA cancer-predisposition genes. Diamond shapes indicate patients for whom no loss-of-heterozygosity was observed.

https://doi.org/10.1371/journal.pone.0239197.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of comparisons between the BRCA-like reference group and groups of patients who harbored a specific type of aberration in a candidate BRCA-like gene.

We evaluated whether somatic-mutation signatures from patients who harbored a given type of aberration (e.g., BARD1 germline mutation) were more similar to the BRCA-like reference group than expected by random chance. The numbers in this table represent empirical p-values from our resampling approach. In cases where no patient had a given type of aberration in a given gene, we list “N/A”. The “Any” group represents individuals who harbored any type of aberration in a given gene. We used Holm’s method to correct for testing multiple hypotheses.

https://doi.org/10.1371/journal.pone.0239197.t002

Lastly, we evaluated whether the following data types were correlated with being BRCA-like: 1) unusually low mRNA expression in a given gene, 2) demographic, histopathological, and surgical observations, and 3) patient drug responses. First, we calculated the median Euclidean distance—based on somatic-mutation signatures—between each patient and the BRCA-like reference group. Then we used a two-sided Pearson correlation test to assess the relationship between these median distances and each candidate variable. In determining whether a tumor exhibited unusually low mRNA expression for a given gene, we used an outlier-detection technique (see Methods). Unusually low expression of BRCA1 (rho = 0.22, p = 0.0024) and RAD51C (rho = 0.20, p = 0.016) showed the strongest positive correlation with the reference group, whereas CDH1 (rho = -0.19, p = 0.023), PIK3CA (rho = -0.19, p = 0.025) and BARD1 (rho = -0.19, p = 0.028) showed the strongest negative correlation (S35 and S36 Figs). Triple-negative status and infiltrating ductal carcinoma histology were the most positively correlated clinical variables (S37 Fig). No chemotherapy treatment was significantly associated with being BRCA-like, though availability was limited for the drug data (n = 211; S38 Fig). Relationships among low mRNA expression and these factors are likely interwoven. For example, CDH1 expression has been associated with molecular subtypes such as triple-negative status [98]. Larger studies will be necessary to disentangle these effects.