Can the Pareto optimality theory reveal cellular trade-offs in diffuse large B-Cell lymphoma transcriptomic data?

doi:10.21203/rs.3.rs-3467629/v1

Download PDF

Article

Can the Pareto optimality theory reveal cellular trade-offs in diffuse large B-Cell lymphoma transcriptomic data?

https://doi.org/10.21203/rs.3.rs-3467629/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

One of the main challenges in cancer treatment is the selection of treatment resistant clones which leads to the emergence of resistance to previously efficacious therapies. Identifying vulnerabilities in the form of cellular trade-offs constraining the phenotypic possibility space could allow to avoid the emergence of resistance by simultaneously targeting cellular processes that are involved in different alternative phenotypic strategies linked by trade-offs. The Pareto optimality theory has been proposed as a framework allowing to identify such trade-offs in biological data from its prediction that it would lead to the presence of specific geometrical patterns (polytopes) in e.g. gene expression space, with vertices representing specialized phenotypes. We tested this approach in diffuse large B-cell lymphoma (DLCBL) transcriptomic data. As predicted, there was highly statistically significant evidence for the data forming a tetrahedron in gene expression space, defining four specialized phenotypes (archetypes). These archetypes were significantly enriched in certain biological functions, and contained genes that formed a pattern of shared and unique elements among archetypes, as expected if trade-offs between essential functions underlie the observed structure. The results can be interpreted as reflecting trade-offs between aerobic energy production and protein synthesis, and between immunotolerant and immune escape strategies. Targeting genes on both sides of these trade-offs simultaneously represent potential promising avenues for therapeutic applications.

Health sciences/Oncology/Cancer/Cancer models

Biological sciences/Cancer/Cancer models

Biological sciences/Cell biology/Mechanisms of disease

Biological sciences/Molecular biology/Transcriptomics

Biological sciences/Systems biology/Genetic interaction

Pareto theory

transcriptomics

lymphoma

oncology

archetypes

trade-offs

optimality

systems biology

Oncology has made tremendous progress over the last few decades (Cancer Progress Timeline | ASCO). New treatment modalities have been developed and existing ones have been refined, combined and optimized to the extent that clinical outcome, including overall survival, has improved significantly for many tumour types (e.g.¹). Yet, in many cases, despite these advances, and although response to therapies is often initially good for many types of cancer, resistance develops and leads to relapse/refractory tumours ². For too many patients, a cancer diagnosis is still synonymous with premature death and cancer still represents one of the main causes of mortality and morbidity in human populations.

With the development of “omics” sciences, the last few decades have also seen the accumulation of a huge amount of publicly available data, especially DNA and RNA sequences, as well as data about numerous signalling pathways and cellular processes. We are therefore facing a situation where the amount of available data is no longer the main limiting factor in advancing cancer therapy. Unfortunately, the availability of such “big data” has yet to propel cancer research to dramatic new therapeutic advances or to cure cancer. Indeed, immunotherapy is probably the most significant innovation in cancer treatment of the last decade and its development does not directly derive from the omics revolution but rather from a more conventional incremental increase in our knowledge of immunology ³.

Arguably, one of the main reasons explaining this inability to take full advantage of the information contained within these public databases stems from the lack of a theoretical framework, organizing principles and methodology allowing to target only their most biologically relevant elements (thus turning data into information). Such a theoretical framework should provide a way to understand the most crucial and essential tasks that malignant cells must perform in order to survive and the genes used to perform these tasks, with the ultimate goal of selectively eliminating malignant cells.

Typically, cancer is understood to be the result of the accumulation of random mutations in genes involved in processes such as cell cycle control, apoptosis and DNA repair, eventually producing the malignant phenotype recognized as cancer ⁴. This view, often called the “somatic”, “reductionist” or “somatic mutation theory (SMT)” does not offer any specific guiding principle to navigate omics data ⁵. Under this classical paradigm, basic cancer research aims at providing increasingly detailed descriptions of various genetic alterations and their impact on cellular pathways, but with no reference to a deeper organising theoretical framework to increase our understanding of the phenomenon and thus becomes an essentially descriptive exercise ⁶.

In order to solve what is currently recognized as one of the main problems in cancer therapy, namely the development of resistance to previously effective treatments ^{7, 8}, an alternative approach might be to try gaining insight into cancer cell vulnerabilities by looking at the fitness of different genotypes through analysis of the genomic fitness landscapes. In principle, this exercise should allow the prediction of cancer evolution under various circumstances. The genomic fitness landscape allows not only the determination of the selection coefficient and gradient associated with any particular allele, but also the identification of critical biological functions in the malignant phenotype and the prediction of response to therapy. In practice however, this is a daunting task as it requires within patient serial measurements at single-cell genomic resolution ⁹, which is currently out of reach for most research teams.

Interestingly, however, an approach from systems biology may alleviate the need to characterize tumors’ fitness landscapes to identify biologically critical tasks and functions performed by tumor cells. This idea relies on the fact that living organisms and cells need to optimize multiple tasks simultaneously and inevitably face trade-offs in doing so. Mathematically, it can be shown that such constraints imply that the data will form polytopes in trait space (e.g. gene expression space) where vertices of the polytope each represent the optimal specialist phenotype for performing a given task ¹⁰. Known as Pareto task inference (implemented in the ParTI package in MatLab available at: Pareto Task Inference (ParTI) method | Uri Alon (weizmann.ac.il)), this idea comes from the general Pareto optimality theory that has been used in fields like engineering and economics for many decades in order to solve multi-task optimization problems involving trade-offs between different objectives. Briefly, ParTI allows the use of algorithms to find the best fitting simplex (the simplest polytope) on the data after dimensionality reduction by principal component analysis (PCA), and to test the statistical significance of the simplex by means of the t-ratio test. This test computes the ratio between the area encompassed by the fitted polytope to the area defined by the convex hull of the data (the t-ratio) and then repeat this process after randomizing the data in order to obtain an empirical null distribution for the t-ratio. Pareto task inference has so far been applied to a few biological datasets, including solid tumor gene expression ^{11, 12}. Most importantly, by circumventing the need to reconstruct the fitness landscape, this method has the potential to directly unveil the most relevant cellular processes for therapeutic intervention by allowing one to infer the fitness consequences of disrupting the cell’s ability to carry out specific tasks.

