Systematic Bias in Genomic Classification Due to Contaminating Non-neoplastic Tissue in Breast Tumor Samples

Effects of normal contamination in tumors with paired normal tissue

Normal tissue in tumor specimens altered the prediction of tumor subtype or prognosis. Figure 1 shows that in 55 pairs of tumor and normal tissue, the genomic prediction changed as normal tissue contribution to gene expression increased. In the case of the PAM50 subtype predictor (Figure 1A), the samples moved from more aggressive subtypes to less aggressive subtypes as the normal contents increased. The same is true of the ROR-S (Figure 1B). The bias induced by normal contents increased the likelihood of underestimating patient risk. However, the PAM50 was the only assay that showed predictable direction of misclassification.

The other two predictors did not show predictable changes in risk score with increasing percentage normal. With the 70-gene good/poor prognosis and 21-gene recurrence score assays (Figure 1C and 1D respectively), normal contamination biased subtype classification in either direction. For example, 9 samples that were classified as good prognosis using the tumor samples were classified as poor prognosis as normal percentage increased, while for 11 samples increasing normal changed the prediction from poor to good (Figure 1C). The 21-gene recurrence score also caused some samples to move from high to intermediate to low risk, while others moved from low to intermediate to high risk as normal percentage increased (Figure 1D). There was also substantial variation in the genomic classification of the 100% histologically normal tissue. In the final column of Figure 1C and 1D, normal tissue in the 70-gene signature is called poor prognosis in 67 percent of cases and normal tissue in 21-gene signature is high or medium risk in 17 percent of cases.

The differences in predictability of bias direction and the prevalence of high risk phenotype in normal tissue between the PAM50 and the other predictors can be explained by qualitative differences in variability of each gene set across normal samples. Due to the bioinformatics methods used to develop the PAM50 (i.e. intrinsic subtype algorithms[2]), the variation and interquartile range for the PAM50 genes is very low in normal tissues, while it is high across different tumor specimens (Figure 2A & Additional File 1, Figure S1). High variation in tumor samples allows good dynamic range in distinguishing the different tumor specimens from one another, while the low variation in the normal samples means that normal tissue contamination tends to lead to similar effects across samples. Normal tissue appears to have greater variability in the 70-gene (Figure 2B) and 21-gene assays (Figure 2C), even showing similar variability to the tumor themselves in the 21-gene assay.

Effects of Normal Contamination in Public Datasets

Biospecimens of histologically normal tissue adjacent to tumor are rarely available for paired analyses like those performed in the current study, so a generic or prototypical normal profile that characterizes the average effect of normal contamination on tumor predictors could be useful in assessing bias or in adjusting tumor classifications to account for normal contamination. It was possible to define a prototypical normal profile because approximately 90% of samples had normal-like subtype by the PAM50 assay (see Figure 1A) and showed highly correlated gene expression (93% of correlation coefficients greater than 0.6) and low variability across normal samples (Figure 2A). It has previously been reported[11] that a small percentage of normal tissue adjacent to cancer can express a profile similar to invasive ductal carcinomas, and this was also observed in our dataset, but for only a minority (11%) of samples. We excluded these tumor-like normal samples in developing a prototypical normal contamination profile and used 48 normal-like samples to calculate median gene expression. We compared results using this median or protypical cancer-adjacent normal to those obtained using the paired cancer-adjacent normal, and found that subtype classifications under various normal contamination scenarios ranging from 0 to 50% are similar using either the median normal or paired normal tissue (Additional File 2, Figure S2).

Correcting for normal percentages in application of the PAM50

The tumor sample profile is a mixture between the malignant cell profile and its adjacent or contaminating normal profile. Using the prototypical normal profile to deconvolute these two profiles for the only assay with predictable direction of bias (PAM50), we applied Equation 2 (see Methods) to correct the subtype of 24 patients where the percentage of tumor cells, a, was carefully estimated. Five of the 24 samples had original subtype classifications that differed from those after correction for percentage normal (Additional File 3, Table S1). Applying the prototypical normal to these specimens with known percentage tumor did result in changes to subtype with one normal-like sample (with 90% normal contamination) and 3 Luminal A samples (with 70%, 30%, and 25% normal) being reclassified as Luminal B after correction. One Luminal B sample (50% normal) was reclassified as a Her2-enriched (HER2E) subtype. Looking at the ROR-S score, 8 of 24 patients changed in score, with 4 patients switching from Low to Medium or High risk and 4 patients switching from Medium to High risk after adjustment. No patients moved to a less aggressive ROR-S score.

