Unifying heterogeneous expression data to predict targets for CAR-T cell therapy.

Chimeric antigen receptor (CAR) T-cell therapy combines antigen-specific properties of monoclonal antibodies with the lytic capacity of T cells. An effective and safe CAR-T cell therapy strategy relies on identifying an antigen that has high expression and is tumor specific. This strategy has been successfully used to treat patients with CD19 + B-cell acute lymphoblastic leukemia (B-ALL). Finding a suitable target antigen for other cancers such as acute myeloid leukemia (AML) has proven challenging, as the majority of currently targeted AML antigens are also expressed on hematopoietic progenitor cells (HPCs) or mature myeloid cells. Herein, we developed a computational method to perform a data transformation to enable the comparison of publicly available gene expression data across different datasets or assay platforms. The resulting transformed expression values (TEVs) were used in our antigen prediction algorithm to assess suitable tumor-associated antigens (TAAs) that could be targeted with CAR-T cells. We validated this method by identifying B-ALL antigens with known clinical effectiveness, such as CD19 and CD22. Our algorithm predicted TAAs being currently explored preclinically and in clinical CAR-T AML therapy trials, as well as novel TAAs in pediatric megakaryoblastic AML. Thus, this analytical approach presents a promising new strategy to mine diverse datasets for identifying TAAs suitable for immunotherapy.

Introduction

Pediatric acute myeloid leukemia (AML) is a challenging disease associated with high relapse and mortality rates.[1-3]The reported overall survival (OS) of children with AML is 65%–70%, which lags behind that for pediatric acute lymphoblastic leukemia (ALL).[4] In recent years, improvement in outcomes of patients with AML has been largely attributed to a reduction in treatment-related mortality (TRM). Intensive chemotherapy regimens have resulted in increased TRM, leading to decreased OS. Also, patients with certain AML subtypes, such as megakaryoblastic AML (AMKL), carry an inherently poor prognosis.[5] Thus, there is an urgent need to develop novel targeted therapies for pediatric AMKL.

Genetically modified T cells that express chimeric antigen receptors (CARs) combine antigen-specific properties of monoclonal antibodies with the lytic capacity of T cells. CAR-T cells recognize and kill tumor cells in a major histocompatibility complex – independent manner by associating with a specific cell surface antigen, which induces T-cell activation and antigen-positive cellular lysis.[6-8] The ideal antigen target for immunotherapy is an antigen that is expressed on the cell surface of tumor cells, but has low to no expression on normal tissues. CD19, for example, is an antigen expressed on B-cell derived acute lymphoblastic leukemia (B-ALL), with its expression on normal tissues limited to B-cells, limiting its potential toxicities.[9] Because of this, immunotherapeutic strategies using CAR-T cells to target CD19+ malignancies have demonstrated clinical effectiveness with an acceptable toxicity profile.[10-12] Identifying an optimal antigen to target AML has proven difficult, given the marked overlap between antigens expressed on leukemic blasts and normal tissues.[13-16] Examples of antigens currently being targeted include CD33, CD123 and CLL-1. CD33 or siglec 3 is a transmembrane receptor that is expressed on AML blasts but is also expressed on normal HSPCs, lymphoid cells and in Leydig cells.[17] CD123 or interleukin 3 receptor alpha (IL3Ra), is an antigen highly expressed on AML blasts and leukemia stem cells (LSC) that also has low level expression on endothelial cells, and normal hematopoietic progenitor cells.[18] CLL1 or CLEC12A, belongs to the C-type lectin domain superfamily but has limited expression on mature myeloid cells.[19] However, clinical trials targeting CD33, CLL-1 or CD123 in adult and pediatric AML patients are ongoing. As myelotoxicity is a concern, CAR T cell therapy is used a bridge to bone marrow transplantation in certain instances.[20,21]

