Immune-focused multi-omics analysis of prostate cancer: leukocyte Ig-Like receptors are associated with disease progression

ABSTRACT

Prostate cancer (PCa) immunotherapy has shown limited efficacy so far, even in advanced-stage cancers. The success rate of PCa immunotherapy might be improved by approaches more adapted to the immunobiology of the disease. The objective of this study was to perform a multi-omics analysis to identify immune genes associated with PCa progression to better characterize PCa immunobiology and propose new immunotherapeutic targets. mRNA, miRNA, methylation, copy number aberration, and single nucleotide variant datasets from The Cancer Genome Atlas PRAD cohort were analyzed after filtering for genes associated with immunity. Sparse partial least squares-discriminant analyses were performed to identify features associated with biochemical recurrence (BCR) in each type of omics data. Selected features predicted BCR with a balanced error rate (BER) of 0.20 to 0.51 in single-omics and of 0.05 in multi-omics analyses. Amongst features associated with BCR were genes from the Immunoglobulin Ig-like Receptor (LILR) family which are immune checkpoints with immunotherapeutic potential. Using Multivariate INTegrative (MINT) analysis, the association of five LILR genes with BCR was quantified in a combination of three RNA-seq datasets and confirmed with Kaplan-Meier analysis in both these and in an independent RNA-seq dataset. Finally, immunohistochemistry showed that a high number of LILRB1 positive cells within the tumors predicted long-term adverse outcomes. Thus, tumors characterized by abnormal expression of LILR genes have an elevated risk of recurring after definitive local therapy. The immunotherapeutic potential of these regulators to stimulate the immune response against PCa should be evaluated in pre-clinical models.

KEYWORDS:

• Prostate cancer
• biochemical recurrence
• multi-omics analysis
• leukocyte Ig-Like Receptors
• immunotherapy

Introduction

Prostate cancer (PCa) immunotherapy has been mostly attempted with therapeutic anti-cancer vaccines using either dendritic cell-based, whole cell-based, or vector-based vaccines as well as with some other approaches but always with limited efficacy.^1–3 In recent years, the development of immune checkpoint inhibitors has revolutionized cancer immunotherapy. Immune checkpoints are a series of receptors/ligands that either inhibit or activate the function of immune cells.⁴ CTLA-4 and PD-1 or its ligand PD-L1 are the most well-known immune checkpoints but several others have been identified.^5,6 The inhibition of CTLA-4 and PD-1/PD-L1 has shown impressive anti-tumor activity in cancers such as melanoma, lung, kidney, and bladder cancers and efforts are being made to improve their efficacy, notably through the identification of biomarkers for the selection of patients more likely to respond or through combination with other drugs or therapies.⁷

Initial attempts of PCa immunotherapy using immune checkpoint inhibitors have been however rather disappointing. Two Phase III trials testing Ipilimumab (anti-CTLA-4) for the treatment of metastatic castrate-resistant prostate cancer (CRPC) showed no effect on overall survival compared to placebo although it had some positive impact on progression-free survival.^8,9 Phase I and II trials testing PD-1 or PD-L1 inhibitors have also been conducted. These studies showed that a higher anti-tumor activity could be observed in subsets of patients with tumors showing DNA repair abnormalities, inactivating CDK12 mutations or when the inhibitors were used in combination with other drugs such as enzalutamide or olaparib.^10–12 Although these results are encouraging, PCa immunotherapy is still suboptimal, and more effective approaches must be identified. A better understanding of the antitumor immune defects associated with PCa progression will help to develop immunotherapies more adapted to this disease.

To fill this gap, we performed a multi-omics analysis using The Cancer Genome Atlas (TCGA) PCa datasets to identify immune-related features associated with biochemical recurrence (BCR; i.e a rise in PSA level after local therapy) with the presumption that the identification of features common to different types of omics data would support their relevance and importance in the immunobiology of PCa. These analyses led us to identify some leukocyte immunoglobulin-like receptors (LILR) as candidate biomarkers of BCR.