Although some authors have raised the possibility that the ParTI algorithm might produce spurious results due to phylogenetic, ancestral or other population structure ^{13, 14}, Adler et al. ¹⁵ tested whether such problems could affect cancer gene expression data and found no evidence that this was the case. The risk of false-positives due to “p-hacking” through data preprocessing has also been raised ¹³. Amid these concerns, rather than speculating on possible artefactual polytopes or on the possibility of statistical false-positives, a more fruitful way of addressing the validity of the ParTI method is to derive specific predictions of the Pareto optimization theory besides the geometrical prediction of polytopes in trait space. In particular, the biological interpretability and underlying logic of the results obtained from that method appear to be of particular relevance in that regard.

For instance, if the most specialized cellular phenotypes (the “archetypes”) identified as the vertices of the polytope in gene expression space are the results of artefactual data structure unrelated to biological specialization and optimization, one would not expect to find statistically significant and distinct enrichment in particular biological functions at these archetypes. Moreover, if polytope formation results from optimization trade-offs, as per the theory, the various tasks defining the archetypes should all be considered important, if not essential, for cell survival. Indeed, trade-offs imply that some of these tasks cannot be simply eliminated completely from the cell’s functional repertoire, otherwise there should be no trade-offs constraining the optimization process. Thus, if the existence of the archetypes identified by the ParTI algorithm is the results of optimization trade-offs as predicted by the Pareto optimization theory, we should expect these archetypes to be statistically enriched in specific biological functions, and to be characterized by varying proportions of some shared tasks/genes in different combinations.

Because of its potential as a useful theoretical principle allowing meaningful analysis of omics and other biological data, and thus to provide important information about both basic biology and therapeutically actionable knowledge, the Pareto task inference approach should be carefully and independently evaluated and validated. So far, the ParTI package has been used on a few solid tumor datasets but, to our knowledge, never on hematological malignancies. Because of their cell type (malignant lymphocytes) and phylogenetic (from clonal proliferation) high level of homogeneity, data from lymph node biopsies of lymphoma patients are particularly well suited to study the validity of this method. In order to test the above predictions, we analyzed transcriptomic data from diffuse large B-cell lymphoma (DLBCL) samples for 481 patients using ParTI code. Archetypes identified by the algorithm were then tested for functional enrichment and compared in terms of gene composition and representation. Results were ultimately evaluated for their interpretability, potential biological meaning and usefulness for generating hypotheses.

Gene expression data

Data were obtained through the National Cancer Institute’s Genome Data Commons Data Portal (GDC (cancer.gov)). Gene expression quantification expressed as transcripts per millions (TPM) from RNASeq was available for 517 DLBCL patient-lymph node biopsy samples (cases.disease_type in ["mature b-cell lymphomas"] and cases.primary_site in ["lymph nodes"]) from projects NCICCR-DLBCL, TCGA-DLBC, and CTSP-DLBCL1. Of these, 481 (93%) were from the NCICCR-DLBCL project, while only 29 (5.6%) and 7 (1.4%) were from the TCGA-DLBC and CTSP-DLBCL1 projects respectively. In order to ensure maximum homogeneity in sampling and sequencing protocol and hence between-sample comparability of data ¹⁶, we analyzed the 481 cases from the NCICCR-DLBCL project. We kept all informative genes/isoforms, defined as having standard deviation and variance ≥ 1 TPM, for downstream analyses. To avoid p-hacking issues, we refrained from performing any preprocessing transformations/manipulations of data and used the raw TPM gene expression quantification. TPM is considered a better unit of RNA abundance than RPKM/FPKM since it respects the invariant-average property (the average TPM is a constant equal to 10⁶ divided by the number of annotated transcripts) and is proportional to the average RNA molar concentration (rmc). It has thus been adopted by the latest computational algorithms for transcript quantification such as RSEM, Kallisto and Salmon ¹⁶. The resulting data matrix had 481 sample lines by 22 250 genes/isoforms columns.

PCA and Pareto Task Inference analyses

Because the ParTI code package relies on principle component analysis (PCA) for dimensionality reduction, as proposed by Vieira ¹⁷, we first tested the validity of performing PCA on the data by computing the Ψ and ϕ statistics, the number of significant principal components, as well as the number of genes/isoforms with significant correlations with each of the principal components, by permutations and bootstraping using PCAtest in R ¹⁸.

Following Mikami and Iwasaki ¹⁴, the simplex best fitting the data was determined by the SDVMM algorithm ¹⁹. Five such fitting algorithms are available in ParTI. The SISAL algorithm is not recommended for datasets of less than about 1000 data points (or more precisely of less than \({10}^{N}\) points, with N being the number of PCA dimensions), and will estimate the archetypes outside of the convex hull of the data, while potentially leaving important points outside of the fitted polytope and t-ratio test analysis (PartiCode homepage: Pareto Task Inference (ParTI) method | Uri Alon (weizmann.ac.il) accessed 11/03/2023), which may generate false positive or false negative results of the t-ratio test. The MVA and MVE algorithms, which also locate the archetypes outside of the convex hull, could be used with smaller datasets of < 1000 points, but are not robust to noise and outliers. Contrary to the previous algorithms, the PCHA algorithm estimate the archetypes within the data but also suffer from a susceptibility to noise and outliers. Thus, we agree with Mikami and Iwasaki ¹⁴ that the SDVMM algorithm is the only option combining a strict data constraint with a robustness to noise and outliers and seems clearly preferable for statistical testing. The default value of eight PCA dimensions, as per the ParTI’s example file “exampleCancerRNAseq.m”, was kept as input.

