Computational selection of antibody-drug conjugate targets for breast cancer

Abstract

The selection of therapeutic targets is a critical aspect of antibody-drug conjugate research and development. In this study, we applied computational methods to select candidate targets overexpressed in three major breast cancer subtypes as compared with a range of vital organs and tissues. Microarray data corresponding to over 8,000 tissue samples were collected from the public domain. Breast cancer samples were classified into molecular subtypes using an iterative ensemble approach combining six classification algorithms and three feature selection techniques, including a novel kernel density-based method. This feature selection method was used in conjunction with differential expression and subcellular localization information to assemble a primary list of targets. A total of 50 cell membrane targets were identified, including one target for which an antibody-drug conjugate is in clinical use, and six targets for which antibody-drug conjugates are in clinical trials for the treatment of breast cancer and other solid tumors. In addition, 50 extracellular proteins were identified as potential targets for non-internalizing strategies and alternative modalities. Candidate targets linked with the epithelial-to-mesenchymal transition were identified by analyzing differential gene expression in epithelial and mesenchymal tumor-derived cell lines. Overall, these results show that mining human gene expression data has the power to select and prioritize breast cancer antibody-drug conjugate targets, and the potential to lead to new and more effective cancer therapeutics.

INTRODUCTION

Personalized cancer therapies are expected to be more effective than conventional treatment and to minimize detrimental effects on normal cells [1]. Different strategies have been used in the development of cancer therapeutic monoclonal antibodies (mAbs) including direct and immune-mediated cell killing, and targeting of the tumor neovasculature [2]. To improve clinical efficacy and to overcome some limitations of first generation mAb-based therapeutics, the industry is currently shifting towards innovative and more powerful modalities such as bi-specific antibodies and antibody-drug conjugates (ADCs) [3, 4]. Antibody-drug conjugates offer the ability to deliver potent cytotoxic drugs specifically to tumor cells [5]. Most current ADC development efforts are focusing on cell surface proteins, the binding of which leads to ADC internalization by receptor-mediated endocytosis, and subsequent release of the cytotoxic payload inside tumor cells [6]. The choice of therapeutic target is a critical success factor in this endeavor [7].

The computational selection of therapeutic targets is a multifaceted process that generally starts with mRNA or protein quantitative analysis to identify targets that are overexpressed in tumor cells. Depending on the type of target and the chosen drug and delivery strategy, candidates are then prioritized using various approaches [8–11]. These approaches include filtering based on subcellular localization, molecular interactions and network modeling, analysis of scientific literature and patents, association of expression with survival, genotype-phenotype analysis, and integration of knowledge from drug and clinical databases. The main criterion for ADC target selection is tumor specificity, to avoid toxicity in vital organs and tissues [7, 12]. To evaluate this property using gene expression data, the analysis of “as many samples as possible” is necessary given the naturally important phenotypic variation within and between individuals [13] and additional noise resulting from sample handling and other experimental procedures [14]. The Gene Expression Omnibus (GEO) [15] is a major functional genomics data repository, currently offering public access to over 100,000 samples analyzed with the Affymetrix Human Genome U133 Plus 2.0 GeneChip [16]. Compiling homogeneous, reliable information using variably granular, semi-structured sample annotations is difficult, while critical for the overall quality of findings in studies reusing public gene expression data [17]. A GEO metadata SQL database was created in 2008 [18], which made records amenable to programmatic analysis. Although sample annotations provide a basis for analysis, various technical problems such as inaccuracy in receptor testing [19] or a low percentage of tumor cells may affect the relevance of a given sample to a given class. Classification techniques [20] are thus essential to validate and refine annotation-based class labels.