LILR is a family of immune receptors that either activate (LILRA members) or suppress (LILRB members) immune cell functions. These Type 1 transmembrane glycoproteins are composed of two or four extracellular Ig-like domains that bind ligands and either a short cytoplasmic tail with an immunoreceptor tyrosine-based activation motif (ITAM), for the LILRA members of the family, or a long cytoplasmic tail with an immunoreceptor tyrosine-based inhibitory motif (ITIM), for the LILRB members of the family.^13,14 These receptors are widely expressed in hematopoietic-lineage cells but their exact function is still not well understood. LILRA and LILRB have been shown to bind to various ligands including membrane-bound proteins such as MHC class I molecules (most strongly with HLA-G molecules) and soluble proteins such as angiopoietin-like proteins (ANGPTLs), myelin inhibitors, and S100A8/9 proteins.¹⁴ LILR up- or downregulation was shown to impact the response to bacterial and viral infections as well as to influence the outcomes in diseases such as autoimmunity, inflammatory diseases, and cancer.¹⁵

We report here the details of our multi-omics analyses and discuss the potential of LILR and especially LILRB1 as targets for PCa immunotherapy.

Results

Eligibility and preparation of data

The TCGA PRAD project comprises 498 men treated with radical prostatectomy. From these cases a large number are inadequate in terms of quality of RNA sequencing as indicated by the TCGA PRAD team.¹⁶ As our objective was to identify features associated with BCR which may happen several years after radical prostatectomy, we selected cases with a minimum of 5 years of clinical follow-up and combined them with those cases where BCR happened before 5 years. We, therefore, discarded several cases that did not encounter our eligibility criteria (see details in Material and Methods). This imposed an important selection as only 45 cases were conserved for this analysis (Supplementary Figure S1). The resulting cohort is enriched in high-risk PCa which may help to identify features associated with BCR (Supplementary Table S1).

mRNA, miRNA, methylation, CNA, and SNV datasets from the TCGA PRAD project were downloaded and curated. Non-informative data in each dataset were discarded. RNA-seq data were completely reanalyzed. The numbers of CNA and SNV data were considerably reduced by the selection of features. Following this first step, the number of RNA, miRNA, methylation, CNA, and SNV features were 29,820, 1,211, 20,112, 13,925, and 928, respectively. We next applied our filter to select features associated with a set of 812 immune genes. After filtering for immune genes the resulting number of RNA, miRNA, methylation, CNV, and SNV features were 768, 1,211, 768, 138, and 6, respectively (Supplementary Figure S2).

Features associated with BCR

To identify immune gene-related features associated with BCR, prediction modeling was performed using sparse partial least squares-discriminant analysis (sPLS-DA). These analyses were first performed on each data type separately. Supplementary Table 2 provides the list of features identified by sPLS-DA in each omics dataset. Overall, 51 mRNA, 44 miRNA, 36 methylation loci, 32 CNA, and 6 SNV were identified. The selected mRNA, miRNA, and methylation loci predicted well occurrence of a BCR with balanced error rate (BER) of 0.20, 0.23, and 0.26, respectively, while the selected CNA and SNV, with BER of 0.46 and 0.51, respectively, did not predict as well BCR (Figure 1). We next merged all those features into one single set of data and performed a general sPLS-DA analysis. This resulted in an almost perfect prediction of BCR with a BER of 0.05 (Figure 2).

Figure 1. Results of sPLS-DA in individual omics datasets from TCGA PRAD. To identify immune-related features associated with BCR, prediction modeling was performed using sPLS-DA. Overall, 51 mRNA, 44, miRNA, 36 methylation loci, 32 CNA and 6 SNV were identified. The selected mRNA, miRNA and methylation loci predicted well BCR with BER of 0.20, 0.23 and 0.26, respectively. With BER of 0.46 and 0.51, the selected CNA and SNV, respectively, did not predict well BCR