Archetype coefficients in original gene expression space were then extracted from ParTI’s output and genes ordered according to coefficient value. Genes with the highest positive and negative coefficient values for each archetype (genes with the highest PCA loadings at the archetype location and thus defining the archetype in terms of most positively/negatively correlated expression values) were identified by plotting ordered coefficient values and visually identifying the abrupt change in slope between the vast majority of genes with coefficients close to zero and the few genes with highly positive (highly positively correlated) or negative (highly negatively correlated) values, and taking this inflexion point as cut-off. These archetype-defining gene lists were then tested for functional enrichment using Gene Ontology tools (GO aspect: biological process) ^{20, 21, 22}. Functions with a Fisher’s exact test enrichment p-value < 0.05 (after adjustment for false-discovery rate (FDR)) were considered significantly enriched at the archetype. The composition of significantly enriched functions was compared and the level of shared vs unique functions among archetype analyzed using Circos ²³.

Data Availability

All data are publicly available on the Genome Data Commons Data Portal at GDC (cancer.gov).

PCA validation

PCAtest was run for 100 permutations and bootstrapping replicates and estimated an empirical Ψ value of 51848193.065, with a max null Ψ = 979604.238 and a min null Ψ = 978848.881, for a p-value < 0.00001. Likewise, the ϕ value was estimated at 0.337, with a max null ϕ = 0.047 and min null ϕ = 0.047, for a p-value < 0.00001. These results indicate highly significant non-random correlations in the data, justifying the use of PCA and confirming that the analysis is biologically meaningful. Moreover, the analysis showed that the first eight principal components (PC) explained 80% of the total original variance (from 1.8–33% for individual PC). Finally, the number of genes with significant correlations with each of these eight principal components ranged from 2830 to 7546. Together, these results strongly support the use of PCA on this dataset.

Polytope fitting

The elbow method applied by the PartiCode package suggested the presence of four archetypes, thus defining a tetrahedron (3-d simplex) (Fig. 1). The identified polytope was highly significant with the t-ratio test indicating a p-value < 0.00001.

Archetype-defining genes selection

For each of the four archetypes, ordered coefficients were plotted and revealed a clear pattern of abrupt change from coefficients close to zero for the vast majority of genes, to a few genes showing strongly positive or negative coefficients (Supp. Figure 1). Archetype-defining genes were therefore selected by considering a horizontal line (slope = 0) for genes with a coefficient close to zero and a line with a slope of 1 (vertical) for genes with strongly positive or negative coefficients. These lines were extrapolated and the line bisecting these right-angled extrapolated horizontal and vertical lines at 45 degrees was determined. The mid-point of this 45-degree line, when positioned so it just touched the data, was considered as the inflexion point (elbow point) and used as the coefficient cut-off to select the positive and negative archetype defining genes (Supp. Figure 2).

Enrichment analysis

The above procedure identified between 23 and 132 archetype-defining genes with positive loadings, and between 12 and 85 archetype-defining genes with negative loadings depending on the archetype (Supp. Table 1). For each archetype, both positive and negative gene lists were submitted to the Panther database for enrichment analysis. All eight archetype-defining gene lists (one positive and one negative loading list for each of the four archetypes) were significantly enriched (after FDR adjustment) in some functions. The 10% functions with the most significant enrichment p-value for each of the eight lists are presented in supplementary table 2. GO terms associated with each individual archetype-defining genes are listed in supplementary table 3.

The first archetype was positively enriched in aerobic energy production (three most significant GO terms: oxidative phosphorylation/aerobic respiration/aerobic electron transport chain) and negatively enriched in protein translation (three most significant GO terms: cytoplasmic translation/translation/peptide biosynthetic process). The second archetype was positively enriched in protein translation (three most significant GO terms: cytoplasmic translation/translation/peptide biosynthetic process) and negatively enriched in energy production and immune functions (three most significant GO terms: oxidative phosphorylation/adaptive immune response/cellular respiration). The third and fourth archetype were positively enriched in a mixture of aerobic energy production and immune functions, and negatively enriched in immune functions. The 10% lowest FDR-adjusted p-values of the enrichment analysis were very small, ranging from 0.0025 to 5.67x10^− 95, indicating that, from a functional point of view, archetype-defining genes represented highly non-random genomic sub-samples. As predicted, a significant proportion (49%) of archetype-defining genes were shared between at least two archetypes, albeit in different combinations and proportions (Fig. 2).

Biological interpretability

One of the most striking features of the functional enrichment analysis was the contrast between archetype 1, which showed specialization for oxidative phosphorylation (OXPHOS) at the expense of protein synthesis, and archetype 2 which was specialized in protein synthesis at the expense of both OXPHOS and immune functions (Fig. 3a), b)). A detailed analysis of the genes with expression positively correlated with archetype 2 and negatively correlated with archetype 1 revealed 34 genes involved in protein synthesis. These genes consisted of various ribosomal subunits and two translation elongation factors, as well as CHCHD2. Interestingly, it was recently shown that in stress conditions produced by carbonyl cyanide m-chlorophenylhydrazone (CCCP) treatment, a decoupling agent known to induce oxidative stress ²⁴, CHCHD2 knockdown triggered the integrative stress response (ISR) in cultured HeLa cells ²⁵. The role of the ISR is to maintain homeostasis under stress conditions, including oxidative stress, and is known to slow down protein synthesis via the phosphorylation of eIF2α ²⁶. Indeed, the ISR was recently observed to inhibit protein synthesis in the context of oxidative stress from mitochondria-derived production of reactive oxygen species (ROS) in a cardiac ischemia/reperfusion model ²⁷. Because of the high rate of energy production via OXPHOS and the associated ROS generated, archetype 1 cells may down-regulate CHCHD2 in order to restore homeostasis via the ISR, at the expense of protein synthesis. In contrast, archetype 2 cells may overexpress CHCHD2 in order to benefit from its anti-apoptotic effect ²⁸. Recently, the importance of CHCHD2 in cancer and its potential as drug target was highlighted by Gundamaraju et al. ²⁹. Consistent with the hypothesis of a driving role for increased ROS production in archetype 1, the expression of ROMO1 was negatively correlated with this archetype. ROMO1 is known to increase ROS production ³⁰ and thus might be downregulated by archetype 1 as an adaptation to compensate for the high ROS production generated by energy production via the electron transport chain.