Publicly available expression data have been used to identify tumor-specific antigens in previous studies. For example, microarray data have been used to mine potential immunotherapy targets in pediatric cancers such as B-lineage ALL (B-ALL) and solid tumors.[22] Similarly, RNA-sequencing (RNA-seq) and proteomics data have been used to predict antigen targets in adult AML patients.[23] An integrative approach that can use data generated from microarray and RNA-seq can maximize the power of discovering immunotherapy targets by increasing the sample size, which is critical for pediatric cancer, a rare disease with limited number of patient samples. Such an approach has not yet been explored due to heterogeneities in sample acquisition and assay platforms. For example, data on the expression of non-disease tissues by transcriptome sequencing (RNA-seq) can be compiled from publicly available resources such as the Genotype-Tissue Expression (GTEx) project.[24] Data on disease tissues such as cancer in children, however, have been generated by RNA-seq through other initiatives such as the Pediatric Cancer Genome Project (PCGP) or the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) project.[25,26] Further, more than 3,000 microarray datasets have been made publicly available in the Gene Expression Omnibus in 2019 alone, which shows that microarray data remain a valuable resource for gene expression profiling despite the emergence of RNA-seq as the main platform for gene expression quantification in recent years. Significant challenges remain in unifying and normalizing these two data resources, given that microarray data measure relative expression whereas RNA-seq data measure absolute read counts.

To enable the comparative analysis of diverse expression datasets in the public domain, we developed a method to transform heterogeneous expression data to maximize the power of predicting tumor-associated antigens (TAAs). We show the application of this approach in identifying candidate TAAs in pediatric AMKL. Our approach to TAA prediction is generally applicable across cancer types and disease states.

Results

Tumor-associated antigen prediction algorithm

Transformed expression values (TEVs) were analyzed using RNA-seq data from primary tumors to identify suitable TAAs: genes absent in physiologically critical tissues and highly expressed in 102 patients with AMKL enrolled in the Pediatric Cancer Genome Project (PCGP).[25] Six additional patients with AMKL enrolled in the NCI TARGET study were also included to confirm that tumor-specific expression of predicted TAAs was not cohort specific. RNA-seq and microarray data from normal human tissue and progenitor cells were used to assess undesirable genes for CAR-T therapy. Genes were implicated as potential antigens according to exposure to the extracellular matrix, which was identified using the definition in the Human Protein Reference Database and Ensembl (Figure 1, Supplemental table 1).[27]

Figure 1.

Workflow for tumor-associated antigen (TAA) prediction

To identify genes with significantly higher expression in the tumor than in normal tissue, Wilcoxon–Mann–Whitney tests were performed using TEVs. Potential TAAs were initially ranked based on the degree of statistical significance associated with tumor-specific expression relative to control tissues. The presence or absence of expression for candidate TAAs was predicted in each primary tumor, provided by the unique detection calling predictions inferred from the control microarray data. However, this information was used only to estimate the percentage of tumors expressing the antigen rather than as a restrictive filter for declaring TAA candidates.

Protein expression was then assessed via immunohistochemistry data by using the following scale: not detected, low expression, medium expression, or high expression, as defined by the Human Protein Atlas.[28] TAAs were excluded as potentially suitable if protein expression was considered to be “low” or greater in more than seven normal tissues. Germline tissues were excluded from this filter, as cancer-testis antigens have been considered effective targets of immunotherapy since their function is generally not critical for survival.[29-32]

Data transformation to improve heterogeneous data comparison

Representing gene expression values as FPKM is standard practice for RNA-seq analysis,[33] and their distribution can vary significantly across datasets such as TARGET and GTEx. However, microarray expression data are distinct from RNA-seq data in that they measure the hybridization value intensities accompanied by detection calls, which represent the probability of absence or presence of gene expression in each sample.[34,35] By normalizing RNA-seq and microarray data such that they can be compared, additional information can be obtained, such as defining gene detection calls derived from microarray data, which can then be applied to RNA-seq TEVs.

We developed a method to transform gene expression values in each sample to a range of 0 to 1 values in both RNA-seq and microarray data to compare relative expression (Figures 1 and 2B). Transformed values in each sample were normalized by quantile normalization, using microarray TEVs as the target distribution, to make the distribution of TEVs equivalent between tumor and control samples (Figure 2c, Supplemental Table 2). In doing so, TEVs and the associated microarray-derived gene detection calls could be compared across platforms.

