Better outcomes in lung-avid/liver-averse metastatic PDAC

From a de-identified dataset of patients treated for PDAC at our institution with a complete set of disease-relevant computed tomography (CT) scans, we identified 35 patients who developed lung metastases but never developed evidence of liver metastases (hereafter referred to as the ‘lung cohort’); within this cohort, the shortest follow-up for patients alive at the time of data freeze was 760 days after resection and 984 days after diagnosis. We identified an additional 130 patients who developed liver metastases (referred to as the ‘liver cohort’), of which 28 also developed lung metastases. Consistent with previous reports1,2,3,4, we observed that lung cohort patients in our dataset fare significantly better by median overall survival (OS) than patients who developed liver metastases, regardless of whether they also developed lung metastases (819 (lung without liver) days versus 450 (liver without lung) or 537 (liver with lung) days; Fig. 1a). Median survival was also significantly longer for patients in the lung versus liver cohorts when limiting our analysis to patients treated by surgical resection (876 days versus 549 days, respectively; Fig. 1b). Patients with disease recurrence in sites other than liver or lung fared similarly to patients in the liver cohort (median survival, 693 days) and patients with no documented recurrence survived longer (median survival, 869 days; Fig. 1b).

Fig. 1: Survival outcomes and the primary organotropism gene set distinguish liver or lung recurrence independent of subtype.

a, Kaplan–Meier (K–M) estimates of OS of all patients with documented liver (n = 102 patients (pts.)) and/or lung recurrence (n = 28 and 34 pts.), P = 0.0005 and P = 0.0007. b, OS of patients treated by resection stratified by metastatic cohort; documented liver metastases (n = 84 pts.) or lung metastases without liver metastases (n = 30 pts.; P = 0.0002), recurrent disease at nonliver/lung (other) sites (n = 73 pts.) or no documented recurrence (n = 103 pts.; P = 0.003); K–M estimates (left), CPH single-variable modeling (right). c, CPH multivariable modeling of OS for patients treated by resection stratified by metastatic cohort; lung metastases (P = 0.005), liver metastases (P = 0.27), no documented recurrence (P = 0.81) and recurrent disease at nonliver/lung (other) sites (P = 0.19) combined with clinical covariates significant in single-variable modeling (n = 160 pts. with clinical covariate data). d, K–M estimates of days between resection and recurrence for metastatic cohorts; liver metastases (n = 83 pts.), recurrent disease at nonliver/lung (other) sites (n = 73; P ≤ 0.0001), or lung metastases (n = 29 pts.; P = 0.0005). e, CPH multivariable modeling of days between resection and recurrence, stratified by metastatic cohort; liver metastases (P = 0.0001), lung metastases (P = 0.005) and recurrent disease at nonliver/lung (other) sites (P = 0.14) combined with clinical covariates (n = 104 pts. with clinical covariate data). f, PurIST subtyping scores for primary and metastatic tumor specimens from patients in the liver (n = 85 pts.) and lung (n = 28 pts.; P = 0.025) cohorts. Black bars represent means. P value from two-tailed t-test. g, K–M estimates of OS for patients categorized by PurIST subtype; basal-like (n = 63 pts.) or classical (n = 206 pts.; P = 0.0003) and liver/lung cohorts; liver classical (n = 61 pts.) or lung classical (n = 19 pts.; P = 0.002). h, CPH multivariable modeling of OS for classical subtype lung cohort versus classical subtype liver cohort patients (n = 39 pts.; P = 0.0041) combined with clinical covariates. i, GSVA scores for the pORG (left; liver or lung (n = 76, P = 1.6 × 10−8), basal-like or classical (n = 218, P = 0.38)) and pSUB gene sets (center; liver or lung (n = 76, P = 0.22), basal-like or classical (n = 218, P = 7.1 × 10−27)) and PurIST scores (right; liver or lung (n = 76, P = 0.17), basal-like or classical (n = 218, P = 1.8 × 10−115)) calculated from primary tumors. j, GSVA scores for the pORG (left; liver or lung (n = 37, P = 0.91), basal-like or classical (n = 71, P = 0.39)) and pSUB gene sets (center; liver or lung (n = 37, P = 0.0013), basal-like or classical (n = 71, P = 1.1 × 10−8)) and PurIST scores (right; liver or lung (n = 37, P = 0.043), basal-like or classical (n = 71, P = 5.1 × 10−34)) calculated from metastatic tumors. k, GSVA scores for primaries versus metastases for pORG (top (n = 289, P = 0.91)), pSUB (middle (n = 289, P = 0.39)) and PurIST scores (bottom (n = 289, P = 0.39)). Patients who died <30 days after resection were omitted (a–e,g,h). P values between groups indicated with brackets determined by log-rank test, shaded regions represent 95% confidence intervals (CIs), and HR, P value and n are from CPH single-variable modeling (a,b,d,g). HR and associated P value for recurrence site variable was determined by CPH modeling; squares mark the HR estimates, and the horizontal bars represent the 95% CI (b,c,e–h). Patients with complete information on covariates were included in CPH multivariable analysis. Black bars represent means; P values were derived from one-way analysis of variance (ANOVA) tests and corrected with the Benjamini–Hochberg method and n indicates number of tumors (i–k). FU, follow-up; LN, lymph node; LV, lymph/vascular.