Another noticeable contrast between archetypes 1 and 2 was the positive correlation of the expression of CKS2 and TPT1 with archetype 2, while the expression of these two genes were negatively correlated with archetype 1. CKS2 and TPT1 have both been implicated in cancer, including as potential therapeutic targets ^{31, 32, 33}. TPT1 encodes the translationally controlled tumor protein (TCTP) and CKS2 encodes the cyclin-dependent kinase regulatory subunit 2, which is known to bind to and be necessary for the activity of the cyclin B1-CDK1 protein kinase, an essential factor for cells to progress past the G2 phase of the cell cycle ³⁴. Interestingly, both TPT1 and CKS2 have also been associated with OXPHOS-induced oxidative stress ^{35, 36}, as well as with cell-cycle regulation and with protein synthesis/degradation ^{37, 38}. In an hematopoietic cells (HSC) mouse model, CKS2 knockout was associated with an accelerated cell cycle ³⁹ (which may contribute to the malignant phenotype) but also with an increase ROS production ³⁸. Moreover, the same study also found that CKS2 was involved in proteostasis of HSC, which may be related to the trade-off in protein synthesis identified for archetype 1 ³⁸. The role of TCTP in protein synthesis regulation is known to involve its interaction with elongation factors eEF1A and eEF1B ³⁷. The fact that both eEF1A1 and eEF1B2 followed the same correlation of expression patterns between archetypes 1 and 2 as TPT1 and CKS2 supports this interpretation.

In general, TPT1/TCTP is thought to protect cells against apoptosis and oxidative stress ⁴⁰. It might thus seem surprising that its expression should be negatively correlated with archetype 1 given the proposed trade-offs involving increased OXPHOS-induced ROS and ISR activation. However, while mild oxidative stress was found to upregulate TCTP, strong oxidative stress was found to downregulate its expression ³⁶, and both CKS2 and TPT1 were downregulated in a butyrate resistant cell line conferring tumorigenesis, apoptosis and stress resistance in a colon adenocarcinoma model ⁴¹. Thus, the level of oxidative stress in archetype 1 might reach levels associated with reduced TCTP expression.

TCTP and CKS2 expression patterns are linked through the TCTP/CDC25C/CDK1 pathway, which was shown to be dysregulated in hepatocellular carcinoma ⁴². Overexpression of TPT1 has been associated with reduced CDK1 activity via ubiquitin-proteasome degradation of CDC25C, which is necessary for the dephosphorylation and activation of CDK1 ⁴². By contrast, CKS2 is thought to promote CDK1 expression ³¹ and to be required for its function. Thus, the reciprocal correlation patterns seen between archetype 1 and 2 might reflect the need to compensate for the reduced activation of CDK1 by CDC25C from the increased expression of TPT1 by an increase in CKS2 expression and vice versa.

In addition, TCTP and CKS2 were both found to exhibit reciprocal repression with p53. In the case of TPT1, p53 is repressed via TCTP ubiquitin-mediated degradation of p53 while p53 directly represses TPT1 transcription ^{33, 43}. Likewise, CKS2 expression was found to be repressed by p53 ⁴⁴, while the overexpression of CKS2 was associated with reduced p53 protein abundance in gastric cancer ⁴⁵. However, the mechanism by which TPT1 and CKS2/CDK1 are repressed by p53 has been questioned and might be due to the indirect DREAM pathway rather than direct interaction with p53 ⁴⁶. Interestingly, CDC25C was also found to be repressed by p53 ⁴⁷ and thus complex TPT1/CDC25C/CKS2/CDK1/p53 interactions might be behind the TPT1-CKS2 opposite correlation pattern seen in these two archetypes. A schematic and synthetic representation of the hypothetical model outlined above for the metabolic trade-offs suggested by the analysis of archetypes 1 and 2’s defining genes is provided on Fig. 3.

A breakdown of the genes characterizing archetype 3 and 4 was also conducted (Fig. 3c),d)). Many of these genes turned out to be known for their involvement in cancer in general and/or lymphoma in particular. For instance, LMO2 expression was negatively correlated with archetype 4, while IGHM expression was positively correlated with this archetype. This expression pattern has previously been associated with the activated B-cell (ABC) DLBCL subtype ⁴⁸. LMO2 expression reduces double-strand break DNA repair mechanisms and has been associated with a better prognosis in DLBCL patients treated with poly(ADP-ribose) polymerase (PARP) inhibitors ⁴⁹. Interestingly, the expression of HLA-A,B,C,E, B2M and HLA-DRA/DRB1 was positively correlated with archetype 3, while expression of immunoglobulins (Ig) IGHM, IGHV4-34, IGHV5-51, IGKV3-20, IGLC2, IGLV1-47, IGLV3-1, IGLV3-19, IGLV3-21, JCHAIN were negatively correlated. In contrast, the expression of several Ig genes was positively correlated with archetype 4, while only the invariant HLA-DRA was positively correlated with the archetype. This pattern would be consistent with an immune system escape from MHC loss ⁵⁰ via a partial plasmablast cell differentiation pathway ⁵¹ in archetype 4. Such a strategy would also be consistent with the many pro-tumorigenesis effect of cancer derived Ig that have been identified, including proliferation, migration, invasion, survival, and immune evasion through inhibitory effect on antibody-dependent cell-cytotoxicity (ADCC) from NK cells ⁵².

