We performed bulk RNA-seq extracted from FFPE samples of 16 patients who also had bulk RNA-seq performed on a different section from the same tumor and deposited in TCGA. The pipeline and workflow details for the TCGA can be found on their website. (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). The second sample for these16 patients was processed using the NovaSeq S4 PE100 with Illumina’s TruSeq Total Stranded RNA prep reagents (https://www.illumina.com/products/by-type/sequencing-kits/library-prep-kits/truseq-stranded-total-rna.html) and used a second DNase treatment to minimize any potential DNA contamination. Table 1 shows select patient characteristics for our study (Table 1). The Institutional Review Board of Mayo Clinic gave approval for this work.

For alignment and quantification of the Mayo RNA-Seq data, we used a pipeline based on HiSat2 and subread (CHURP (Baller et al., 2019), HiSat2 (Kim et al., 2019), and featureCounts (Liao et al., 2014)). The TCGA and Mayo data sets were analyzed independently for deconvolution. For each data set, raw counts were normalized to transcripts per million (tpm). Cell type reference signatures were created from two different published single cell RNA-seq studies: sig1 from GSE229413 and sig2 from GSE205049. Annotations from the original publications were used to define cell types and to subset cells to those from tumor samples (not normal tissue). A final, consolidated reference signature matrix consisting of 23 distinct cell-type profiles was used. We filtered each signature set to remove genes with all zeros in all samples and then filtered each data set to remove any genes with mean expression across all samples that waere less than the overall median gene expression value. This was done to ensure that the most robust genes were kept for the signatures. We kept only the genes that were found in both scRNA-seq signatures to have the same set of genes. We then used granulator (Pfister S et al., 2024) to run deconvolution for each bulk RNA-seq data set with each scRNA-seq reference signature. The granulator R package is designed to run multiple deconvolution algorithms, and we report the non-negative algorithm results for the methods defined by the granulator package as dtangle, nnls, and qprogwc. We used these three different deconvolution methods to compare cell type proportion estimates for both runs of RNA-seq separately.

The use of 23 established cell-type signatures with a limited sample size of 16 tumor samples was a necessary methodological choice. It is important to clarify that cell type deconvolution is not a standard regression model. The 23 cell-type signatures (features) represent fixed, biologically validated reference profiles derived from large external single-cell studies. The deconvolution algorithms (nnls, dtangle, qprogwc) apply constrained linear models (e.g., non-negative and sum-to-one constraints) to estimate the mixing proportions of these known, fixed profiles within the bulk data. These models do not learn or fit the 23 signatures within our samples. We acknowledge that the low sample size limits the power to detect small differences in cell fractions and that the results are susceptible to single-sample variation. Therefore, our findings are presented as preliminary and hypothesis-generating, aiming to establish consistency across three independent deconvolution algorithms and two technical replicates (Mayo and TCGA RNA-seq runs) to demonstrate the robustness of the core immune/stromal shifts.

We evaluated the impact of batch correction and normalization on median expression. We present in the main manuscript the log2 normalized study-specific gene expression values per gene by calculating tpm across all shared genes and then compare median expression values of key pancreatic cancer genes between studies. We use the ComBat-seq function in the R package sva to adjust for batch effect using an empirical Bayes framework. EdgeR to was used to normalize the batch corrected counts within gene and transformed the values by the log2 function. We selected KRAS, TP53, SMAD4, and CDKN2A, as these are the genes which are commonly mutated in PDAC (Hu J. X. et al., 2021). We also selected CTNNB1, JUN, SMAD3, SMAD7, and TCF7, as these transcription factor enrichment identified genes have been observed to have altered gene expression in PDAC compared with adjacent normal tissue (Atay, 2020).

Using a reference scRNA-seq study (GSE205049), we evaluated median differences in cell type proportion between the PDAC sample and the adjacent normal paired tissue samples from 9 patients using a Wilcoxon Rank Sum Test for 23 cell types. After deconvolution, we created scatter plots to show the correlation across the two RNA-seq tumor samples for each patient and evaluated significant differences in median cell type proportion using a paired Wilcoxon rank sum test for 21 identified cell types across the two paired samples. We used the ESTIMATE algorithm implemented in the tidyestimate package in R to estimate tumor purity (stromal and immune cell signatures) using the normalized RNA-seq data from each study separately. Additionally, we evaluated median differences using a paired Wilcoxon rank sum test across the selected nine genes previously determined to be important in PDAC. We visualized differences in log2 normalized gene expression using violin plot and statistically evaluated the significance of the median difference using a Wilcoxon Rank Sum Test with Holm adjusted p-value with values < 0.05 indicating significance. We also explored the effects of batch correction, normalization, and transformation on median differences. A plot was also generated to show correlations between each of the genes based on Spearman correlation coefficients.

We observed that average cell type proportion varied the most between studies (i.e., between samples for each patient) for NK and macrophages (Figures 1A,B). We visually observed across all plots that Macrophages, NK cells, T cells, cycling, and dendritic cells are among the cell types that vary the most. Based on the two single cell reference samples, we observed that NK cells are significantly (adj p-value <0.05) enriched in PDAC tissue compared to adjacent normal tissue while Dendric cells (DC1 and 2) and CD16⁺ monocyte populations are significantly reduced (Figure 2A). Using the Wilcoxon test in our paired samples, we observe the most significant difference in the deconvolution-based cell types are for CD4 (adjusted p-value = 0.00613), endothelial (adjusted p-value = 0.00000172), fibroblasts (adjusted p-value = 0.000161), granulocytes (adjusted p-value = 0.000895), neural (adjusted p-value = 0.000378), NK (adjusted p-value = 0.00224), and monocytes (adjusted p-value = 0.00816; Figure 2B).

We estimated tumor infiltration scores in each study separately using normalized RNA-seq data in the ESTIMATE function in R. In summary, this plot (Figure 3) demonstrates the positive relationship between the estimated stromal and immune cell infiltration and the ESTIMATE score, suggesting that tumors with high stromal and immune cells also likely have high infiltration of the tumor. This is expected because a higher presence of non-tumor cells will naturally lower the proportion of tumor cells in the sample. The blue-green (TCGA) and salmon (Mayo) circles highlight specific tumor samples with distinct characteristics in terms of their microenvironment. It is important to note that ESTIMATE scores can only be interpreted relatively, and it cannot be inferred across studies.

Violin plots of gene expression in the paired Mayo and TCGA sample data for the four PDAC specific genes comparing the effect of batch correction, normalization, and transformations are included (Supplementary Figure S1a,b). The plots highlight the importance of adjusting both within sample (normalization) and across sample (batch correction) to avoid bias. When looking at the first four PDAC specific genes, we did not observe significant median expression differences across paired samples. However given our small sample size and low power, these plots suggest there may be potential differences for KRAS (adj. p-value = 0.18) and SMAD4 (adj. p-value = 0.59, Figure 4A). When evaluating the five transcription factor enriched genes, there were no statistically significant results. Visually, CTNNB1(adj. p-value = 0.30) and TCF7 (adj. p-value = 0.56) varied the most between the paired TCGA and Mayo sample data (Figure 4B).

We used the median as a cut point for each gene and calculated one score for each patient separately to determine concordance (Supplementary Table S1). We observed a substantial concordance (kappa = 0.75) for JUN, moderate concordance for KRAS, TP53, SMAD3 and CDKN2A (Kappa = 0.5) and low concordance value for TCF7, SMAD7, SMAD4, CTNNB1(Kappa = 0.25) across the two studies suggesting that sufficiently representative genomic data cannot be collected with one sample. Likewise, the correlation of the gene expression values for the nine selected genes in general vary significantly across the paired samples (Figure 5). In this correlation matrix, each row and column represent a gene and the number in the cell represents the Spearman correlation coefficient with colors ranging from highly positively correlated values (blue) to highly negatively correlated values (red). The coefficients which reached statistical significance ranged from <-0.24 to >0.24 and are shown in the upper triangle. The hierarchical clustering method was used to order the genes in the matrix. We observed moderate positive correlations between TP53/SMAD4 (0.47) among only TCGA samples; SMAD7/CTNNB1 (0.69), KRAS/CTNNB1 (0.36), and SMAD4/CTNNB1 (0.31) using only Mayo samples. We observed a moderate negative correlation between TCF7/CDKN2A (−0.24) within the TCGA samples only. We observed moderate positive correlation across study genes for CDKN2A/CDKN2A (0.52), SMAD4/TCF7 (0.47), and SMAD3/CDKN2A (0.40). Moderate negative correlations were observed for the following genes across studies: SMAD4/CDKN2A (−0.64), TCF7/SMAD3 (−0.46), SMAD7/TP53 (−0.31), and TCF7/CDKN2A (−0.24). These results could suggest that the TME is likely different across the two samples.