Feature selection is another critical aspect of computational target selection. In a recent review, Saeys et al. [21] divided bioinformatics feature selection techniques into three categories depending on if and how the feature search is combined with the classification model. The most common approach to select features in microarray data consists in ranking and filtering features using the Student t-test [22] or the analysis of variance (ANOVA) F-test [23]. Limma, a popular software package used for differential expression analysis, fits a linear model to expression data for each gene, and variance estimates are adjusted by borrowing information across genes [24, 25]. Problems associated with parametric methods include distributional assumptions and the dependence of p-values on sample size [26, 27]. Another class of supervised feature selection techniques makes use of weights acquired in the construction of classifiers such as random forests [28] or support vector machines [29]. These methods, though very effective, do not yield representations that are directly interpretable. In this paper, we used the coefficient of overlap of kernel densities, a concept previously used in social statistics [30], adapted here using a locally adaptive form of the kernel density estimate [31, 32], as a bioinformatics feature selection algorithm. The algorithm is easy to interpret, does not depend on sample size, accommodates various distributions, and is shown to perform equally or better than the above methods in breast cancer and tumor-derived cell line classification problems.

In this study, we have applied feature selection and classification methods to identify candidate therapeutic targets in breast cancer, the most common cancer in women and a heterogeneous disease in nature [33]. Breast cancer is categorized in three basic therapeutic groups associated with distinct molecular subtypes, based on the status of the estrogen receptor (ER), the progesterone receptor (PR), and the receptor tyrosine-protein kinase erbB-2 (ERBB2, a.k.a. HER2) [34]. Although further subdivisions could have been made in each group, we focused our analysis on the molecular subtypes associated with these three basic therapeutic groups (luminal, HER2+ and triple-negative). Over 4,500 breast cancer samples were collected and classified into these three molecular subtypes. For the selection of candidate ADC targets overexpressed in each breast cancer subtype, differential gene expression analysis was performed against over 3,500 samples from a range of vital organs and tissues. Although ADC strategies generally rely on their internalization by cancer cells, a recent report [35] suggests that non-internalizing ADCs targeting the tumor microenvironment may also be effective. For this reason, and also to provide candidate targets for alternative modalities such as antibody-radionuclide conjugates [36], we included both cell surface and extracellular proteins in the analysis. We also aimed to prioritize targets linked with metastasis, since this is the main cause of mortality in patients with solid tumors including breast cancer [37]. Metastasis involves a series of steps where specific tumor cells break through the basement membrane and invade subjacent stromal cell layers, and traverse the endothelium into blood microvessels where they travel to and infiltrate distant sites [38]. The first step in this series of events involves phenotypic changes in subpopulations of cells at the invasive margins of carcinomas, which acquire traits that are important for motility and dissemination, a conversion called the epithelial-to-mesenchymal transition (EMT) [39]. Resistance to therapy and recurrence have been linked with stem cell properties of mesenchymal cells including self-renewal, motility, resistance to apoptosis, cell cycle arrest, suppression of immune responses and enhanced drug transport [40, 41]. Many of the phenomena surrounding EMT and metastasis have been studied in cell line models [42, 43]. Here, we performed classification and differential gene expression analysis in a large collection of tumor-derived cell lines [44, 45], to further prioritize targets linked with the mesenchymal phenotype and metastasis.

RESULTS

Our approach for target selection and prioritization is schematized in Figure 1. In brief, breast cancer samples were classified into three molecular subtypes. Differential gene expression analysis was performed against normal tissues to identify genes overexpressed in each subtype. Subcellular localization information was used in conjunction with gene expression data to select a primary list of cell surface and extracellular candidate targets. In parallel, differential gene expression analysis was performed in epithelial against mesenchymal tumor-derived cell lines to identify, among selected targets, those also potentially linked with EMT.

Figure 1: Overview of the approach for target selection and prioritization. ADC, antibody-drug conjugate.

Breast sample classification