Source data

We performed multivariable analysis to account for clinical covariates that significantly correlated with survival in our dataset as single variables, including lymph/vascular invasion, grade, stage and lymph node positivity. Assignment to the lung cohort independently predicted longer survival for patients treated by resection in multivariable analysis, but assignment to other cohorts was not independently predictive of survival (Fig. 1c). Compared to patients in the liver cohort, lung cohort patients survived longer recurrence-free after resection (median 303 days versus 167 days, respectively; Fig. 1d) and survived longer overall after resection (median 784 versus 498 days, respectively; Extended Data Fig. 1a). By multivariable analysis, days from resection to recurrence for both liver and lung cohorts was significant independent of clinical covariates (Fig. 1e). Lung cohort patients generally survived longer after metastatic recurrence than the liver cohort (397 days versus 302 days, respectively, P = 0.053; Extended Data Fig. 1b) and survival after metastatic recurrence correlated with survival after resection (Extended Data Fig. 1c), but not with days from resection to recurrence (Extended Data Fig. 1d), suggesting biological differences in disease progression in the liver and lung cohorts between these two clinical time periods: before and after metastatic recurrence.

Lung cohort survival advantage independent of tumor subtype

We generated gene expression data by performing RNA-seq on histologically confirmed tumor regions inclusive of integrated stroma from formalin-fixed paraffin-embedded (FFPE) primary (n = 218) and metastatic (n = 71) PDAC tumors (Extended Data Fig. 1e), and then used PurIST12 to assign consensus subtypes of PDAC (basal-like or classical) to each tumor. We found that tumors from lung cohort patients skewed significantly more classical than liver cohort tumors (Fig. 1f); and, as others have reported, patients with classical subtype tumors survived longer and had longer times to recurrence than patients with basal-like tumors (600 versus 394 days; Fig. 1g; and 250 versus 153 days; Extended Data Fig. 1f)9,11. When restricted to only patients with classical subtype tumors, the lung cohort survived longer and had later recurrence than the liver cohort (1,681 versus 520 days; Fig. 1g; and 303 versus 167 days; Extended Data Fig. 1f). These results highlight a subtype-independent survival benefit for patients in the lung cohort relative to the liver cohort, which is also independent of clinical covariates (hazard ratio (HR) = 0.15, P = 0.0041; Fig. 1h).

Clinical comparisons reveal inflammation in lung cohort

