Reduced immunogenicity of MYC amplified, metastatic prostate cancer

MYC CNV in primary and metastatic datasets

To assess potential survival probability distinctions based on the amplification of oncogenes, we first assessed the CNV for several oncogenes using the precision-guided method of Mauro, Varkhedi and colleagues [11, 12], whereby we established ratios of tumor and blood read counts for each case ID in the TCGA-PRAD and WCDT-MCRPC datasets (Methods). For the TCGA-PRAD dataset, the case IDs representing the top 20 read count ratios versus all remaining cases (representing all of the lower read count ratios) were assessed with a KM analysis for progression-free survival (PFS). In the case of MYC CNV, a survival difference at the 70-month timepoint was observed, where the cases with a greater number of MYC copies represented a lower PFS probability (Figure 1A, two-proportion test p-value = 0.0113). When comparing the cases representing the top 20 and bottom 20 read count ratios with a KM analysis, a similar result was seen at the 70-month timepoint, where the cases with a greater number of MYC copies represented a worse PFS (Figure 1B, two-proportion test p-value = 0.0057). We next assessed whether the apparent MYC copy numbers (CNs) correlated with MYC gene expression levels, for the TCGA-PRAD dataset, which revealed a statistically significant positive correlation for the MYC CNs and the MYC RNAseq values (Figure 1C, Pearson’s correlation coefficient, R = 0.128, Pearson’s correlation p-value = 0.005). The assessment of a correlation of MYC copies and MYC RNAseq values was also performed for the WCDT-MCRPC dataset (Methods), which revealed a significant positive correlation for the MYC CNs and the MYC RNAseq values (Figure 2, Pearson’s correlation coefficient, R = 0.442, Pearson’s correlation p-value < 0.001).

Figure 1. Kaplan-Meier (KM) survival and gene expression analyses of TCGA case IDs representing MYC amplification (A) Progression-Free Survival (PFS) of case IDs representing WXS-based CNs of MYC, comparing top 20 CNV ratios (black line, n = 20) and all other CNV ratios (grey line, arrowheads, n = 455), with statistical analysis for the KM plot and the 70-month survival timepoint (log-rank p-value = 0.15; two-proportion time, single time point test p-value = 0.0113) (Supplementary Table 4). (B) PFS of case IDs representing WXS-based CNs of MYC, comparing top 20 CNV ratios (blue line, n = 20) and bottom 20 CNV ratios (gray line, n = 20), with statistical analysis for the KM plot and the 70-month survival timepoint (log-rank p-value = 0.19; two-proportion time, single time point test p-values = 0.0057) (Supplementary Table 5), (C) Correlation of MYC tumor:blood sequencing read ratios (CNV) with the RNAseq values for MYC (r = 0.128, p = 0.005) (Supplementary Table 6).

Figure 2. Correlation of RNAseq values with tumor: blood sequencing read ratios for MYC, representing case IDs from the WCDT-MCRPC dataset (r = 0.442, p < 0.001) (Supplementary Table 6).

MYC CN ranges compared across different datasets, to assess potential differences in amplification levels across different stages of prostate cancer

The WCDT-MCRPC dataset, representing metastatic prostate cancer, exhibited a significantly higher average MYC copy number compared to the TCGA-PRAD dataset, which represents primary prostate cancer (Table 1, p < 0.001). Similarly, the CMI-MPC dataset (Methods), also representing metastatic prostate cancer, demonstrated a notably elevated average MYC copy number when contrasted with the TCGA-PRAD dataset (Table 1, p < 0.001).

Recovery of adaptive IR recombination reads from MYC amplified and non-amplified cases

To determine whether there was a difference in the recovery of adaptive IR recombination reads for MYC-amplified versus nonamplified prostate cancer, the immunoglobulin heavy (IGH) and light (IGK, IGL) chain recombination reads from exome and RNAseq files were extracted from the TCGA-PRAD dataset (Methods) [13–15]. For these preliminary considerations, average recombination read counts for the MYC-amplified cases were not significantly different than for MYC-nonamplified cases (data not shown). This process was repeated for the WCDT-MCRPC dataset, for both WGS and RNAseq files. A cutoff tumor:blood sequencing read ratio of ≥5.0 was used to distinguish MYC amplified from non-amplified cases. For the RNAseq file-derived recombination reads, MYC amplified cases represented significantly reduced number of IGH recombination reads (p = 0.0056), IGK recombination reads (p = 0.0149), IGL recombination reads (p = 0.0039), and IGH+IGK+IGL recombination reads combined (p = 0.0021) (Table 2). For WGS file-derived recombination reads, cases with MYC amplification represented significantly reduced numbers of IGK recombination reads (p = 0.0238), IGL recombination reads (p = 0.0011), and IGH + IGK + IGL recombination reads combined (p = 0.0024) (Table 2).

Next, the T-cell receptor alpha (TRA), beta (TRB), gamma (TRG), and delta (TRD) recombination read counts from WGS and RNAseq files were extracted from the TCGA-PRAD dataset. For all these comparisons, average recombination read counts for the MYC-amplified cases were not significantly different than for MYC-nonamplified cases (data not shown). This process was repeated for the WCDT-MCRPC dataset (Methods), for both WGS and RNAseq files. A cutoff tumor:blood sequencing read ratio of ≥5.0 was used to distinguish MYC amplified from non-amplified cases, as in Table 2. For WGS file-derived recombination reads, cases with MYC amplification revealed significantly reduced numbers of TRA recombination reads (p = 0.0077), TRB recombination reads (p = 0.0167), and TRA + TRB recombination reads combined (p = 0.0073) (Table 3). No significant difference was detected in recombination read counts for RNAseq files (data not shown). For all comparisons, there was no significant difference detected in TRG or TRD recombination read counts.