Breast samples (total of 5,379) were initially assigned to one of four classes: normal, luminal, HER2+ and triple-negative, based on sample annotations and receptor status. Class labels were validated using repeated cross-validation combining three feature selection methods, six classification algorithms and two multiclass classification strategies (Figure 2). The performance of all approaches was compared using analysis of variance. The kernel-based feature selection technique slightly surpassed the other two algorithms (p<1E-3). The other factors (multiclass classification strategy, classification algorithm and number of features) all affected performance (p<1E-10). The accuracy under one-against-one (OAO) classification was higher than under one-against-all (OAA) classification. The best performing classification algorithms were: support vector machines (SVM), random forests (RFO) and bagging, followed by k-nearest neighbors (KNN), J48 and naïve Bayes. Accuracy increased as the number of features increased in and remained relatively stable in . Classification accuracy was close to 90% for two of the best performing classification algorithms (SVM and RFO) in combination with the three feature selection methods, with feature numbers in , under OAO classification. This information was used in model selection for iterative ensemble classification. For the final classification of breast samples, filter-based feature selection was performed by selecting the top ranking features using three statistics (q-value from Limma, overlapping coefficient of kernel densities and weight of SVM), and at each iteration classifiers were trained on all labeled instances. In total, 36 predictions were made for each sample (three feature selection methods, two classification methods, and six increasing number of features). Labels assigned with high confidence (>95% of votes) by the ensemble of experts were fed back into the data and used for subsequent feature selection and training of the classifiers. Complete convergence was achieved after 15 iterations. At this point, 5,107 samples were assigned class labels unanimously among experts (100% of votes), and 70 samples were assigned class labels with high confidence (>95% of votes). An additional 82 samples were labeled with reasonable confidence (>75% of votes), and 120 samples were left unlabeled. As a final result, a total of 5,259 samples were labeled, in the following classes: 549 normal, 3,085 luminal, 479 HER2+ and 1,146 triple-negative breast cancer. Among labeled samples, 4,808 (91%) retained original labels assigned using sample annotations and receptor status. Figure 3 shows representative gene expression of markers used for classification in the above classes. The top 64 genes for each of six binary classification problems were selected and hierarchical clustering was performed within each class, and within two groups of genes, the first one consisting of features used for the binary classification problems involving breast cancer against normal samples, (Figure 3, top gene dendrogram), and the second group comprising genes used for classification problems involving breast cancer molecular subtypes (Figure 3, lower gene dendrogram). Distinct gene clusters with low and high expression characteristic of each class are clearly visible.

Figure 2: Repeated cross-validation of class labels assigned using receptor status in breast cancer and normal breast tissue samples. A) One-against-one feature selection and classification. B) One-against-all feature selection and classification. LIMMA, linear models for microarray; OVL, overlap of locally adaptive kernel densities; WSV, weight of support vectors; BAG, bagging; J48, C4.5 decision tree; KNN, k-nearest neighbors, NBA, naive Bayes; RFO, random forests; SVM, support vector machine.

Figure 3: Heatmap and intra-class hierarchical clustering of 5,259 breast cancer and normal breast tissue samples. Top and lower horizontal color bars show sample class and receptor status (ER, PR and HER2), respectively. The top gene dendrogram represents features used for the binary classification problems involving breast cancer against normal samples, and the lower gene dendrogram comprises features used for binary classification problems involving the different breast cancer molecular subtypes.

Cell line classification