If the true subtype is identifiable by adjusting for normal contamination, then the sensitivity and negative predictive value should increase after correction. Figure 3 shows the change in sensitivity and negative predictive value (NPV) as a function of correction rate in the NKI, Naderi et al., and Wang et al. datasets. In all three datasets, the sensitivity in detecting patients with relapse-free survival < 5 years increased after adjustment for normal contamination and the negative predictive value increased. For each dataset, the maximum % correction shown is consistent with the maximum percentage of normal allowed (for NKI and Wang et al.) or the median % normal reported (Naderi et al.). Sensitivity and negative predictive value are also given for 10-year and relapse-free survival at the end of follow up for all three datasets in Additional File 4, Supplementary Data 1 (excel file). In predicting 5-year relapse-free survival, NPV and sensitivity tended to increase with increasing normal adjustment.

Different correction rates were needed to attain maximum NPV with different datasets. This reflects the fact that contamination rates differed for each of these populations; different quality control criteria were stated for each dataset, and the optimum percentage correction correlated with the stated tumor purity. The NKI dataset had the most inclusive tumor quality control criteria (> 50% tumor nuclei) and required the most aggressive corrections (30% adjustment), while the Naderi et al. dataset with median 60% tumor nuclei required 20% correction and the Wang et al. dataset had the most stringent quality controls (> 70% tumor nuclei) and required the least adjustment (10%). To further confirm that samples with lower contamination rates produce more stable estimates, we used a fourth public dataset of microdissected samples. This fourth dataset had too few samples to perform survival analyses and NPV calculations, but using 48 samples with > 90% tumor cellularity[14], we applied 5% and 10% normal correction. As expected, very few samples changed subtypes with normal correction in this dataset. Only one sample (2% of samples) changed subtype at 5% normal correction and four samples (8% of samples) changed at the more aggressive 10% correction rate. Thus, the misclassification evident using normal adjustment was correlated with reported contamination severity for these datasets.

Normal-contamination correction improves the relapse-free survival curve separation for node-negative patients (Figure 4). Unadjusted Kaplan-Meier curves are shown along with curves adjusted for normal contamination. Additional File 5, Figure S3 shows overall survival in node-negative patients from NKI and Naderi et al., but overall survival was unavailable for the Wang et al. dataset. Across all of these survival analyses, the qualitative trends are the same: (1) normal contamination in tumor specimens causes more aggressive tumors to be misclassified as less aggressive and (2) adjustment for normal contamination improves negative predictive value.

Samples and clinical data

Microarrays

Total RNA was extracted from tumor and paired tumor-adjacent normal tissue by homogenizing tissue and then using an RNeasy RNA extraction kit (Qiagen) to isolate total RNA. RNA quality was assessed using an Agilent bioanalyzer and only samples with RIN > 7 were used for subsequent microarrays. Agilent custom 244 k or catalog 4 × 44 k arrays were performed with linear amplification and labeling according to manufacturer protocol. Two-color Agilent protocols were used for all samples, with Cy3-labeled reference produced from total RNA from Stratagene Universal Human Reference (spiked 1:1000 with MCF-7 RNA and 1:1000 with ME16C RNA to increase expression of breast cancer genes) and Cy-5 labeled patient specimens. New data are publicly available through the Gene Expression Omnibus (GSE22384).

Tumor-Normal Mixing

The impact of increasing proportions of adjacent normal tissue was assessed for each sample using linear combinations of tumor gene expression and paired normal gene expression according to the following:

(1)

where a is the percentage of tumor signal and 1-a is the percentage of normal signal. Combinations were made with increments of a equal to 0.10. The assumption that RNA level has a linear and symmetrical influence on transcript abundance for the PAM50 genes was tested experimentally. Specifically, expression data from pure populations of cancer cells grown in monoculture were linearly combined with expression data from pure populations of breast fibroblasts; expected expression values for ratios of epithelium to stroma equal to 3:1, 2:1 or 1:1 were computed. These computed values were compared to observed gene expression for cocultured cells where the actual percentage of epithelial content was 75%, 66% purity, or 50%, respectively. All data for these analyses are available in the Gene Expression Omnibus, GSE26411) and described in Camp et al. [32]. Additional File 6, Table S2 shows that the correlations between observed and computed were high (0.80-0.92) and that the slopes were very close to 1, suggesting linear and symmetrical influence of epithelium and stroma on transcript abundance.

Normalization and Classification with Genomic Predictors

Data were Lowess-normalized and genes that had a signal of < 10 dpi in either channel were excluded as missing, and genes that had more than 30% missing data across all samples were excluded from further analysis. For genes that passed this filter, missing data were imputed using 10 nearest neighbors. The following genomic predictors were applied to these data based on published methods: the good/poor prognosis signature of Van't Veer et al.[4, 5] that is currently the basis for the Mammaprint assay by Agendia, the recurrence score 21-gene assay that is the basis for Oncotype DX[6, 7], and the recurrence risk (ROR-S) and subtype predictions based on the PAM50 assay[8]. The versions of these assays applied in this paper are based on the published scientific literature and may differ from the patented versions of these assays that are in clinical use. However, subtle differences between patented and published versions of the assay are unlikely to substantially reduce the validity of the sensitivity analyses performed in this work.