Figure 2.

Comparison of gene expression value distributions across data cohorts

Table 2.

TAA candidate gene expression status prediction in RNA-seq data

CAR targets	Cancer(% Present)	GTEx(% Present)	WilcoxonP-value
B-ALL
CD19	99	4	5.00E-127
CD22	99	28	2.11E-114
AMKL
ITGA2B (CD41)	100	20	4.48E-70
GP9 (CD42a)	94	5	2.09E-68
CTSW	100	5	1.52E-71
GP1BA (CD42b)	97	2	9.50E-72
ITGB3 (CD61)	94	28	1.08E-58
KCNN4	99	18	2.00E-70
MLC1	100	40	4.78E-38
GPR174	90	1	2.86E-69
RASAL3	100	11	3.42E-69
PRSS21	80	5	2.03E-59
ABCC4	100	17	2.54E-68
AGER	99	18	1.53E-57
PRAME	72	7	4.02E-50
IL2RG (CD132)	100	29	1.32E-64
CD7	85	8	5.70E-58
CD69	100	11	1.18E-71
CD83	100	9	1.30E-71
GP6	80	2	1.85E-68
PTPRCAP (CD45)	100	7	1.16E-69
CD96	80	5	1.41E-61
MS4A2 (CD20)	72	4	7.52E-52
CMTM5	86	31	4.59E-34
GYPA (CD235a)	77	2	6.07E-65
CARD11	88	10	7.30E-53
CD33	89	7	4.94E-59
CD38	87	12	2.10E-61
ITGAX (CD11c)	88	32	2.10E-31

Note: Results of the logistic regression model gene status prediction across transformed FPKM values from patients with AMKL and GTEx control tissue. Assessment of TAA expression status may be indicative of the percentage of observed patients who may or may not be candidates for CAR-T cell therapy by using a specific antigen.

AMKL, megakaryoblastic AML; B-ALL, B-cell acute lymphoblastic leukemia; CAR, antigen receptor; chimeric GTEx, Genotype-Tissue Expression project.[24]

Model-based assessment of gene expression

To confirm our algorithm’s ability to infer gene detection calls strictly based upon TEVs, we built a logistic regression model to predict whether a gene is present or absent based on microarray TEVs. To assess the performance of this logistic regression model, we used a 10-fold cross-validation approach to test 648 bone marrow (BM) microarray samples.[37] Groups 1–9 included 65 samples per group (n = 585), and Group 10 included the remaining 63 samples. The logistic regression model predicted the presence or absence of the gene in the 10-fold cross validation with 90.8 ± 0.36% accuracy relative to microarray gene status prediction based on transformed probe hybridization values (Figure 3a).

Figure 3.

Accuracy of gene expression status prediction based on transformed microarray hybridization values

To further assess the validity of the logistic regression model, a distinct set of six samples derived from the same microarray platform was analyzed to assess the predictive ability of the logistic regression model across data cohorts. Samples were cell sorted into CD133+ and CD133 – cell populations.[38] The present or absent call for each gene from microarray was used as ground truth. The logistic regression model was applied to TEVs, and its prediction was 79.7 ± 1.76% accurate when using TEVs in CD133+ cells (Figure 3a). In CD133 – cells, gene status was predicted with 79.4 ± 0.79% accuracy (Figure 3b). The model correctly identified CD133+ and CD133 – subpopulations solely based on TEVs.

The logistic regression model built on microarray TEVs was then applied to RNA-seq TEVs derived from pediatric patients with B-ALL or AMKL. Common housekeeping genes, as well as those externally defined, were analyzed to assess the model’s ability to predict expected gene presence. The logistic regression model detected gene presence in housekeeping genes with 99.8% accuracy in 304 pediatric patients with B-ALL and 108 pediatric patients with AMKL. To assess the model’s ability to predict the expected absence of gene expression, we also checked the predicted status of tissue-specific antigens. Less than 1% of tissue-specific antigens – which were not expected to be expressed in pediatric B-ALL or AMKL – had predicted presence of gene expression based on TEVs (Table 1).

Table 1.

Modeling of gene expression using TEVs