Figure 2. Results of sPLS-DA of the combined omics datasets from TCGA PRAD. Following the sPLS-DA of individual omics datasets, we merged those features into one single set of data and performed a general sPLS-DA analysis. This resulted in an almost perfect prediction of BCR with a BER of 0.05

Analysis of mRNA selection showed that many features were associated with leukocyte activation, cell activation, regulation of catalytic activity, immune system process, intracellular signal transduction, etc. (Supplementary Figure S3). Among these were some features associated with antigen processing and presentation that were also retrieved in the other types of omics data. These comprised predominantly HLA molecules, killer-cell immunoglobulin-like receptor (KIR), and leukocyte immunoglobulin-like receptor (LILR) genes (Supplementary Table S2). Since LILR are a family of immune regulators that are showing some potential as targets for cancer immunotherapy, we further analyzed the association between 30 LILR gene-associated features and BCR^17,18 (Supplementary Table S3). As shown in Figure 3, the 30 LILR gene-associated features alone could predict BCR in sPLS-DA analysis with a BER of 0.34 suggesting a role of this family of receptors in the progression of PCa. Similar analyses were performed with features associated with the KIR genes or the HLA genes, alone or in combination with those of LILR genes. These sPLS-DA resulted in BER that were higher than 0.34 indicating that genes of the LILR family alone were the most strongly associated with BCR (results not shown).

Figure 3. Results of sPLS-DA of LILR gene-related features in the combined omics dataset of TCGA PRAD. The analysis of the 30 LILR gene-related features in the combined dataset resulted in the prediction of BCR with a BER of 0.34

To further validate the association of LILR genes with BCR and because of the paucity of PCa multiomics datasets besides TCGA, we sought to validate this association in a combination of RNA-seq datasets of PCa that would represent more than 150 patient tumors to ensure statistical power. Therefore, we selected the GSE54460 RNA-Seq dataset from Long et al.¹⁹ and an RNA-seq dataset from VPCC²⁰ with the objective to combine them with RNA-seq from the 52 selected cases of the TCGA PRAD project. We reanalyzed these new data in the same way to ensure proper assembly of the datasets (Supplementary Figure S4). Thereafter, we used the Multivariate INTegrative (MINT) approach to assess whether the LILR genes were associated with BCR in the combined RNA-seq dataset of 171 tumors (Supplementary Table S4). Figure 4 shows that five genes from the family of LILR genes (LILRB1, LILRB2, LILRB3, LILRB5, and LILRA3) were associated with BCR with a BER of 0.34 confirming an association between LILR genes and BCR.

Figure 4. Results of the MINT analysis. Using the combined TCGA-GSE54460-VPCC RNA-seq dataset, five LILR genes were found to be associated together with BCR. The five LILR genes were associated with BCR with a BER of 0.34

To further characterize the association of these LILR genes with BCR, we performed Kaplan-Meier analyses. Figure 5 shows that, as revealed by the MINT analysis, LILRB1 is a gene that is strongly associated with BCR as a high level of LILRB1 mRNA is associated with shorter BCR-free survival (Figure 5(a); log-rank p < .0001). High level of LILRB2 mRNA is also associated with shorter BCR-free survival but this association is less significant (Figure 5(b); log-rank p = .04). A high level of LILRB5 mRNA is nearly significantly associated with a shorter BCR-free survival (Figure 5(d); log-rank p = .06) while that of LILRB3 mRNA alone is not associated with BCR-free survival (Figure 5(c); log-rank p = .27). When the mRNA levels of these four genes are summed, the association with BCR-free survival is not better (Figure 5(e); log-rank p = .008) than that of LILRB1 mRNA alone indicating that LILRB1 is the main driver of the association. Moreover, removing LILRB1 from the combination greatly affects the significance of the association (not shown), supporting the importance of LILRB1. At the opposite, a high level of LILRA3 mRNA was significantly associated with better BCR-free survival (Figure 5(e); log-rank p = .003).