For the good/poor prognosis predictor we used the set of probes from the original 70 probes that were available in our platform and in the 55 paired sample dataset. The 70-probe signature contains 60 unique genes, and 95% of these (57 of 60) were on our array. Implementation of the 70 gene good/poor prognosis and 21-gene recurrence assays on NKI data resulted in highly significant prognosis predictions (Additional File 7, Figure S4) suggesting that the current dataset represents tumors with similar clinical characteristics as previous tumor datasets. Applying the PAM50 to the 55 tumor samples gave representatives of all five subtypes, while the paired normal tissue was predominantly of normal-like subtype. A small percentage of non-normal-like subtype specimens was expected based on previous studies that have shown up to 10% of histologically normal specimens from cancer patients may have invasive-like signatures [11]. Of 55 histologically normal samples, 48 (87%) were Normal-like, 6 (11%) were Luminal A, and 1 (2%) was Basal-like. Jackknife distances revealed that the sample with Basal-like subtype was a statistical outlier, suggesting that this normal sample may have had occult contamination with malignant cells that were not identified by pathologist review. The clinical relevance of normal tissue with invasive-like signatures is not established, but it is clear that some genes that are expressed in tumors are also variable in normal tissue.

We next defined a prototypical normal signature that could be used in the absence of a paired normal genomic profile to study the robustness of PAM50 genomic prediction for a given sample. To be used as a prototypical signature, expression must represent the majority of normal samples, thus we selected the 48 samples (representing approximately 90% of samples) that were classified as Normal-like; the prototypical normal was defined as the median gene expression across these Normal-like samples. The observation that some histologically normal tumors may be classified as invasive tumors using genomic predictors has been reported previously [11], so we also considered a prototypical normal based on the median of all 55 samples; at the 30% correction rate, the PAM50 subtype and ROR-S was identical for all samples in the NKI295. The correlation between the two prototypical normal profiles was 0.999. Thus, because the meaning of an invasive-like signature in adjacent non-neoplastic tissue is still poorly understood, we used a prototypical signature based on the 48 normal-like samples to perform adjustments. Using a value of a that represents the percentage of malignant cells in the tumor sample, we can rearrange equation 1, and subtract the effects of normal tissue (1-a) to obtain a corrected tumor signature:

(2)

For the UNC24 dataset, a was set to equal the percentage of malignant nuclei in the sample. For the public datasets, we applied a single value of a for all tumors in a given dataset, with a chosen to maximize negative predictive value, while also ensuring that all ROR-S categories had sufficient sample size to provide stable survival analyses. Analyses of the UNC24 dataset were performed in an attempt to identify gene sets that could be used to predict percent normal in independent data and we evaluated a published percent normal signature [23], but neither had predictive accuracy in independent test sets, suggesting that while gene sets can be identified that distinguish normal from tumor, identification of gene sets that quantitatively estimate percentage normal is a more challenging problem. In the absence of a genomic method for estimating percent normal, we used pathologist assessment where available, or selected a uniform value of a across each dataset/population. With the latter method there is potential to 'over-adjust' or 'under-adjust' individual tumors, however the number of tumors that are adjusted 'more than necessary' will be approximately equal to those that are adjusted 'less than necessary'. Thus, despite the lack of individual-level tumor purity data in most public dataset, the overall population effects and survival estimates are well approximated and the population-level results illustrate the vulnerability of tumor prognostics to normal contamination. Furthermore, the estimated a values selected for each dataset are correlated with the reported quality control standards (see Results) for each dataset, suggesting validity of this approach.

Survival analysis

Among the 295 patients with primary breast carcinoma in the NKI dataset, 151 were lymph node-negative and 36% of these node-negative individuals developed distant metastasis. For the Naderi et al.[13] and Wang et al.[12] datasets, 86 and 286 patients had lymph node-negative disease, respectively, from which 27.9% and 37.4% developed distant metastasis. We used PAM50 to classify all node-negative patients as low, medium or high risk of relapse (ROR-S) and compared these classifications to their metastasis event occurrences. For the NKI dataset 50 of 50 genes in the PAM50 were on the array. Likewise for Naderi et al.[13] and for Wang et al.[12], 48 and 44 of the PAM50 genes were available. We used the Survival package in R to perform univariate survival analyses (Kaplan-Meier Analyses).

Sensitivity and negative predictive value

We used the sensitivity and negative predictive value (NPV) to evaluate the performance of PAM50 as a binary classifier of risk of relapse. We defined sensitivity as the probability of receiving a high or medium risk of relapse score (ROR-S), given relapse occurred by 5 years, by 10 years, or by the end of follow up. Negative predictive value was defined as probability of remaining disease free at a particular time point or at the end of follow up (overall) given a low risk score. The focus of the study was on sensitivity and NPV because these tests are used to guide decisions about chemotherapy and sensitivity (to detect those who will relapse) and negative predictive value (to identify patients who will not relapse) are the highest priority [33].