We did not observe significant differences in sex, age, stage at diagnosis, tumor grade, lymph-vascular invasion or lymph node positivity between patients in the lung and liver cohorts (Table 1 and Extended Data Fig. 1g,h). Patients in the lung cohort were more likely to be treated by resection than patients in the liver cohort (89% versus 65%, respectively; Extended Data Fig. 1i); however, the survival advantage in the lung cohort is still evident when only comparing patients treated by resection (Fig. 1b). A small fraction of patients in this dataset were treated with standard-of-care neoadjuvant chemotherapy in both cohorts (Extended Data Fig. 1i), but neoadjuvant treatment did not influence OS (Extended Data Fig. 1j). By histopathology, significantly more lung cohort tumors had chronic inflammation and plasmacytoid inflammation (Table 1). Inflammatory scores were not different between the two cohorts when comparing only resected primary tumors (Extended Data Fig. 2a). More lung cohort metastases had tertiary lymphoid structures/lymphoid aggregates (TLSs/LAs) (Extended Data Fig. 2b). Perineural invasion, angiolymphatic invasion and desmoplasia were not significantly different in liver versus lung cohort primaries or metastases (Extended Data Fig. 2c,d).

Table 1 Patient demographics, disease characteristics and tumor specimen histology parameters for all patients in study and subsets categorized into liver and lung cohorts or high and low pORG score in primary tumorsSubtype-independent organotropism gene set predicts survival

We sought to identify gene expression in primary tumors associated with liver-avidity versus lung avidity/liver aversion without being influenced by the higher percentage of basal-like tumors in the liver cohort (Fig. 1f). We ran a two-factor analysis with DESeq2 (ref. 22) to identify differentially expressed (DE) genes in primary tumors from the liver cohort versus lung cohort (organotropism) and from the basal-like versus the classical subtype. To focus on the biology of metastatic organotropism independent from subtype12, we excluded the top DE genes for subtype from the DE genes for organotropism to generate a primary organotropism gene set termed pORG (55 upregulated genes). We also applied this process to the DE genes from basal-like versus classical, subtracting the top DE organotropism genes to generate a gene set termed pSUB (primary tumor subtype; 51 upregulated genes).

We used Gene Set Variation Analysis (GSVA)23 to generate activity scores of our primary tumor samples for both the pORG and pSUB gene sets. As expected, pORG scores for liver cohort primary tumors were significantly higher than those from the lung cohort, but pORG score did not significantly separate primary basal-like tumors from classical tumors (Fig. 1i, left). Conversely, pSUB scores were significantly higher for basal-like than classical tumors but not different between liver and lung cohort primary tumors, similar to PurIST (Fig. 1i, center and right). The pORG scores for metastatic samples did not distinguish liver cohort from lung cohort metastases (Fig. 1j, left). The pSUB score distinguished basal-like from classical metastases (Fig. 1j, center) and also distinguished metastases from the liver and lung cohorts, similar to PurIST (Fig. 1j, right). GSVA scores for all specimens showed a similar spread in scores between all primary and all metastatic tumors for pORG, pSUB or PurIST (Fig. 1k).

We found significant differences in OS between patients with tumors scoring high versus low for pORG, pSUB and PurIST (Fig. 2a), as well as significantly different recurrence-free survival (RFS; Extended Data Fig. 2e). Using the same high/low risk cutoffs for pORG, pSUB and PurIST scores determined in our dataset (Fig. 2a), pORG, pSUB and PurIST similarly predicted survival in two external datasets: OS in the pancreatic adenocarcinoma patient dataset (PAAD)24, reported by The Cancer Genome Atlas (TCGA) (cBioPortal) (Fig. 2b), and RFS in the Australian Pancreatic Cancer Genome Initiative (APGI)9, part of the International Cancer Genome Consortium (ICGC) study (Extended Data Fig. 2f). Low pORG primary tumors were more likely to be early stage and treated with neoadjuvant chemotherapy (Table 1); however, neoadjuvant treatment did not affect OS (Extended Data Fig. 1j). Multivariable analysis indicated that both pORG and pSUB predicted survival independently from other clinical covariates, but PurIST was influenced by grade (Fig. 2c).