Figure 5. Kaplan-Meier analysis of dichotomized LILRB1 (a), LILRB2 (b), LILRB3 (c), LILRB5 (d) and LILRA3 (f) mRNA levels in the combined TCGA-GSE54460-VPCC RNA-seq dataset. High levels of LILRB mRNA have a tendency to be associated with shorter BCR-free survival but only LILRB1 and LILRB2 can predict BCR with a significant p value. The sum of the levels of these mRNA (e) was also associated with a significantly shorter BCR-free survival. At the opposite, higher level of LILRA3 mRNA (f) is associated with a higher BCR-free survival

As the combined TCGA-GSE54460-VPC cohort is composed of high-risk tumors, we next sought to determine whether the association of these genes with BCR was maintained in a cohort of intermediate-risk PCa samples. We, therefore, performed new Kaplan-Meier analyses using a sub-cohort of the CPC-GENE project²¹ comprising 144 men operated at our institution (Supplementary Table S5). Figure 6(b-d) shows that high levels of LILRB2, LILRB3, and LILRB5 mRNA were significantly associated with shorter BCR-free survival. However, LILRB1 mRNA level was not significantly associated with BCR in this cohort although a trend can be observed. LILRA3 mRNA level was not quantified in the CPC-GENE dataset; therefore, the association of this gene with BCR could not be assessed. When taken together, the sum of the levels of LILRB1, LILRB2, LILRB3, and LILRB5 mRNA was associated with BCR (Figure 6(e); log-rank p = .009).

Figure 6. Kaplan-Meier analysis of dichotomized LILRB1 (a), LILRB2 (b), LILRB3 (c), and LILRB5 (d) mRNA levels in the CPC-GENE RNA-seq dataset. High mRNA levels of LILRB2, LILRB3, LILRB5 genes, but not that of LILRB1, were significantly associated with shorter BCR-free survival. The sum of these mRNA levels (e) was also associated with a significantly shorter BCR-free survival (HR = 2.51, p = 0,009). The association of LILRA3 mRNA level with BCR could not be analyzed as there was no LILRA3 mRNA data in the CPC-GENE dataset

The absence of a statistically significant association of LILRB1 mRNA level with BCR in the intermediate-risk cohort suggests that LILRB1 gene expression could be associated with grade. Spearman correlation analysis showed indeed that LILRB1 mRNA level was correlated with grade in the TCGA-GSE54460-VPCC cohort (r_s = 0.53, p = .01; Supplementary Table S6), while the mRNA levels of LILRB2, LILRB3, LILRB5, and LILRA3 were not associated with grade in the combined cohort. To further assess the association of the LILRB1 gene with BCR, we analyzed the expression of LILRB1 protein in a series of 20 high-risk prostate tumors by immunohistochemistry. LILRB1 protein was found on immune cells scattered between tumor glands (Supplementary Figure S5 and Table S7). No tumor cells or stromal cells expressed the protein. In Kaplan-Meier analyses, a high level of LILRB1⁺ cells infiltrating the tumor was found to be associated with poor clinical outcomes such as the need for definitive androgen deprivation therapy (Figure 7(b); log-rank 0 = 0.009) and having lethal PCa defined as PCa that has already led to death or metastatic castration-resistant PCa that will eventually lead to death by PCa.

Figure 7. Kaplan-Meier curves of BCR-free (a), ADT-free (b) and lethal PCa-free (c) survival according to high (level 3) vs low (levels 1–2) levels of LILRB1+ cells in the tumor area of high-risk PCa samples. The number of LILRB1+ cells in the tumor area of 20 formalin-fixed and paraffin-embedded T2-T3 stage prostate tumor samples with long clinical follow-up was analyzed using immunohistochemistry. Expression of LILRB1 was classified as level 1, 2 and 3 corresponding to low, intermediate and high number of positive cells. Statistical significance was determined by the log-rank test