However, since Ig are common neo-antigens in B-cell malignancies and Ig-derived neoantigen presentation by MHC is a general phenomenon in lymphomas including DLBCL ⁵³, archetype 3 may limit the production of Ig in the context of retained MHC expression to avoid displaying neo-antigen Ig to the immune system ⁵⁴. Moreover, archetype 3 was also positively correlated with the expression of CCL18, CCL19, CXCL9 and CXCL10. CCL18 is known to increase the proliferation of B-cell lymphoma ⁵⁵ and may contribute to immune evasion from its effect on immune surveillance mediated by macrophages and DC and by simultaneously favoring T cell-tolerance ^{55, 56, 57}. CCL19 directs B-cell migration after activation via antigen binding and is known to be upregulated in both gcb and abc DLBCL subtypes. Recently, autocrine CCR7-CCL19 signaling was proposed to significantly contribute to lymphomagenesis under malignant conditions via a stronger activation of the survival pathways ⁵⁸. More generally, CCR7 signalling upon binding to its ligands (CCL19/21) is associated with many pro-tumorigenic effects in hematological malignancies, including migration, proliferation, survival and immune evasion ⁵⁹. Interestingly, increased expression of CCL19, CXCL9 and CXCL10 is associated with recruitment of immature CD56^bright NK cells with low perforin content in the tumor microenvironment (TME), which is thought to protect tumor cells from NK cells ⁶⁰. Additionally, MHC expression is a well-known way for tumor cells to avoid immune detection and destruction by NK cells ⁶⁰. This correlation pattern involving MHC/CXCL9-10/CCL18-19 might thus represent the signature of an alternative immune evasion strategy for archetype 3, distinct form the one displayed by archetype 4. Finally, both archetypes were positively correlated with the expression of CD74, a gene involved in B-cell differentiation, proliferation and survival ⁶¹. An integrative schematic summary of the immune evasion trade-offs suggested by the above analysis of archetypes 3 and 4 is provided on Fig. 4.

By looking at some of the genes defining each archetype, trade-offs and expression patterns noted by other studies using different approaches could be detected. Besides the cases already highlighted, the proportion of archetype-defining genes found in nine comparative DLBCL or non-Hodgkin lymphoma gene expression studies ^{62, 63, 64, 65, 66, 67, 68, 69, 70} was assessed. Overall, 23% of the archetype defining genes identified by ParTI were part of the significantly differently expressed genes reported by these studies. Of note, the recently comparative transcriptomic study of Rapier-Sharman et al. ⁶⁶ compared RNASeq gene expression data from 322 samples, including 134 B-cell lymphoma samples (of which 123 were LBCL/DLBCL) and 188 healthy B-cell controls. Among the 20 most differently expressed genes between B-lymphoma and normal control samples, 7 (35%) were found among the archetype defining genes identified here by the ParTI algorithm (CXCL9, CXCL13, C1QA, C1QB, C1QC, CCL18, CCL19). An additional three (15%) archetype defining genes identified by ParTI were among the 20 genes showing the most significant differences in the presence of splice variants (APOE, COL1A1 and RPL5). Although the tasks and trade-offs identified by ParTI are not necessarily expected to match differently expressed genes between malignant and normal cells, finding a substantial overlap in genes identified by these two kinds of analyses is convincing evidence that the Pareto optimality theoretical framework is uncovering biologically meaningful and interpretable information at the scale of systems organization, with potential therapeutic relevance.

Uncovering cancer cell vulnerabilities in the form of trade-offs represents a promising avenue to avoid the emergence of treatment resistance. The presence of trade-offs implies that certain tasks cannot be completely avoided by the cells, and yet cannot be simultaneously optimized to the maximum level allowable in principle by the genetic potential. The genes underlying these trade-offs are thus attractive therapeutic targets that could make resistance more difficult to acquire for the malignant cells. The Pareto task inference approach predicts that whenever such trade-offs exist, they should produce detectable geometrical structures in the data. We further predicted that besides geometrical patterns in trait space, phenotypic optimums (archetypes) under trade-offs would be significantly enriched in biological functions, and characterized by patterns of different combinations of a certain proportion of shared tasks/genes.

The data analyzed here do indeed confirm those three predictions at a high level of statistical significance. The t-ratio test empirically calculates the probability of observing a polytope providing as good or better fit to the data as compared to the best possible fit defined as the convex hull. To this end, the data are randomized and the best fitting polytope and its ratio to the convex hull for each replicate data set are re-estimated. The value of the observed ratio is then compared to the distribution of ratios from the simulated randomized data to derive its probability. This test strongly supports the presence of a polytope defined by 4 vertices (tetrahedron) in this transcriptomic dataset.

According to the theory, the vertices of this polytope should represent optimal specialist phenotypes in trait space. Thus, these archetypes are predicted to be significantly enriched in certain particular functions, some of which being either different or performed by different genes, and others shared among different archetypes. FDR-adjusted Fisher’s exact-tests on archetype defining gene lists clearly show that these genes are non-random genomic sub-samples, being statistically highly significantly enriched in particular functions. Some of the most significantly enriched functions were different among archetypes, although the third and fourth archetypes showed substantial overlap in broadly defined functions, but little overlap in the genes underlying these functions within each archetype (the overall percentage of positively and negatively correlated genes unique to either archetype in pairwise comparison was 72.5%). The fact that different archetypes are statistically significantly enriched in different functions and/or gene combinations is inconsistent with an enrichment simply reflecting a general “B-lymphocyte” functional category. Given the typical complete effacement of lymph node architecture and extensive infiltration by malignant B-lymphocytes in DLBCL, these specialized functional sub-categories are likely to reflect, at least partly, various B-cell malignant phenotypic strategies, although some contribution from other cell types from the TME such as dendritic cells, macrophages or T-lymphocytes cannot be completely excluded.