Expression	B-ALL(% present)		AMKL(% present)	GTEx(% present)
Housekeeping genes
GAPDH	100	100		100
ACTB	100	100		100
CHMP2A	100	100		100
EMC7	100	100		100
GPI	100	100		100
PSMB2	100	100		100
PSMB4	100	100		100
RAB7A	100	100		100
REEP5	100	99		100
Tissue-specific antigens(non–B-ALL and non-AML)
EGFR[39,40,41,42]	0	0		95
MUC16[41]	0	1		9
CEACAM5[42,43]	0	0		7

Note: Presence or absence of gene expression was predicted using a model built from transformed microarray values and their associated gene status predictions. Commonly observed and empirically defined housekeeping genes were monitored to demonstrate the model’s ability to accurately detect gene presence by solely using RNA-seq TEVs.[44] Gene absence was assessed by observing genes that were tissue specific and not expected to be present in B-ALL or AML.

Abbreviations: AML acute myeloid leukemia; AMKL, megakaryoblastic AML; B-ALL, B-cell acute lymphoblastic leukemia; GTEx, Genotype-Tissue Expression project.

When the logistic regression model is applied to primary tumor TEVs, presence or absence of candidate TAAs can be assessed on a patient-by-patient basis. Assessing the presence of a TAA in a cohort of patients may provide insights into the range of patients who may be suitable for receiving CAR-T therapy targeting a specific TAA (Table 2). Conversely, gene status information can exclude potential target antigens or identify a set of patients for whom a particular antigen is not recommended. For example, a TAA that has demonstrated widespread clinical effectiveness as a CAR-T target might not be present in one patient or a small subset of patients, perhaps due to antigen escape in a patient who has relapsed.[45] Our algorithm can detect the presence or absence of a single gene in a single patient, which allows us to detect this type of case before administering treatment.

Performance assessment of suitable antigen targets in pediatric B-ALL

Examples of suitable TAAs for B-cell malignancies include CD19 and CD22.[10,12,46,47] Therefore, our TAA prediction algorithm was first applied to data derived from 196 pediatric B-ALL patients for performance assessment. TAAs that have already demonstrated clinical effectiveness, CD19 and CD22, ranked seventh and sixth, respectively, among the 28 total TAAs predicted by our algorithm.[11,12,47] Several other known B-cell surface markers – CD11a, CD38, CD45, CD69, CD72, CD79a, CD79b, CD132, and CD179b – were also identified as potential antigen targets due to their overexpression in B-ALL.[48-55] Identification of TAAs demonstrating clinical effectiveness in B-ALL validated of our algorithm and provided the proof of concept necessary to apply the algorithm to other diseases (Figures 1 and Figures 4A-B).

Figure 4.

Gene expression profiles of predicted TAAs in B-ALL and AMKL

Identification of tumor-associated antigens in pediatric AMKL

Clinically effective antigen targets with a minimal adverse-effects profile have not yet been identified in pediatric patients with AMKL. We focused on identifying TAAs in 108 patients with AMKL in the PCGP and TARGET databases, given that patients with AMKL tend to have a worse prognosis.[13,14,56,57] Gene expression in normal tissue was tested using 11,688 RNA-seq samples across 53 normal human tissues in the GTEx database.[24] Further, gene expression was tested in progenitor cells using 648 normal human BM samples generated from the Microarray Innovations in Leukemia study, as antigens presented on progenitor cells are not desirable targets for CAR-T therapy.[37]

Our algorithm identified 41 potential TAAs in patients with AMKL (Figure 1, Supplemental table 3), and their profiles were similar to those of clinically effective antigen targets for B-ALL (Table 2). After a manual review of the 41 identified targets, 14 genes were removed as they were expressed in vital tissues or not expressed on the cell surface. Of the remaining 27 AMKL TAAs, 15 were members of cluster of the differentiation (CD) family.[58] We observed a similar distribution of TEVs of established antigens in CAR-T therapy, CD19 and CD22, as well as antigens that have demonstrated promise in AML, CD33 and CD38. Further, we identified a similar, advantageous expression profile in two cancer-testis antigens, PRSS21 and PRAME (Figure 4).[12,17,59-61]