Epithelial and mesenchymal cell lines were identified in a collection of 359 cell lines, to perform differential expression analysis and identify candidate targets linked with EMT and metastasis. We used results from a previous study in the NCI60 panel [46] to identify five groups of cell lines having homogeneous gene expression patterns in respect to four major gene clusters. These groups of cell lines were assigned the following class labels: epithelial (10), mesenchymal (23), mixed (12), melanoma (8) and leukemia (6) (Supplementary Data, Table S1). Class labels were validated using repeated cross-validation (Figure 4). The accuracy for each combination of factors was similar to that obtained in the classification of breast samples, although in this case the three feature selection methods performed equally well, and the best performing classification algorithms were SVM and KNN, followed by naïve Bayes, RFO, bagging and J48. Classification accuracy was close to 100% for the two best performing classification algorithms in combination with the three feature selection methods, with feature numbers in , under OAO classification. For the iterative classification of 300 additional cell lines, feature selection was performed independently in three sets of NCI60 replicates, classifiers were trained on each set of replicates, and predictions were made on each unlabeled cell line replicate. In total, 324 predictions were made for each cell line (three sets of labeled replicates, three unlabeled replicates, three feature selection methods, two classification methods, and six increasing number of features). In this case, convergence was achieved after seven iterations. At this point, 266 of the initially unlabeled cell lines were assigned class labels unanimously among experts (100% of votes), and six cell lines were assigned class labels with high confidence (>95% of votes). An additional 20 cell lines were labeled with reasonable confidence (>75% of votes), and the remaining 8 cell lines were left unlabeled. Of the total 359 cell lines, 108 were labeled as epithelial, 88 as mixed, 66 as mesenchymal, 15 as melanoma and 74 as leukemia (Supplementary Data, Table S2). Representative gene expression of markers used for classification in the above cell lines classes is shown in Figure 5. The top 64 genes for each of ten binary classification problems were selected, hierarchical clustering was performed within each class, and within two groups of genes, the first one consisting of features used for the binary classification problems involving the epithelial, mixed and mesenchymal classes only (Figure 5, top gene dendrogram), and the second group comprising genes used for classification problems involving the melanoma and the leukemia classes (Figure 5, lower gene dendrogram). Gene expression was highly homogeneous and characterized by major gene clusters with low and high expression within each class. Two classes (leukemia and melanoma) were dominated by specific tissues and the other classes (epithelial, mesenchymal and mixed) contained cell lines from many different tissues (Figure 6).

Figure 4: Repeated cross-validation of assigned class labels in NCI60 cell lines. A) One-against-one feature selection and classification. B) One-against-all feature selection and classification. LIMMA, linear models for microarray; OVL, overlap of locally adaptive kernel densities; WSV, weight of support vectors; BAG, bagging; J48, C4.5 decision tree; KNN, k-nearest neighbors, NBA, naive Bayes; RFO, random forests; SVM, support vector machine.

Figure 5: Heatmap and intra-class hierarchical clustering of 351 tumor-derived cell lines. Top and lower horizontal color bars show cell line class and tissue origin, respectively. The top gene dendrogram represents features used for the binary classification problems involving the epithelial, mixed and mesenchymal classes only, and the lower gene comprises genes used for binary classification problems involving the melanoma and the leukemia classes. CNS, central nervous system; CT, connective tissue; HLS, hematopoietic and lymphatic system; SM, synovial membrane.

Figure 6: Distribution of tissues of origin and classes in tumor-derived cell lines. CNS, central nervous system; CT, connective tissue; HLS, hematopoietic and lymphatic system; SM, synovial membrane.

Target selection and prioritization

Among 19,674 probesets defined in the Entrez Gene custom CDF [47], 18,282 were matched to human genes from the HUGO Gene Nomenclature Committee (HGNC) database of human genes [48], and 16,811 were of type “protein-coding gene”. From this list of genes, a total of 1,713 genes respected the following conditions and were included in a list of potential cell membrane targets: gene ontology annotation [49] contained GO:0005886 (plasma membrane) or Uniprot subcellular location annotation [50] contained the following term: “cell membrane”, and Uniprot topological domain annotation contained at least one transmembrane and one extracellular domain. In addition, 1,369 genes respected the following conditions and were included in a list of potential extracellular targets: gene ontology annotation contained GO:0005576 (extracellular region) or Uniprot subcellular location annotation contained one of the following terms: “extracellular matrix”, “extracellular space” or “secreted”. To identify genes overexpressed in cancer, gene expression in breast cancer samples was compared with gene expression in normal tissues. For each breast cancer molecular subtype (luminal, HER2+ and triple-negative), genes were ordered according to the mean rank of differential expression in cancer versus normal tissues, and the top 50 membrane and extracellular protein-coding genes were retained (Tables 1 and 2). Batch effects were evaluated using series (experiments) as a surrogate, and were found to be random and relatively inert as compared with the tissue and status effects (Supplementary Data, Figures S1 and S2). Epithelial and mesenchymal marker genes were identified by comparing gene expression in corresponding cell line classes, and filtering genes with a differential expression ratio greater than 2 or smaller than 0.5, and an overlap of kernel densities smaller than 0.6. Four membrane and 13 extracellular targets were also in the mesenchymal gene set (Tables 1 and 2).