However, besides these different functional characteristics, archetypes also displayed significant proportions of shared genes. As predicted if archetypes result from optimization in the face of trade-offs, these shared elements were distributed in different combinations and proportions among different archetypes. This pattern is consistent with archetypes having to reconcile various functional constrains by shuffling certain genetic toolkits, turning on and off the expression of genes in a way that preserves certain core functions and functional combinations, at the expense of other less essential and potentially dispensable elements. Whether the archetypes identified by the ParTI approach represent irreversible commitment to certain phenotypic pathways or phenotype through which cells can cycle sequentially remains an open question. The former possibility could thus represent stages in malignant progression, while the latter would include the possibility that different archetype might represent temporary adaptations to transient intra or extra-cellular environmental circumstances. Our analysis also confirmed that the archetypes identified by the Pareto approach are biologically interpretable and can be used to generate hypotheses about possible mechanisms underlying the identified correlation patterns. Thus, taken together, these results broadly confirm the predictions of the Pareto optimality theory as applied to these transcriptomic data.

If trade-offs can be identified, therapeutic approaches could be tailored to exploit them in order to minimize the risk of resistance. Thus, the metabolic and immune evasion trade-offs suggested by the data may represent therapeutic opportunities that deserve further study. The first trade-off supports recent findings suggesting an important role for mitochondria, OXPHOS and ROS in tumorigenesis ^{71, 72, 73}, including in relation with the stress responses induced by increased ROS production and their impact on macromolecule synthesis ⁷⁴. The second trade-off is in line with recent suggestions that in follicular lymphomas, loss of MHCII may be selectively acquired in cells that have accumulated immunogenic mutations in their idiotype sequences in order to avoid displaying Ig neoantigens to T-cells ⁷⁵. In that study, MHCII expressing lymphoma cells were associated with a TME rich in a CD4⁺ T-cell population with a high cytotoxicity expression profile signature (CD4_CTL), while the reverse was observed for cells expressing low levels of MHCII ⁷⁵. The pattern revealed by ParTI is also consistent with earlier findings that MHC loss can occur through a partial plasmablastic phenotypic differentiation which could be associated with high levels of Ig production ⁵¹. As mentioned previously, the retention of MHC associated with the expression of CXCL9/10 and CCL18/19 in archetype 3 might represent an alternative immune evasion strategy to avoid NK-cells detection (from the expression of MHC and the action of CXCL9/10 and CCL19) and promote T-cell tolerance from CCL18 ^{55, 56, 76}. Uncovering the best way to exploit such trade-offs in energy production, protein synthesis and immune evasion, will require additional detailed studies focused at testing the effect of disrupting the function of specific archetype defining genes on both side of the trade-offs simultaneously.

Acknowledgements

We thank Uri Alon and his research group for useful advices on running the PartiCode package.

Author contributions

JB carried out the Pareto task inference analysis, the statistical functional enrichment analysis, results interpretation and drafted the manuscript. JJ provided the bioinformatic pipeline for data formatting, PCAtest analysis, circos plots and associated functional annotations, and contributed to the writing of the final manuscript.