Discussion

As the number and diversity of publicly available gene expression data continue to increase, it is critical to develop analytical approaches that can effectively normalize datasets. Having the capability to compare heterogeneous data can increase the potential and improve the quality of studies such as biomarker discovery. To enable the comparison of gene expression data generated from different datasets/batches and platforms, we reasoned that absolute values may be different among batches or platforms but relative values are conserved between different batches or platforms. We developed an expression value transformation method to reduce data heterogeneity and use it to predict suitable TAAs in pediatric AMKL.

We show that apart from comparisons of TEVs, information on gene presence or absence derived from microarray data is also an effective metric to assess the degree of gene expression in RNA-seq data. We generated a logistic regression model to predict the status of gene expression strictly on the basis of TEVs. This model predicted gene expression status with consistent accuracy by using microarray TEVs in a 10-fold cross-validation as well as assessment in an independent microarray dataset (Figure 3). We then applied this model to determine gene status using RNA-seq TEVs of TAAs for each patient. Our predicted AMKL TAAs had high gene expression in AMKL blasts and low gene expression in control samples, which is similar to what is seen in clinically effective antigen targets in B-ALL (Table 2 and Figure 4). The AMKL TAAs identified had strong expression in at least 72% of patients with AMKL (Table 2). An objective assessment of sample-specific gene expression, as described herein, by using microarray TEVs and the corresponding prediction of gene expression status could be useful to confirm experimental conditions, define an empirical threshold to filter out low gene expression before gene expression analysis, or even select patient cohorts for targeted therapy.

This study demonstrates an application of the transformation methods that can be used to compare gene expression values across data cohorts and even experiment types. We assessed TAAs that could serve as suitable therapeutic candidates for CAR-T cell therapy in B-ALL to confirm that our algorithm can effectively identify clinically effective antigen targets such as CD19 and CD22. We then applied this algorithm to data from pediatric patients with AMKL to identify potentially suitable antigen targets.

Most early phase clinical trials for AML have focused on adult patients and targeted genes in the CD family. Our algorithm successfully identified three genes – CD7, CD33 and CD38 – that are currently being tested in AML clinical trials for CAR-T therapy.[62-65] CD32 and CD123, which are being explored preclinically in AML, were excluded due to the stringency of parameters described herein.[66,67] CD32 was excluded as a candidate TAA due to detectable protein expression in 14 control tissues. CD123 exceeded the maximum median RNA expression value cutoff (0.35) in GTEx control tissues and exhibited detectable protein expression in 45 control tissues.

We identified two novel, potentially suitable cancer-testis antigens, PRSS21 and PRAME, in pediatric AMKL which have a similar expression profile to successfully targeted antigens using CAR-T cell therapy. These genes were overexpressed in primary tumor samples from patients with AMKL, and their RNA and protein expression was almost exclusively limited to testes in normal tissues.[28] Therapies targeting cancer/testis antigens have demonstrated anti-tumor efficacy.[29] Previous research supports that PRSS21 and PRAME are targetable antigens on the plasma membrane.[68,69] Both PRSS21 and PRAME have been associated with tumor progression in testicular tumors and melanomas, respectively, which may explain their AMKL-specific presence.[70] Further, PRAME is known to be recognized by cytolytic T lymphocytes, which is a critical checkpoint in defining antigen potential in CAR-T therapy.[71] Another candidate identified by our algorithm that has proven therapeutic success, potassium-activated channel KCNN4, is currently a Food and Drug Administration–approved drug target in sickle cell anemia.[72]

In conclusion, we successfully generated a data transformation method that improves the ability to probe datasets from multiple gene expression platforms. We identified relatively tumor-specific overexpression of 28 genes in pediatric AMKL (Table 2). Our outlined strategy can also be extrapolated to identify suitable surface antigens as immunotherapy targets in other contexts.

Materials and methods

Data sources