Table 1: Membrane ADC targets, mean expression in breast cancer and normal tissues

Table 2: Extracellular ADC targets, mean expression in breast cancer and normal tissues

DISCUSSION

“Ideal” targets, with very high expression in one or more tumor types and very low expression in all normal tissues are rare. The best example of such targets is ERBB2, which is incidentally the only target for which an ADC (trastuzumab emtansine) is currently commercialized for the treatment of breast cancer [51]. In this study, to select candidate ADC targets, we performed gene expression analysis in three breast cancer subtypes versus a range of normal organs and tissues. Our results show that metadata mining and sample classification are instrumental in the assembly of large datasets representative of patient populations, and that feature selection methods and the incorporation of biological knowledge are essential for the selection of clinically relevant targets. Although currently available ADC target data may be too scarce for a formal discussion about sensitivity and specificity, the selection of targets for which ADCs are in clinical development is nevertheless a good indicator of the validity of our approach.

Our list of cell membrane candidates contained one target for which an ADC is already in clinical use and six additional targets for which antibody-drug conjugates are in clinical trials for the treatment of breast cancer and other solid tumors (Table 3). Combining a recent review on ADCs in clinical trials [52] and a search of the https://clinicaltrials.gov/ database [53, 54] revealed that our list of cell membrane targets contained the majority of target antigens for ADCs in clinical development for the treatment of breast cancer, with the possible exception of the tumor-associated calcium signal transducer 2 (TACSTD2), a recently identified ADC target for triple-negative breast cancer [55]. In our analysis, we found that this target, although displaying a high level of expression in breast cancer, also had a relatively high expression in a number of normal organs and tissues including the skin, lungs, and kidneys and did not score high for this reason. According to the list presented in [6], the only candidate target in our list for which an ADC was previously discontinued is mucin-1 (MUC1), for lack of efficacy in ovarian cancer therapy [56, 57]. Clinical efficacy, however, does not depend only on the selected target but also on the design and components of the ADC (drug, linker and antibody). In fact, there is still interest in MUC1: an ADC targeting a specific glycol-epitope of MUC1 (SAR-566658), as well as an anti-MUC1 chimeric antigen receptor (CAR) T cell therapy are currently in clinical trials [54]. Other selected targets in clinical trials for CAR T cell therapy, an approach that also requires highly tumor-specific targets [58–60], include ERBB2, mucin-16 (MUC16), prominin 1 (PROM1) and the prolyl endopeptidase FAP (FAP) [54]. Overall, the high proportion of clinically relevant targets suggests that our target selection method is valid and that our list of targets may contain new candidates with a high potential for ADC development. A number of these are at various stages of pre-clinical research and development, which further validates our findings. Selected examples are discussed below, with a focus on triple-negative breast cancer.

Table 3: Selected targets for ADCs in clinical trials for the treatment of breast cancer and other solid tumors

Target	Drug	Company
ERBB2	T-DM1	Roche/Genentech
GPNMB	CDX-011	Celldex Therapeutics
SLC39A6	SGN-LIV1A	Seatle Genetics
TPBG	PF-06263507	Pfizer/Oxford Biomedica
MUC1	SAR-566658	Sanofi
MUC16	DMUC-5754A	Roche
PVRL4	AGS-22M6E	Astellas