Arnold M, et al. Progress in cancer survival, mortality, and incidence in seven high-income countries 1995-2014 (ICBP SURVMARK-2): a population-based study. Lancet Oncol 20, 1493-1505 (2019).
Vasan N, Baselga J, Hyman DM. A view on drug resistance in cancer. Nature 575, 299-309 (2019).
Paucek RD, Baltimore D, Li G. The Cellular Immunotherapy Revolution: Arming the Immune System for Precision Therapy. Trends Immunol 40, 292-309 (2019).
Sonnenschein C, Soto AM. Over a century of cancer research: Inconvenient truths and promising leads. PLoS Biol 18, e3000670 (2020).
Selvarajoo K, Giuliani A. Systems Biology and Omics Approaches for Complex Human Diseases. Biomolecules 13, 1080 (2023).
Monti N, Verna R, Piombarolo A, Querqui A, Bizzarri M, Fedeli V. Paradoxical Behavior of Oncogenes Undermines the Somatic Mutation Theory. Biomolecules 12, (2022).
Bukowski K, Kciuk M, Kontek R. Mechanisms of Multidrug Resistance in Cancer Chemotherapy. International Journal of Molecular Sciences 21, 3233 (2020).
Saleh R, Elkord E. Acquired resistance to cancer immunotherapy: Role of tumor-mediated immunosuppression. Seminars in Cancer Biology 65, 13-27 (2020).
Salehi S, et al. Clonal fitness inferred from time-series modelling of single-cell cancer genomes. Nature 595, 585-590 (2021).
Hart Y, et al. Inferring biological tasks using Pareto analysis of high-dimensional data. Nat Methods 12, 233-235, 233 p following 235 (2015).
Hausser J, et al. Tumor diversity and the trade-off between universal cancer tasks. Nat Commun 10, 5423 (2019).
Hausser J, Alon U. Tumour heterogeneity and the evolutionary trade-offs of cancer. Nat Rev Cancer 20, 247-257 (2020).
Sun M, Zhang J. Rampant False Detection of Adaptive Phenotypic Optimization by ParTI-Based Pareto Front Inference. Molecular Biology and Evolution 38, (2020).
Mikami T, Iwasaki W. The flipping t ‐ratio test: Phylogenetically informed assessment of the Pareto theory for phenotypic evolution. Methods in Ecology and Evolution 12, (2021).
Adler M, et al. Controls for Phylogeny and Robust Analysis in Pareto Task Inference. Mol Biol Evol 39, (2022).
Zhao S, Ye Z, Stanton R. Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols. Rna 26, 903-909 (2020).
Vieira V. Permutation tests to estimate significances on Principal Components Analysis. Computational Ecology and Software 2, 103-123 (2012).
Camargo A. PCAtest: testing the statistical significance of Principal Component Analysis in R. PeerJ 10, e12967 (2022).
Chan TH, Liou JY, Ambikapathi A, Ma WK, Chi CY. Fast algorithms for robust hyperspectral endmember extraction based on worst-case simplex volume maximization. In: 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)) (2012).
Ashburner M, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25, 25-29 (2000).
Aleksander SA, et al. The Gene Ontology knowledgebase in 2023. Genetics 224, (2023).
Thomas PD, Ebert D, Muruganujan A, Mushayahama T, Albou LP, Mi H. PANTHER: Making genome-scale phylogenetics accessible to all. Protein Sci 31, 8-22 (2022).
Krzywinski MI, et al. Circos: An information aesthetic for comparative genomics. Genome Research, (2009).
Park JS, Kang DH, Bae SH. p62 prevents carbonyl cyanide m-chlorophenyl hydrazine (CCCP)-induced apoptotic cell death by activating Nrf2. Biochem Biophys Res Commun 464, 1139-1144 (2015).
Ruan Y, et al. CHCHD2 and CHCHD10 regulate mitochondrial dynamics and integrated stress response. Cell Death Dis 13, 156 (2022).
Bilen M, Benhammouda S, Slack RS, Germain M. The integrated stress response as a key pathway downstream of mitochondrial dysfunction. Current Opinion in Physiology 27, 100555 (2022).
Zhang G, et al. Integrated Stress Response Couples Mitochondrial Protein Translation With Oxidative Stress Control. Circulation 144, 1500-1515 (2021).
Jiang T, Wang Y, Wang X, Xu J. CHCHD2 and CHCHD10: Future therapeutic targets in cognitive disorder and motor neuron disorder. Frontiers in Neuroscience 16, (2022).
Gundamaraju R, Lu W, Manikam R. CHCHD2: The Power House's Potential Prognostic Factor for Cancer? Front Cell Dev Biol 8, 620816 (2020).
Amini MA, Talebi SS, Karimi J. Reactive Oxygen Species Modulator 1 (ROMO1), a New Potential Target for Cancer Diagnosis and Treatment. Chonnam Med J 55, 136-143 (2019).
You H, Lin H, Zhang Z. CKS2 in human cancers: Clinical roles and current perspectives (Review). Mol Clin Oncol 3, 459-463 (2015).
Lee H-J, et al. Targeting TCTP sensitizes tumor to T cell-mediated therapy by reversing immune-refractory phenotypes. Nature Communications 13, 2127 (2022).
Acunzo J, Baylot V, So A, Rocchi P. TCTP as therapeutic target in cancers. Cancer Treat Rev 40, 760-769 (2014).
Wang Z, Zhang M, Wu Y, Yu Y, Zheng Q, Li J. CKS2 Overexpression Correlates with Prognosis and Immune Cell Infiltration in Lung Adenocarcinoma: A Comprehensive Study based on Bioinformatics and Experiments. J Cancer 12, 6964-6978 (2021).
Jonsson M, Fjeldbo CS, Holm R, Stokke T, Kristensen GB, Lyng H. Mitochondrial Function of CKS2 Oncoprotein Links Oxidative Phosphorylation with Cell Division in Chemoradioresistant Cervical Cancer. Neoplasia 21, 353-362 (2019).
Lucibello M, et al. TCTP is a critical survival factor that protects cancer cells from oxidative stress-induced cell-death. Exp Cell Res 317, 2479-2489 (2011).
Bommer UA, Telerman A. Dysregulation of TCTP in Biological Processes and Diseases. Cells 9, (2020).
Grey W, et al. The CKS1/CKS2 Proteostasis Axis Is Crucial to Maintain Hematopoietic Stem Cell Function. Hemasphere 7, e853 (2023).
Grey W, et al. The Cks1/Cks2 axis fine-tunes Mll1 expression and is crucial for MLL-rearranged leukaemia cell viability. Biochim Biophys Acta Mol Cell Res 1865, 105-116 (2018).
Bommer UA. The Translational Controlled Tumour Protein TCTP: Biological Functions and Regulation. Results Probl Cell Differ 64, 69-126 (2017).
López de Silanes I, et al. Acquisition of resistance to butyrate enhances survival after stress and induces malignancy of human colon carcinoma cells. Cancer Res 64, 4593-4600 (2004).
Chan TH, Chen L, Guan XY. Role of translationally controlled tumor protein in cancer progression. Biochem Res Int 2012, 369384 (2012).
Amson R, et al. Reciprocal repression between P53 and TCTP. Nature Medicine 18, 91-99 (2012).
Rother K, et al. Gene expression of cyclin-dependent kinase subunit Cks2 is repressed by the tumor suppressor p53 but not by the related proteins p63 or p73. FEBS Lett 581, 1166-1172 (2007).
Kang MA, et al. Upregulation of the cycline kinase subunit CKS2 increases cell proliferation rate in gastric cancer. Journal of Cancer Research and Clinical Oncology 135, 761-769 (2009).
Fischer M, Steiner L, Engeland K. The transcription factor p53: not a repressor, solely an activator. Cell Cycle 13, 3037-3058 (2014).
Liu K, et al. The role of CDC25C in cell cycle regulation and clinical cancer therapy: a systematic review. Cancer Cell International 20, 213 (2020).
Blenk S, et al. Germinal center B cell-like (GCB) and activated B cell-like (ABC) type of diffuse large B cell lymphoma (DLBCL): analysis of molecular predictors, signatures, cell cycle state and patient survival. Cancer Inform 3, 399-420 (2007).
Parvin S, et al. LMO2 Confers Synthetic Lethality to PARP Inhibition in DLBCL. Cancer Cell 36, 237-249.e236 (2019).
de Charette M, Houot R. Hide or defend, the two strategies of lymphoma immune evasion: potential implications for immunotherapy. Haematologica 103, 1256-1268 (2018).
Wilkinson ST, et al. Partial plasma cell differentiation as a mechanism of lost major histocompatibility complex class II expression in diffuse large B-cell lymphoma. Blood 119, 1459-1467 (2012).
Cui M, Huang J, Zhang S, Liu Q, Liao Q, Qiu X. Immunoglobulin Expression in Cancer Cells and Its Critical Roles in Tumorigenesis. Front Immunol 12, 613530 (2021).
Khodadoust MS, et al. B-cell lymphomas present immunoglobulin neoantigens. Blood 133, 878-881 (2019).
Han G, et al. Follicular Lymphoma Microenvironment Characteristics Associated with Tumor Cell Mutations and MHC Class II Expression. Blood Cancer Discov 3, 428-443 (2022).
Korbecki J, Olbromski M, Dzięgiel P. CCL18 in the Progression of Cancer. International Journal of Molecular Sciences 21, 7955 (2020).
Cardoso AP, et al. The immunosuppressive and pro-tumor functions of CCL18 at the tumor microenvironment. Cytokine Growth Factor Rev 60, 107-119 (2021).
Kidani Y, et al. CCR8-targeted specific depletion of clonally expanded Treg cells in tumor tissues evokes potent tumor immunity with long-lasting memory. Proc Natl Acad Sci U S A 119, (2022).
Uhl B, et al. Distinct Chemokine Receptor Expression Profiles in De Novo DLBCL, Transformed Follicular Lymphoma, Richter's Trans-Formed DLBCL and Germinal Center B-Cells. Int J Mol Sci 23, (2022).
Cuesta-Mateos C, Terrón F, Herling M. CCR7 in Blood Cancers - Review of Its Pathophysiological Roles and the Potential as a Therapeutic Target. Front Oncol 11, 736758 (2021).
Castriconi R, et al. Molecular Mechanisms Directing Migration and Retention of Natural Killer Cells in Human Tissues. Frontiers in Immunology 9, (2018).
Zhao S, et al. High frequency of CD74 expression in lymphomas: implications for targeted therapy using a novel anti-CD74-drug conjugate. J Pathol Clin Res 5, 12-24 (2019).
Steen CB, et al. The landscape of tumor cell states and ecosystems in diffuse large B cell lymphoma. Cancer Cell 39, 1422-1437.e1410 (2021).
Dybkær K, et al. Diffuse large B-cell lymphoma classification system that associates normal B-cell subset phenotypes with prognosis. J Clin Oncol 33, 1379-1388 (2015).
de Groot FA, et al. Biological and Clinical Implications of Gene-Expression Profiling in Diffuse Large B-Cell Lymphoma: A Proposal for a Targeted BLYM-777 Consortium Panel as Part of a Multilayered Analytical Approach. Cancers (Basel) 14, (2022).
Monti S, et al. Molecular profiling of diffuse large B-cell lymphoma identifies robust subtypes including one characterized by host inflammatory response. Blood 105, 1851-1861 (2005).
Rapier-Sharman N, Clancy J, Pickett BE. Joint Secondary Transcriptomic Analysis of Non-Hodgkin's B-Cell Lymphomas Predicts Reliance on Pathways Associated with the Extracellular Matrix and Robust Diagnostic Biomarkers. J Bioinform Syst Biol 5, 119-135 (2022).
Davies A, et al. Gene-expression profiling of bortezomib added to standard chemoimmunotherapy for diffuse large B-cell lymphoma (REMoDL-B): an open-label, randomised, phase 3 trial. Lancet Oncol 20, 649-662 (2019).
Michaelsen TY, et al. A B-cell–associated gene signature classification of diffuse large B-cell lymphoma by NanoString technology. Blood advances 2, 1542-1546 (2018).
Kotlov N, et al. Clinical and Biological Subtypes of B-cell Lymphoma Revealed by Microenvironmental Signatures. Cancer Discov 11, 1468-1489 (2021).
Tripodo C, et al. A Spatially Resolved Dark- versus Light-Zone Microenvironment Signature Subdivides Germinal Center-Related Aggressive B Cell Lymphomas. iScience 23, 101562 (2020).
Liu Y, Shi Y. Mitochondria as a target in cancer treatment. MedComm (2020) 1, 129-139 (2020).
Ghosh P, Vidal C, Dey S, Zhang L. Mitochondria Targeting as an Effective Strategy for Cancer Therapy. Int J Mol Sci 21, (2020).
Vasan K, Werner M, Chandel NS. Mitochondrial Metabolism as a Target for Cancer Therapy. Cell Metabolism 32, 341-352 (2020).
Jin P, et al. Mitochondrial adaptation in cancer drug resistance: prevalence, mechanisms, and management. Journal of Hematology & Oncology 15, 97 (2022).
Han G, et al. Follicular Lymphoma Microenvironment Characteristics Associated with Tumor Cell Mutations and MHC Class II Expression. Blood Cancer Discovery 3, 428-443 (2022).
Seliger B, Koehl U. Underlying mechanisms of evasion from NK cells as rationale for improvement of NK cell-based immunotherapies. Frontiers in Immunology 13, (2022).

There is NO Competing Interest.

SupplementaryFigure1.pdf
Supplementary Figure 1
SupplementaryFigure2.pdf
Supplementary Figure 2
Supplementarytable1.xlsx
Supplementary table 1
Supplementarytable2.xlsx
Supplementary table 2
Supplementarytable3.xlsx
Supplementary table 3

Download PDF

Version 1

posted

You are reading this latest preprint version

Can the Pareto optimality theory reveal cellular trade-offs in diffuse large B-Cell lymphoma transcriptomic data?

Status:

Version 1

Abstract

Figures

INTRODUCTION

MATERIALS AND METHODS

RESULTS AND DISCUSSION

DECLARATIONS

REFERENCES

Additional Declarations

Supplementary Files

Status:

Version 1