Expression of an immune signature gene set in MYC amplified and non-amplified cases

To support the above observations related to the adaptive IR recombinations, the expression of an immune marker gene set [10] was assessed for MYC amplified and non-amplified cases. The RNAseq values for the immune marker gene set for the WCDT-MCRPC MYC amplified and non-amplified cases were extracted (Methods). Both “tpm_unstranded” and “fpkm_unstranded” values were extracted. This approach revealed statistically significantly lower levels of expression of many of the immune marker genes among the MYC amplified cases (Figure 3A–3D and Tables 4, 5).

Figure 3. Box and whisker plots of selected immune marker genes for comparison of expression levels (tpm and fpkm) for MYC amplified (shaded) and non-amplified WCDT-MCRPC cases (A) tpm: CD19, p-value = 0.004; CD22, p-value <0.001; CD33, p-value = 0.004; CD68, p-value = 0.006; CD79B, p-value=0.007; MS4A1, p-value = 0.002; SPHK1, p-value = 0.037. (B) tpm: CD4, p-value = 0.038; CD79A, p-value = 0.005; CIITA, p-value = 0.017; CTSC, p-value = 0.018; ITGAX, p-value = 0.010. (C) fpkm: CD19, p-value = 0.002; CD22, p-value <0.001; CD33, p-value = 0.002; CD68, p-value = 0.004; CD79B, p-value = 0.004; MS4A1, p-value = 0.002; SPHK1, p-value = 0.021. (D) fpkm: CD4, p-value = 0.017; CD79A, p-value = 0.007; CIITA, p-value = 0.023; CTSC, p-value = 0.012; ITGAX, p-value = 0.013. For the full range of RNA-seq data corresponding to the specified immune marker genes, along with data for other immune marker genes refer to Supplementary Table 7.

Access to genomics files and CNV analyses

Precision-guided copy number variation (CNV) assessments were conducted using the approach outlined by Mauro et al. [22] and further developed by Varkhedi et al. [12]. The cancer genome atlas, prostate adenocarcinoma (TCGA-PRAD, phs000178) dataset, accessed via National Institutes of Health database of genotypes and phenotypes (dbGaP) project approval number 6300, was downloaded to USF research computing. The genomic data commons (GDC, https://gdc.cancer.gov/) download tool was used to facilitate the acquisition of whole exome sequence (WXS), binary alignment map (BAM) format file slices, for both blood and tumor samples, with the file slices representing specific genes. The genes and BAM file slices were delineated by their start and end nucleotide numbers, according to the hg38 version of the human reference genome sequence available at genome.ucsc.edu. A GDC-provided reference file (sample sheet) was used to identify the case ID associated with each WXS BAM file slice, as well as to determine whether the file represented a blood-derived normal or a tumor sample. Utilizing this sample sheet, the files were programmatically renamed to reflect their corresponding case IDs and to indicate whether the files were derived from blood or tumor samples. The number of mapped reads in each file was then extracted using a function from the SAMtools package [23, 24]. Following this, a ratio of tumor to blood read counts was computed from the respective files. These ratios were subsequently compiled and outputted to a csv file for further analysis (Supplementary Tables 1–3). This approach was also applied to the count me in-metastatic prostate cancer (CMI-MPC, phs001939) and West Coast dream team-metastatic castration resistant prostate cancer (WCDT-MCRPC, phs001648) datasets, via dbGaP project approval numbers 25670 and 31203, respectively.

Sample numbers

TCGA-PRAD dataset, comprising 499 primary prostate tumor samples, was analyzed as described above. The approaches above were also applied to the CMI-MPC dataset, consisting of 63 metastatic prostate cancer samples, and the WCDT-MCRPC dataset, which included 101 metastatic castration-resistant prostate cancer samples.

Extraction of the adaptive immune receptor (IR) recombination reads from genomics files

The procedure for extracting the IR recombination reads is described in ref. [15] and the latest iteration of the software is freely available at https://github.com/kcios/2021. The entire set of data representing the WCDT-MCRPC dataset, except for the recombination sequencing reads themselves, can be obtained at the following link: https://usf.box.com/s/waqy80pz663wd5lo16midh3qdf4gv11n. (The sequencing reads are controlled access data and can only be accessed by dbGaP approved users.) For the RNAseq-based adaptive IR recombination reads, only the files with the “genomic” suffix were sourced. For the WGS-based adaptive IR recombination reads, only the metastatic tumor files were sourced.

KM analyses

The data for the survival analyses were obtained from cbioportal.org [25, 26] for TCGA-PRAD dataset (Pancancer). The survival distinctions were assessed using a web tool at cbioportal.org and then verified using R software (version 4.3.2) and the survminer package (Supplementary Tables 4 and 5).

RNAseq values

RNAseq values for TCGA-PRAD dataset (Pancancer) were obtained from cbioportal.org [16, 17]. RNAseq files for the WCDT-MCRPC dataset were obtained from the GDC. Specific data (tpm_unstranded and fpkm_unstranded values) were extracted using R software (version 4.3.2) (Supplementary Tables 6 and 7).