Targeted therapies are currently unavailable for the treatment of triple-negative cancer, and patients in this group have a generally poorer prognosis [61]. Of the cell membrane targets in our list, MUC16 and the cancer/testis antigen 83 (CT83) had the most interesting profiles with high expression in triple-negative tumors and lower expression in all normal tissues examined. An ADC against MUC16 is in development by Genentech, primarily for the treatment of ovarian cancer [62]. CT83, on the other hand, is to our knowledge absent from current ADC development pipelines and has only been recently identified as a potential target in lung cancer [63]. Other genes with high differential expression ratio in triple-negative tumors versus all normal tissues included FAP, the disintegrin and metalloproteinase domain-containing protein 12 (ADAM12) and the low density lipoprotein receptor-related protein 8 (LRP8). FAP is a membrane protein of the serine protease family involved in the proteolysis of the extracellular matrix, which contributes to invasiveness in malignant cancers [64, 65]. In xenograft models, FAP5-DM1 induced long-lasting inhibition of tumor growth and complete regressions in different solid tumors with no detectable side effects [66]. LRP8 has only been recently identified as a potential target in triple-negative breast cancer [67] and is not currently, to our knowledge, considered for ADC development. ADAM12 is involved in a variety of biological processes involving cell-cell and cell-matrix interactions, and is also known as a potential drug target in breast cancer [68]. Interestingly, two of the abovementioned potential targets, namely FAP and ADAM12, are also known to be involved in the EMT [69–71].

In this study, gene expression data mining was preferred as a medium for target selection and prioritization because of the near-transcriptome coverage of modern microarray platforms and the public availability of thousands of human gene expression datasets. Mass spectrometry-based proteomics data analysis would, in contrast, provide better estimates of the quantity of interest (the actual protein abundance), although in a lower number of publically available samples, and at a typically lower resolution. Correlation between mRNA and protein levels was found to be relatively poor in a number of studies [72]. This lack of correlation may be due in part to experimental noise and biases unique to each technique. The quantification of analytes (mRNA and peptides) using indirect signals (probe and peak intensities) may also contribute to this discrepancy. In a recent study, Schwanhaüsser et al. [73] quantified cellular mRNA and protein expression levels and turnover in mouse cells and found that the cellular abundance of proteins was mostly controlled at the level of translation. However, Wilhelm et al. [74] recently analyzed mRNA and protein levels in human tissues and demonstrated that the translation rate is a constant characteristic of a transcript, and that the amount of protein in a cell is primarily controlled by transcription. This result is important: it implicates that differential expression at the mRNA level should correspond to differential abundance at the protein level, although not necessarily in a linear fashion.

Apart from tumor specificity, other factors such as tumor-specific aberrant subcellular localization may influence target selection. The glycoprotein nmb (GPNMB), for example, was characterized with high expression in all breast cancer subtypes versus normal tissues in average, but with relatively high expression in some normal tissues (breast, heart, lung and skin). GPNMB subcellular localization, however, tends to be restricted to intracellular compartments in normal cells, while being enriched on the cell surface in tumor cells [75]. In other cases, normal tissues may be considered expandables in some patients. In our list of candidate targets, some genes including, for example, the gamma-aminobutyric acid receptor subunit pi (GABRP) and cadherin 3, type 1 (CDH3), were characterized by high expression in breast cancer versus normal tissues, but with relatively high levels in normal breast. Toxicity in healthy breast tissue could be a concern and such targets may be appropriate only for subsets of patients undergoing complete mastectomy.