Discussion

PCa immunotherapy using the first-generation immune checkpoint inhibitors (i.e. against CTLA-4 and PD-1/PD-L1) has shown poor success in initial clinical trials. Some explanations for this might be the limited immunogenicity of PCa cells or immunosuppressive mechanisms other than those involving these major immune checkpoints.²² To explore the immunobiology of PCa, we performed a bioinformatic analysis using different types of omics data to identify immune-related features associated with the first event of PCa progression, i.e. the BCR consisting in an elevation of serum PSA after local therapy. We hypothesized that the identification of biological features selected in different types of omics would support their biological relevance and might provide candidate molecular targets for immunotherapeutic intervention.

To perform these multi-omics analyses involving large datasets we used for variable selection sPLS-DA, a multivariate exploratory approach, that provides more insight into cell biology, biological pathways, or complex traits than other commonly used approaches such as machine learning approaches.²³ We first used this approach on the TCGA PRAD data within each type of omics to select features related to BCR and then we merged the selected features and applied the same approach on all the selected features. BCR was used as the clinical outcome of PCa progression instead of late outcomes associated with aggressive cancer such as detection of metastases or PCa-specific death. A major limitation of our approach is the size of the cohort available after the selection based on quality criteria. As reported, the TCGA PRAD cohort has a very short median follow-up which limits the correlation analyses between genomic features and clinical outcomes, especially the late ones. In order to increase the accuracy of the association between features and BCR, we selected 5 years as a minimal follow-up knowing that it would considerably reduce the size of the final cohort.

The sPLS-DA led to the identification of an association between a series of genes involved in antigen presentation and regulation of immune cell activity and PCa progression. HLA antigens, KIR, and LILR receptors form a complex system of molecules involved in the recognition of self/non-self antigens which can have an impact on various immunological responses and impact, for example the outcome of viral infections, and diseases such as autoimmunity and cancer.^13,24 Amongst these, LILRs were those that had the smallest BER associated with BCR. The association was further demonstrated in the combined RNA-seq datasets of 171 tumors using the MINT method which revealed an association between the sum of LILRB1, LILRB2, LILRB3, LILRB5, and LILRA3 mRNA levels and BCR. Kaplan-Meier analyses further confirmed this association and identified LILRB1 as the gene with the strongest association with BCR. However, validation in a cohort of intermediate-risk PCa with very few high-grade tumors showed no significant association of LILRB1 mRNA level with BCR while an association was observed for LILRB2, LILRB3, and LILRB5 mRNA levels. The absence of a significant association of LILRB1 mRNA level with BCR in this cohort is concordant with the association of LILRB1 mRNA level with Gleason grade observed in the combined TCGA-GSE54460-VPCC cohort (Supplementary Table S5). This positive relationship with grade might be explained by the fact that LIlRB1 is known to be predominantly expressed in macrophages and higher levels of M2 macrophages have been shown to be positively associated with Gleason grade and worst outcome.^25,26 Immunohistochemistry analysis of 20 high-risk tumors supported this association with adverse long-term outcomes. Multi-parametric analyses would be needed to confirm whether the immune cells expressing LILRB1 are indeed macrophages.