RNA-seq data from 102 patients with AMKL in the PCGP and 6 patients with AMKL in the TARGET data cohort were analyzed to represent the primary tumor population.[56,57] To consider the potential off-target effects of TAAs, 11,688 RNA-seq samples were considered across 53 normal human tissues from the GTEx database.[24] Also, 648 Affymetrix Human Genome U133 Plus 2.0 Array samples generated from the Microarray Innovations in Leukemia (MILE) study in human BM tissue were also used as controls and in building the logistic regression model to predict gene presence in RNA-seq data.[37] Six CD133 cell-sorted Affymetrix Human Genome U133 Plus 2.0 Array samples from GSE7181 were used to evaluate the performance of the logistic regression model when analyzing microarray TEVs and gene expression status.[38]

Potential membrane association was assigned to genes via Gene Ontology terms specified in Ensembl.[27] Relative protein expression data were downloaded from the Human Protein Atlas.

Expression value transformations

Gene expression value transformations were performed on the largest common gene subset among data cohorts being compared. For microarray expression data wherein multiple probes were associated with the same gene, the largest expression value was considered. The transformation methods described below can be applied to gene expression values across platforms to assess relative gene expression. We describe the efficacy of transformations in normalizing heterogeneous gene expression comparisons between data cohorts, as well as between RNA-seq and microarray gene expression data.

Weighted Cumulative Percentage (WCP)

The WCP value transformation method is a fine discriminator of relative gene expression intensity within a sample. The method involves transformation of non-negative log2 raw expression values of descending intensity into a continuous value with the range 0–1, with 0 being the least expression and 1 being the highest expression assigned. The WCP transformation method can be given as

where g represents the cumulative non-negative log2 raw expression value and g represents the transformed value calculated at each iteration. Genes with the same non-negative log2 raw expression value are assigned the same transformed value without increasing g. If not already, the minimum quantile normalized WCP values are set to 0.

Parameter Tuning

Given that CD19 and CD22 have demonstrated clinical effectiveness in B-ALL, we benchmarked transformed values and TAA prediction metrics from these examples to help predict suitable TAAs in AMKL. Thresholds for RNA expression were set to consider only genes with a minimum median TEV of 0.5 in primary cancer samples. In addition, a threshold was set for a maximum median value of 0.3 in normal tissue from the GTEx database. If the gene under consideration is expressed on any tissues that are considered vital for survival, such as the brain or heart, the gene is flagged so that this information is easily accessible upon manual review. However, a strict maximum expression value threshold was not set for the BM dataset, since CD19 expression was high in the BM, likely due to the known role of CD19 in B-cell development, but this did not diminish the suitability of this antigen.[12,73] However, consideration of BM gene expression data was critical in filtering out TAAs such as CCR7 that are not expressed in GTEx normal tissues but are highly expressed in the BM, which contain critical hematopoietic progenitor cells.

Expert Review

After detecting gene targets by using our algorithm, potential TAAs were further analyzed by expert review to confirm their suitability as targets for CAR-T therapy. The default ranking for each gene is defined by an aggregate of negative log10 p-values associated with the degree of significance of tumor expression relative to control tissue, which were multiplied by the difference in median transformed expression value on primary tumors versus controls. However, this order is not a definitive ranking for TAA suitability.

Candidate TAAs were required to be present on the plasma membrane of the intended target cell for recognition by CAR-T cells. Confirmation of extracellular exposure involved extensive literature review and data mining. TAAs were excluded in this analysis if their localization on the plasma membrane was questionable.

Further, literature review helped exclude genes that would likely not serve well as TAAs due to their association with the critical tissue that should not be targeted by CAR-T therapy. For example, RHAG was excluded as a suitable TAA because even though its RNA and protein expression are confined, expression on erythrocytes would likely yield off-target effects that are not clinically tolerable.[28] Another TAA ANK1 initially identified by our algorithm was excluded because direct assay suggested that it was localized on the cytoplasmic side of the membrane, which would likely not be suitable for recognition via CAR-T cells.[74]

Data analyses

Data were analyzed and visualized by using self-developed bash and R scripts. BiocParallel, MASS, pheatmap, and plyr R packages were used within the script for parallel processing, logistic regression modeling, visualization, and object manipulation, respectively.[75-78]

Code Availability

https://github.com/pschreiner/newCAR

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.