Exploratory gene expression analysis was performed to select candidate ADC targets, for further experimental validation at the protein level, in cell line and animal models and ultimately in clinical trials. Any of the genes in our list, EMT-related or not, respects the fundamental criterion for ADC target selection (high expression in tumor cells and lower expression in normal tissues). Although membrane proteins represent more attractive targets for use in ADC internalization strategies, extracellular proteins may also prove useful, given that strategies such as that reported in [35] are developed and tested. Some of the selected ADC targets would be interesting candidates in triple-negative tumors. These deserve even higher attention for further experimental testing and validation, because triple-negative breast cancer patients have currently no targeted therapy options, and have a generally poorer prognosis. The same approach could be used to mine gene expression data in other cancers, and to identify additional targets for ADC development.

MATERIALS AND METHODS

Microarray data collection

Raw microarray data corresponding to 101,334 samples analyzed using the Affymetrix Human Genome U133 Plus 2.0 platform [16] was obtained from the GEO database [15] using custom Perl scripts. Additional, tumor-derived cell line gene expression profiling data was downloaded from the CellMiner database [76] and the caBIG database [77]. To have a complete set with three replicates for each of the 359 unique cell lines, additional samples from GEO (GSM274690, GSM274785, GSM559851, GSM886956, GSM887076, GSM887415, GSM887651) were added to those datasets. Data analysis was performed using R version 3.1.1 [78] and Bioconductor version 3.0 [79]. CEL files were read in R with the affy package [80] using BrainArray Entrez Gene custom chip definition file version 18 [47], and normalized using the MAS 5.0 algorithm [81]. Pre-processed microarray data was stored in indexed binary files for efficient storage and retrieval.

Metadata analysis

Metadata associated with 101,334 samples (sample id, series id, title, description, source, characteristics) and the corresponding 3,643 series (series id, title, summary, overall design, pubmed ids) were retrieved using GEOmetadb [18] and stored in a SQLite database. Publications (pubmed id, year, journal, title, abstract) linked with the experiments were retrieved using Bioperl utilities [82] and stored in the database. A new table was created for the purpose of sample re-annotation. This table was populated using dedicated Java software comprising a search engine translating Boolean queries into SQL statements and a spreadsheet-like interface allowing direct and programmatic editing of annotations. Breast cancer samples (4,853 tumor) and normal tissue samples (1,067 blood, 291 bone marrow, 526 breast, 353 colon, 50 heart, 279 kidney, 287 liver, 478 lung, 85 pancreas, 334 skin and 25 stomach) were selected from over 200 experiments using this tool (Supplementary Data, GEO Series).

Breast sample classification

Gene expression data was collected for a total of 5,379 breast samples. Receptor status data (ER, PR and HER2) was available in sample annotations for 3,500 samples (1,766 complete, 1,734 partial). For each receptor, two locally adaptive kernel densities were estimated in annotated samples (positive, negative), and receptor status in other samples was predicted by assigning the label corresponding to the maximum posterior [83]. Samples were then assigned to one of four classes: normal (annotated as normal), luminal (ER+ and/or PR+), HER2+ (ER-, PR-, and HER2+) and triple-negative (ER-, PR-, HER2-). Classes were compared one-against-one (OAO) and one-against-all (OAA) for multiclass classification [84]. Filter-based feature selection was performed by selecting the top features ranked using three different statistics: q-values derived from linear models for microarray (Limma) moderated t-test [24, 85], the overlapping coefficient of locally adaptive kernel density estimates [30, 31], and the weight of support vector machines (SVM) [29]. Locally adaptive kernel densities and overlapping coefficients were computed using an in-house R package implemented with Rcpp [86]. The weight of SVM were computed using the e1071 R package [87]. To estimate classification error for combinations of feature selection and classification algorithms, repeated (ten times) five-fold stratified cross-validation was performed [88, 89]. Classification was achieved using six algorithms implemented in the RWeka package [90]: bagging [91], J48 (C4.5 decision tree) [92], k-nearest neighbors (KNN) [93], naive Bayes [94], random forests (RFO) [95] and support vectors machines (SVM) [96]. For the OAO classification, class labels were assigned where the maximum label allocation was reached (three in six binary classification problems). For OAA classification, class labels were assigned where only one label was assigned among four classification problems. For each combination of factors in the cross-validation, the accuracy was calculated as the sum of correct predictions divided by the total number of predictions. The performance of combinations of feature selection methods, classification algorithms, number of features and classification strategies was compared using analysis of variance, and the best performing combinations were retained for ensemble classification [97]. An iterative method [98] was used to assign labels to breast samples. Labels assigned with high confidence (>95% of votes) by an ensemble of experts (36 votes from the combination of three feature selection methods, two classification methods, and six increasing number of features) were fed back into the data and used for subsequent feature selection and training of the classifiers. This procedure was repeated until the number of predictions was stable over a number of iterations, or until complete convergence was achieved.