In cancer, LILRBs and especially LILRB1 immunosuppressive activity have been shown to contribute to cancer evasion. For example, Raji cells proliferation was proportionally inhibited by increasing amounts of HLA-G aggregated on nanoparticles and this inhibition was reversed when LILRB1 expression was inhibited using small interfering RNA or antagonistic mAb demonstrating that HLA-G inhibition is depending on LILRB1 expression.¹⁴ The immunosuppressive function of the HLA-G/LILRB1 signaling pathway has led to the identification of LILRB1 and HLA-G as new immune checkpoints that are potential targets for immunotherapy.²⁷ More recently the LILRB1/MHC-I signaling pathway was identified as a second “Don’t Eat Me” signal in tumor-associated macrophages.¹⁷ Studies of the primary “Don’t Eat Me” signal, the CD47/SIRP-$\alpha$- signaling pathway, showed that inhibition of LILRB1 or MHC-I molecules potentiates the phagocytosis of tumor cells by macrophages in a manner that is independent of the inhibition of the CD47/SIRP-$\alpha$ axis.¹⁷ Further analysis of the LILRB1/MHC-1 pathway may lead to the development of therapies to help restore macrophage function in the tumor microenvironment. Such therapies could complement the CD47/SIRP-$\alpha$-based therapies that are already showing great potential in pre-clinical and early clinical studies^28,29.

In conclusion, we performed a multi-omics analysis using PCa datasets that led us to identify a series of immune features that all together were strongly associated with BCR. Further analysis of these features allowed us to identify some candidate molecular targets that could be prioritized for immunotherapeutic intervention in PCa. Our data point toward a role for LILRB molecules and especially LILRB1 and suggest that these receptors could play a role in the resistance of PCa to anti-tumor immune response. Immunotherapeutic interventions aiming at the inhibition of the LILRB1/MHC-I pathway alone or in combination with therapies targeting complementary pathways (e.g. CD47/SIRP-$\alpha$; PD-1/PD-L1, etc.) may provide a more adapted immunotherapeutic treatment to PCa immunobiology and would hopefully lead to better clinical response. Testing this approach in pre-clinical models to assess the immunotherapeutic potential of LILRB1 inhibition to stimulate an immune response against PCa may be of significant promise.

Material and methods

Patients and datasets

This study was approved by the Research Ethics Committee of the CHU de Québec-Université Laval (Project 2018–3670). Next-generation sequencing (NGS) RNA-seq, miRNA, methylation, copy number aberration (CNA), and single nucleotide variants (SNV) data from the TCGA PRAD project (498 samples) along with their associated clinical data were downloaded from the Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov/). The GSE54460 RNA-Seq and clinical data (106 samples) published by Long et al.¹⁹ were downloaded from the Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54460). The Vancouver Prostate Cancer Center (VPCC) RNA-Seq and clinical data (85 samples) were provided by C. Collins (University of British Columbia, Vancouver, BC, Canada).²⁰ A sub-cohort from the Canadian Prostate Cancer Genome Network (CPC-GENE) (n = 144) was used as a validation cohort.¹⁹ This cohort corresponds to men that were operated in our institution. Data were downloaded from https://ega-archive.org/. All patients had localized disease and were treated by radical prostatectomy. For each patient, available clinical data comprised at least the pathological characteristics of the tumor (grade and stage), PSA level at diagnosis, the occurrence of BCR, the time between radical prostatectomy and BCR, the occurrence of death, and date of death as well as the date of the last follow-up.

Eligibility criteria

Eligibility was set by criteria of minimal quality. The TCGA PRAD project comprised of 498 participants. However, according to the TCGA Research Network, 131 participants must be omitted because of excessive RNA degradation.¹⁶ The TCGA cohort is also characterized by a short follow-up. Patients with less than 60 months of follow-up were discarded. We also ignored tumors with less than 40% of the tumor cell content. Patients treated with neoadjuvant or concomitant hormonal therapy were not conserved for the study. The same selection criteria were applied to GSE54460, VPCC, and CPC-GENE cohorts except for the selection based on the percentage of tumor cell content as this information was not provided.

OMIC data processing

RNA-Seq data