Fig. 2: pORG predicts survival independently of clinical and genomic features.

a, K–M estimate of OS for patients with primary tumors having high or low pORG (left; high (n = 101 pts.), low (n = 107 pts.; P = 0.01)), pSUB (middle; high (n = 140 pts.), low (n = 68 pts.; P = 4.2× 10−5)) or PurIST scores (right; high (n = 126 pts.), low (n = 82 pts.; P = 0.00049)) from the OHSU dataset. High/low risk was determined by receiver operating characteristic curve (ROC) and maximum Youden’s index. b, K–M estimate of OS in TCGA pORG (left; high (n = 73 pts.), low (n = 67 pts.; P = 0.032), pSUB (middle; high (n = 99 pts.), low (n = 41 pts.; P = 0.087) or PurIST (right; high (n = 69 pts.), low (n = 71 pts.; P = 0.03)) patients with PDAC. High/low score is defined using cutoff from OHSU dataset. c, CPH multivariable modeling of OS versus primary GSVA score for pORG (top; P = 0.0062), pSUB (middle; P = 0.023) and PurIST (bottom; P = 0.37) with clinical covariates (n = 132 pts.). d, pORG, pSUB and PurIST scores of primaries and metastases (Met) from the same patient, grouped by clinically defined liver cohort (documented liver recurrence, n = 3 pairs) or lung cohort (documented lung recurrence without liver recurrence, n = 4 pairs) showing cohort mean GSVA (point) and 95% CI (error bars). e, pORG and PurIST scores for primaries (circles) and metastases (x) in liver and lung cohorts (n = 113 pts.). f, Fraction of primaries or metastases in each quadrant of the graph in e; liver (n = 85 pts., P = 0.93) or lung (n = 28 pts., P = 0.0012). P values from two-way chi-squared test between primary and metastatic specimens. g, UMAP of Werba et al.25 scRNA-seq, shaded by per-cell scores for pORG in PDAC primaries (top; n = 17 pts.) and PDAC liver metastases (bottom; n = 10 pts.). h,i, Oncoprints of the top ten altered genes and alteration types (n = 271 tumors) in the DNA dataset (h) and top (above, n = 50 pts.) and bottom quartile (below, n = 50 pts.) (i) by pORG primary GSVA score. j, pORG primary GSVA score versus TP53 (left; altered (n = 131 pts.), WT (n = 70 pts.), FDR = 9.3 × 10−13) or CDKN2A (right; altered (n = 55 pts.), WT (n = 146 pts.), FDR = 0.00052) gene alteration. P value from two-tailed t-test calculated for genes with ≥10 alterations in the dataset, corrected with the Benjamini–Hochberg method. k, CPH multivariable modeling of OS versus pORG GSVA score and genomic alterations prognostic in single-variable CPH modeling in primary tumors (left; n = 193 pts., P = 0.04) and all tumors (right; n = 251 pts., P = 0.014). l, Oncoprints of the top ten altered genes and their alteration types in primaries (left; n = 203 tumors) and metastases (right; n = 68 tumors). The log-rank test P values and n per group are indicated with brackets, shaded regions represent 95% CI, and CPH single-variable modeling HRs and associated P values are displayed on plots (a,b). Frequency is indicated at left, top bars indicate variant types by tumor, and right bars indicate variant types by gene (h,i,l). Alteration key (i). HR and associated P value for GSVA or PurIST score was determined by CPH modeling, squares represent HR estimates, and error bars represent 95% CIs (c,k). Patients who died within 30 days after resection are not shown (a,c,k).

Source data

Analysis of pORG, pSUB and PurIST scores in ten matched primary tumors and metastases (Extended Data Fig. 2g) revealed that lung metastases (n = 2 pairs) and metastases in the clinically defined lung cohort (n = 4 pairs) went from low in primaries to high in metastases, whereas liver cohort primaries and metastases stayed high (Fig. 2d and Extended Data Fig. 2h); consistent with pORG not distinguishing between liver and lung cohort metastatic samples (Fig. 1j). In unpaired primaries and metastases, the liver cohort had a similar fraction of low pORG samples in primaries and metastases, whereas lung cohort metastases show a shift to 70% high pORG (Fig. 2e,f). Although the metastatic TME may contribute to this shift in gene expression, analysis of a publicly available single-cell RNA-seq dataset25 demonstrated that both the pORG and pSUB gene sets are enriched in the epithelial cell populations from PDAC primaries and liver metastases (Fig. 2g and Extended Data Fig. 2i–k).

TP53 and CDKN2A alterations are enriched in high pORG tumors

We used a tumor-relevant, 595 gene sequencing panel to analyze DNA alterations from the same specimens used for RNA sequencing (RNA-seq) (271 specimens with DNA data; Fig. 2h) and compared and ranked gene alterations between liver and lung cohorts, and high and low pORG, pSUB and PurIST quartiles (Fig. 2i and Extended Data Fig. 3a–j). TP53 and CDKN2A altered primaries had significantly higher pORG GSVA scores (Fig. 2j), whereas KRAS, CDKN2B and SMAD4 altered primaries trended higher and GATA1 and ELF3 altered primaries trended lower in pORG score (Extended Data Fig. 4a). In metastases, TP53 altered tumors had higher pORG scores (Extended Data Fig. 4b); and MTAP, CDKN2A and CDKN2B altered tumors had higher PurIST scores (Extended Data Fig. 4c). We used Cox proportional hazards (CPH) multivariable modeling for OS against pORG score combined with gene alterations that were prognostic as single variables (TP53, CDKN2A, KMT2D, MTAP and ARID1B) and found that pORG score predicted shorter survival independent of genomic alterations in primaries and all samples (Fig. 2k), as did pSUB and PurIST scores (Extended Data Fig. 4d). We examined alteration differences between all primaries and all metastases in our dataset and in nonmatched samples we found MTAP and SMARCB1 trended toward more alterations in metastases compared to primaries, whereas KDM5C and GATA1 trended toward fewer alterations in metastases (Fig. 2l and Extended Data Fig. 4e). We found no consistent changes across nine paired primary tumors and metastases (Extended Data Fig. 4f). There was an average of 3.4 differences in genetic alterations between paired primaries and metastases from our gene panel, consistent with stochastic changes in matched samples (Extended Data Fig. 4g). DNA analysis indicated higher tumor cell content in basal-like versus classical primary tumors, consistent with another report14, and in high pORG and high pSUB primary tumors, but no significant difference between the liver and lung cohort primary tumors; or metastases in any of the groups (Extended Data Fig. 4h–i).

Distinct pathways enriched by pORG and pSUB gene sets

Gene set enrichment analysis (GSEA) analysis revealed that high pORG and liver cohort primary tumors were enriched (normalized enrichment score (NES) > 1.7, false discovery rate (FDR) < 0.05) in Hallmark pathways26 related to oncogene-mediated RS: G2M checkpoint, E2F targets, mitotic spindle, MYC targets V1, DNA repair, IFN-α response, cell metabolism and mitogenesis (Fig. 3a). We found that high pSUB and PurIST primary tumors were enriched in pathways related to glycolysis, epithelial–mesenchymal transition, apical junctions and hypoxia, whereas high PurIST was de-enriched for bile acid metabolism and pancreas β cells (Fig. 3b). Visualization of GSVA scores for these pathways supported results from GSEA and showed that the no documented recurrence clinical group skewed pORG low (Extended Data Fig. 5a–c). GSEA in metastatic-sample cohorts yielded significant differences only in high versus low pORG, with 7 upregulated pathways overlapping with the 12 found in primaries (Fig. 3c). Thus, separating primary tumor metastatic organotropism and molecular subtype using the pORG and pSUB gene sets identifies unique pathway enrichments.

Fig. 3: High pORG, liver-tropic PDAC is associated with replication stress tolerance and IFN response.

a–c, NES colored by FDR P-adjusted (FDR.q) value (FDR.q.val) (from one-way ANOVA) is shown for Hallmark GSEA pathways if any of the comparisons reached a NES > 1.7 and FDR.q < 0.05 from the cohorts indicated on each plot. a, Solid bars, top versus bottom quartile by pORG (n = 108 pts.); hatched bars, liver versus lung cohort (n = 76 pts.). b, Solid bars, top versus bottom quartile by pSUB (n = 108 pts.); hatched bars, top versus bottom quartile by PurIST (n = 108 pts.). c, Solid bars, top versus bottom quartile by pORG in metastases (mets) (n = 34 pts.). d, Mean differential (diff.) VIPER regulon activity scores colored by FDR.q.val (from one-way ANOVA) in top versus bottom quartile by pORG (solid bars, n = 108 pts.) and liver cohort versus lung cohort (hatched bars, n = 76 pts.) primary tumors for regulons related to cell cycle (left), DNA replication (center) and DNA damage repair (right). e, Example immunostaining of epithelial cells (KRT+), proliferation (Ki67+) and algorithmic detection of pRPA foci in PDAC tissue (n = 55 cores imaged in total). f, Mean RS pRPA foci in epithelial cells (left; high pORG (n = 16 pts.), low pORG (n = 18 pts.), P = 0.033) and Ki67+ proliferating epithelial cells (right; high pORG (n = 16 pts.), low pORG (n = 18 pts.), P = 0.036) in each patient determined by immunostaining a TMA with 34 primary specimens, 1–2 cores each. g, K–M estimate of OS for patients with tumors with high or low pORG GSVA scores stratified by tumors with or without a known pathologic somatic alteration (VUS were excluded) in a DDR-related gene (DDR altered high (n = 23 pts.), DDR intact high (n = 73 pts.), DDR altered low (n = 20 pts.) or DDR intact low (n = 77 pts.), P = 0.018). log-rank P value, and shaded regions represent 95% CI. h, pORG GSVA scores for primary tumors (Pri) and metastases (Met) from patients in the liver cohort categorized by a known pathologic somatic alteration (VUS were excluded) in a DDR-related gene (DDR altered (n = 13 tumors) or DDR WT (n = 74 tumors), P = 0.044). i, Pearson correlation (two-sided) of the indicated VIPER regulon scores and pORG GSVA scores (n = 218 pts.). j, Pearson correlation (two-sided) of pORG and IFN- and immune-related signature GSVA scores for primary tumors (n = 210 pts.). Two-tailed Students t-test P value; black bars represent the mean (f,h). P values from Pearson correlation and corrected with the Benjamini–Hochberg method (i,j). *FDR P adjusted < 0.05, **P adjusted < 0.01, ***P adjusted < 0.001. FDR-adjusted P values were 0.0037, 4.7 × 10−6, 3.5 × 10−8, 0.24, 1.8 × 10−9, 2.6 × 10−8, 5 × 10−11, 8 × 10−12, 7.4 × 10−39, n = 218 patients (for i) and 8.8 × 10−5, 1.8 × 10−5, 0.0047, 2 × 10−7, 4.1 × 10−11 (for j), n = 210 patients.

Source data

Cell cycle, RS and DNA repair up in high pORG, liver-avid tumors

Virtual inference of protein-activity enrichment regulon (VIPER)27,28 analysis followed by Gene Ontology network analyses identified nodes for cell cycle and DNA replication and repair enriched in both high pORG and liver cohort primary tumors (Extended Data Fig. 5d,e). Accordingly, cell cycle, DNA replication and DNA repair proteins demonstrated significantly higher activity in high pORG and liver cohort tumors (Fig. 3d). To further analyze RS, we immunostained a tissue microarray (TMA) prepared from 34 primary tumors using the same FFPE blocks from our RNA and DNA-seq data for foci of phosphorylated replication protein A (pRPA) (Fig. 3e), an indicator of single-stranded DNA exposed during RS. We found a significantly higher mean number of pRPA foci in cytokeratin-positive (KRT+) epithelial cells in high pORG versus low pORG primary tumors (Fig. 3f, left). Additionally, we found that Ki67+ proliferating tumor cells had significantly more pRPA foci in high pORG tumors (Fig. 3f, right). Similar, though not significant, trends were observed in nine liver cohort primary tumors compared to four lung cohort (Extended Data Fig. 6a). The percent of epithelial cells and of proliferating epithelial cells positive for pRPA foci were also higher in high pORG primary tumors and trended higher in liver cohort primaries (Extended Data Fig. 6b-c), and more pRPA+ cells were proliferating in high pORG tissues (Extended Data Fig. 6d). Together, these results support the hypothesis that in high pORG tumors, pRPA+ cells are a viable, expanding part of the tumor despite ongoing RS, likely due to the associated increase in DNA repair.

Low pORG tumors are less tolerant to defects in DNA repair

A report by Dreyer et al. suggested that treatment-agent efficacy may depend on both RS and DNA damage response (DDR) gene alteration status, dividing patients into four categories based on the presence or absence of those two factors16. As our data indicate that liver-avid, high pORG primary tumors are enriched for pathways associated with ongoing RS and DNA repair, we divided patients into four categories by high/low pORG score and the presence/absence of a known DDR gene alteration29. Although patients with high pORG scoring tumors fared poorly regardless of DDR gene status, patients with low pORG tumor scores survived significantly longer if their primary tumors had DDR nonsilent gene alterations, whether or not variants of unknown significance (VUS) were excluded (Fig. 3g and Extended Data Fig. 6e). Additionally, liver cohort tumors with DDR gene alterations had higher pORG scores compared to those without (Fig. 3h), suggesting that the presence of DDR gene alterations may promote mechanisms supporting tumor cell responses to RS and DNA damage to avoid mitotic catastrophe, and a lack of this response, as seen in low pORG tumors, combined with a DDR gene alteration improves patient outcome (Fig. 3g).

Suppressed tumor immunity in high pORG, liver-avid tumors

Consistent with enrichment of the Hallmark IFN-α response in high pORG samples by GSEA (Fig. 3a), VIPER scores for IFN-α/β receptor subunits activity positively correlated with pORG score (Fig. 3i). Chronic IFN signaling in cancer is reported to induce an IFN-related DNA damage resistance gene expression signature (IRDS), associated with tumor cell resistance to DNA damage30,31,32 and escape from tumor immunity33. We found a significant positive correlation between the IRDS gene signature and pORG score in primary tumors (Fig. 3j). Two genes in the IRDS gene set matched VIPER regulons (STAT1 and MX1) and these were both significantly positively correlated with pORG and trended higher in liver cohort tumors (Fig. 3i and Extended Data Fig. 6f).

Consistent with chronic IFN signaling inactivating adaptive immune cells34, we found that high pORG scores negatively correlated with B cell and T cell gene signatures, regulons and marker genes; and positively correlated with response to IFN, macrophage and neutrophil marker genes, signatures, and regulons (Fig. 3i,j and Extended Data Fig. 6g). Similarly, liver versus lung cohort tumors had a trend of lower CD20 B cell VIPER activity scores (FDR = 0.06; Extended Data Fig. 6f). We found similar results with deconvolution algorithms: notably, negative correlations between pORG score and most lymphocyte subsets, as well as endothelial cells and cancer-associated fibroblasts, and positive correlations between pORG score and immune suppressive TH2 CD4+ T cells, macrophages, plasmacytoid dendritic cells and γδ T cells (Fig. 4a and Extended Data Fig. 6h).

Fig. 4: Transcriptomic and multiplex imaging evidence of immune suppression in high pORG, liver-tropic tumors.

a, Pearson correlation (two-sided) of xCell deconvolution scores and pORG GSVA score for primary tumors (n = 204 pts.). P values from Pearson correlation and corrected with the Benjamini–Hochberg method. *FDR P adjusted < 0.05, **P adjusted < 0.01, ***P adjusted < 0.001. FDR-corrected P values are 0.198, 9.06 × 10−6, 0.0752, 0.802, 0.123, 2.31 × 10−9, 5.23 × 10−10, 1.26 × 10−16, 2.58 × 10−22, 0.000701, 0.0112, 0.164, 0.0215, 0.0966, 0.0125, 0.000937, 0.0682, 0.281, 0.326, 1.93 × 10−6, 3.54 × 10−5, 0.0215, 0.362, 0.228, 0.462, 0.326, 0.0147, 0.00857, 0.362, 0.018, 0.422 and 0.227. b, Representative images of mIHC staining of a low pORG, lung cohort patient tumor (left) and a high pORG, liver cohort patient tumor (right). n = 12 tissues imaged, 174 ROIs total. c, Average leukocyte densities for primary tumors from patients in the liver cohort (mean pORG 0.23 s.e.m. = 0.11, n = 9 pts.) and lung cohort (mean pORG −0.51 s.e.m. = 0.09, n = 3 pts.) (top). Average leukocyte densities for primary tumors from patients with high pORG (pORG 0.38 s.e.m. = 0.04, n = 7 pts.) and low pORG GSVA scores (mean pORG –0.43 s.e.m. = 0.08, n = 5 pts.) (bottom). DC, dendritic cell. d, Leukocyte densities in ROIs from liver (n = 121 ROIs) or lung cohort (n = 53 ROIs) primaries (top) and high (n = 84 ROIs) or low (n = 90 ROIs) pORG primaries (bottom). Each dot represents an ROI colored by patient specimen (n = 12 patients). Box represents the median and interquartile range (IQR), and whiskers extend 1.5 × IQR. P values from two-tailed t-test corrected with the Benjamini–Hochberg method. FDR-corrected P values are 0.45, 1.1 × 10−8, 0.45, 0.23, 0.001, 0.07, 0.00034, 7.6 × 10−13, 0.98 (top) and 0.0097, 0.00061, 0.049, 0.33, 0.006, 0.31, 5.4 × 10−6, 2 × 10−5, 0.034 (bottom, n = 174 ROIs). e, pORG score from RNA-seq of liver versus lung cohort tumors in the TCRβ dataset, primaries (left; liver (n = 38 pts.), lung (n = 11 pts.), P = 1.7 × 10−5) and metastases (right: liver (n = 20 pts.), lung (n = 5 pts.), P = 0.47). P values from two-tailed t-test. Black bars represent the means. f, K–M estimation of OS of patients with high (n = 106 pts.) versus low (n = 82 pts.; P = 0.0054) pORG GSVA scores (left; cutoff determined by ROC and maximum Youden’s index in the full dataset in Fig. 2a) and liver (n = 76 pts.) versus lung cohort (n = 16 pts.; P = 0.097) patients (right) in the TCRβ dataset. log-rank test P values and n patients per group are indicated with brackets and shaded regions represent 95% CI. CPH single-variable modeling HR and associated P values are displayed on plots.

Source data

We used a multiplexed immunohistochemistry (mIHC) platform35,36 to measure densities of leukocyte subsets in multiple 1.0-mm2 regions of interest (ROIs) in tissue sections from primary tumor specimens with pORG scores assigned from gene expression data and classified as