Cell line classification

Results from Ross et al. [46] were used to label NCI60 cell lines according to patterns of expression in the epithelial, mesenchymal, melanoma and leukemia gene clusters (Supplementary, Table S1). Cell lines having high levels of expression of genes in the epithelial cluster and low expression in the mesenchymal cluster were labeled as epithelial, and cell lines characterized by the opposite pattern were labeled as mesenchymal. Cell lines having intermediate expression profiles between these two classes were labeled as mixed. Cell lines having levels of expression characteristic of the melanoma and leukemia gene clusters were labeled as melanoma and leukemia, respectively. Assigned class labels were validated using repeated cross-validation as described above, with the difference that feature selection was performed separately in three sets of replicates, and validation was performed in the two remaining sets. For the OAO classification, class labels were assigned where the maximum label allocation was reached (four in ten binary classification problems), and for the OAA classification, class labels were assigned where only one label was assigned among five classification problems. For the classification of unlabeled cell lines, starting with the NCI60 panel, each cell line replicate was randomly assigned to one of three sets, and for each set, filter-based feature selection was performed using three statistics as described above. Classifiers were trained on each set of labeled replicates, and predictions were made on each unlabeled cell line replicate. Iterative ensemble classification method was used as described above to assign labels to the new cell lines using 324 votes from an ensemble of experts (combination of three sets of labeled replicates, three feature selection methods, two classification methods, six increasing number of features, and three unlabeled replicates).

Target selection and prioritization

To identify genes specific to or overexpressed in breast cancer, expression profiles in each molecular subtype were compared with gene expression in major organs and tissues were toxicity would likely be a serious concern (blood, bone marrow, colon, heart, kidney, liver, lung, pancreas, skin and stomach). An initial filtering was done to retrieve genes with a maximum ratio (normal/cancer) of 2, an average ratio (cancer/normal) greater than 2 and an average overlapping coefficient [30, 31] smaller than 0.6. Subcellular localization data was obtained from Uniprot annotations [50] and the Gene Ontology Annotation (GOA) database [49]. For each molecular subtype, protein-coding genes were ranked according to the mean rank of ratios in cancer versus normal tissues, and the top 50 (balanced with respect to the three subtypes) cell membrane or extracellular protein-coding genes were retained. Batch effects were evaluated by extracting the first two principal components (52% of total variance) from the 1,000 genes with highest total variance in breast cancer and normal tissues, and visualizing patterns associated with tissue source, sample status and experiments. In addition, batch effects were quantified using linear mixed models [99] with tissue and status as fixed effects and series as random effect, and were compared with the amplitude of differential expression in the breast cancer subtype(s) in which the target was identified versus normal tissues. To prioritize targets linked with EMT, differential gene expression analysis was performed between epithelial and mesenchymal tumor-derived cell lines. The same parameters (expression ratio > 2 and overlapping coefficient < 0.6) were used for the selection of EMT-related targets.

ACKNOWLEDGMENTS

We thank Bernard Massie, leader of the Biologics and Biomanufacturing program at National Research Council Canada, for his support of this project. We also thank the anonymous reviewers for their suggestions and comments.

FUNDING

This project was funded by National Research Council Canada.

CONFLICTS OF INTEREST

The authors declare that there are no conflicts of interest.