The RNA-Seq data were completely re-analyzed to avoid variability in the data processing. The use of a common pipeline of analysis ensures the accuracy of integrative analyses of transcriptomic datasets. The quality of the raw FASTQ files was assessed using FastQC³⁰ (v0.11.5) and trimmed with Trimmomatic³¹ (v0.32). A threshold quality per base of 30 (based on Phred 33) and a minimal length of 40 bases was necessary otherwise the read was not conserved for analysis. The sequences were then mapped and quantified (pseudo-alignment) on GRCh38.p7 using Kallisto (v0.43.0, default parameters, provided index).³² Kallisto provides isoform counts, adjusted for the amount of bias in the experiment to ensure a coherent non-naive mapping; consequently, gene counts were computed with tximport.³³ The Ensembl Gene ID was converted with Biomart tools^34,35 from transcript ID to gene ID. The RNA-seq counts were then normalized to negative control genes (housekeeping genes) using the RUVg method.^36,37 In order to perform this normalization, we selected from the literature a series of six housekeeping genes that could be candidates for control reference genes in PCa experiments.^38–41 These genes were: RRN18S, ACTB, PPIA, GAPDH, PGK1, and GUSB. The expression of these genes was tested by RT-qPCR in a series of 50 prostate tumors and they were shown to be stably expressed between tumor samples (data not shown). However, we excluded from the final list the ribosomal gene RRN18S because ribosomal RNAs were removed from our RNA-seq datasets. We also excluded PGK1 as it was shown that hypoxia in PCa alters the RNA abundance of this gene.³⁸ Therefore, we finally used GUSB, PPIA, GAPDH, and ACTB as negative control genes for the normalization of the counts. The same process was applied to GSE554460 and VPCC RNA-seq datasets. The dataset corresponding to 144 tumors from the CPC-GENE cohort was used for validation. The mRNA expression data were not processed as were the data from TCGA, GSE54460, and VPCC. FASTQ files were downloaded and directly used for statistical analyses.

miRNA data

The level 3 NGS miRNA data from the TCGA PRAD project were provided as normalized counts in reads-per-million-miRNA-mapped. miRNAs with a normalized count of zero or with no value were removed from the miRNA dataset. mirWalk 2.0,⁴² which relies on different databases (mirTarBase, mirDB, and TargetS), was used to assign genes to miRNA according to predicted target genes.

Methylation data

Genome-wide methylation data from the TCGA PRAD project were generated using the Illumina Infinium HumanMethylation27 and HumanMethylation450 BeadChip platform (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Methylation_LO_Pipeline/). The level 3 methylation data were used and Beta (ß) values selected. ß values (0 for an unmethylated allele to 1 for fully methylated allele) are the estimate of methylation level using the ratio of intensities between methylated and unmethylated alleles. Genes with no ß values were removed from the dataset.

CNA data

The level 3 CNA data from the TCGA PRAD project were preprocessed using Birdsuite²³ from the Broad Institute and the R package DNAcopy (v1.44).⁴³ The data were cleaned, normalized, segmented, and log transformed. From these data, we created an index with all unique regions found. These regions were annotated with ChipSeeker⁴⁴ which associates the closest genomic object to the region’s coordinates.

SNV data

The level 3 SNV data from the TCGA were retrieved and processed with the VCFtools suite.⁴⁵ An index with the mutation found in all patients at the base level was created and all unique mutations were kept as final variables. For each mutation, we kept the information of location (chromosome and coordinates) and mutation type.

Set of immune genes

The multivariate analyses were focused on immune genes. The set of selected immune genes for these analyses is composed of 812 Ensembl genes. This set of genes was derived from a meta-analysis that targeted the anti-genome of tumor cells in interaction with the immune system.⁴⁶

Statistical analyses

To identify features associated with BCR in every set of omics, sparse partial least square-discriminant analysis (sPLS-DA) models were calculated using the mixOmics package.²³ In each sPLS-DA analysis, the BCR was defined as the Y response. The Mfold validation strategy with a fold of 5 and 200 repetitions was used to ensure the stability of the model. To assess the predictive potential of mixed omics data, the selected omics features were merged and again sPLS-DA were calculated as above. The association between selected features and BCR was further analyzed by multivariate regression analysis in three RNA-seq datasets (i.e. TCGA, GSE54460, and VPCC) using the Multivariate INTegrative (MINT)⁴⁷ method from the mixOmics package. Again, the Mfold validation strategy with a fold of 5 and 200 repetitions was used as in the sPLS-DA analyses. Within all R (v3.3 to 3.6) analyses, the seed parameter was defined as 2543 for reproducibility. Code used for the PLS-DA model can be found here http://mixomics.org/methods/pls-da/and the one for the MINT model can be found here http://mixomics.org/mixmint/stemcells-example/.

The balanced error rate (BER = 1–0.5*(sensitivity + specificity)) measured at the centroid distance was used in sPLS-DA and MINT analyses to assess the quality of the association between the BCR and the omic features. A BER of 0 means a perfect classification while a BER over 0.5 means no association with the response variable for binary classification problems. We considered a BER<0.4 as a good score value.

To perform Kaplan-Meier survival analysis, the expression data for each mRNA of interest were optimally dichotomized using the Cutoff finder tool.⁴⁸ The method fitting Cox proportional hazard models to dichotomize variable was used to define a threshold abundance value. Then the survival (v2.41–3) and survminer (v0.4.1) packages were used to perform the survival analysis within R.

Immunohistochemistry

Analysis of LILRB1 expression was performed on prostate tumors obtained from our local biobank URO-1. This analysis was approved by the Research Ethics Committee of the CHU de Québec-Université Laval (Project 2012–1059). Briefly, 5-µm-thick sections of formalin-fixed and paraffin-embedded tumors were deparaffinized and submitted to heat-induced antigen retrieval (97°C, 20 min) in Tris/EDTA, pH 9 (Dako Code K8004: EnVision™ FLEX, High pH buffer) using a PT Link, Pre-Treatment Module for Tissue Specimens (Dako, Burlington, ON, Canada). Endogenous peroxidases were blocked by incubation in 3% peroxide solution for 10 min. Bound antibodies were revealed using the IDetect super stain HRP polymer kit (ID labs, London, Ontario, Canada) as follows. Slides were initially incubated for 10 minutes at room temperature with Super block solution to avoid nonspecific background staining and then with anti-LILRB1 rabbit monoclonal antibody (mAb)(clone EPR21007, dilution 1:500, Abcam, Toronto, ON) during 1 h at room temperature. After washes, slides were incubated for 30 min with HRP-Polymer Conjugate according to manufacturer’s recommendations. Final revelation was performed by a 5 min of incubation with DAB. Finally, slides were rinsed, counterstained with hematoxylin, dehydrated, and mounted with a coverslip using MM 24 low viscosity mounting medium (Leica Microsystems, Durham, USA). Slides were digitized using a Nanozoomer (Hamamatsu Photonics, Bidgewater NJ, USA) and visualized using the NDP.view2 software (Hamamatsu Photonics). Scoring of the relative number of positive cells in the tumor area was performed by a trained technician and was reported on a scale from 1 to 3.

Author contribution

BV, YF, AB, ML, and AD conceived the study. YF, AB, and AD supervised the study. BV performed all bioinformatics and statistical analyses. JL and PCB performed some bioinformatic analyses. OEM performed the analysis of housekeeping gene expression. OEM, XPL, HH, and VP were involved in the immunohistochemistry analyses. AB, ML, MLMM, PB, CC, YF, and AD interpreted data and provided advice. BV, AB, and ML wrote the manuscript. All authors critically reviewed the manuscript.

Acknowledgments

The authors would like to thank Fan Mo and Antonio Hurtado-Coll for the preparation of the data from VPCC.

Disclosure statement

No potential conflicts of interest were disclosed.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This research was realized with internal funds from the Laboratoire d’Uro-Oncologie Expérimentale. The production of RNA-seq data at VPCC was realized with funds from the Terry Fox Research Institute New Frontier Program Project Grant #1062.