Vascular invasion-associated gene expression is detectable in pre-surgical biopsies of stage I lung adenocarcinoma

Domain: biology
Published: 2026-03-24
Editorial outcome: minor_revision

Abstract
Microscopic vascular invasion (VI) predicts recurrence and benefit from lobectomy in stage I lung adenocarcinoma (LUAD) but cannot be accurately predicted before surgery. Thus, biomarkers that identify this aggressive tumor subset are needed. Here, we show that VI in stage I LUAD is associated with reproducible gene expression programs detectable beyond the invasive focus. Using bulk RNA sequencing of 162 resected tumors and spatial transcriptomics in a subset, we identify and characterize a VI-associated gene signature. A predictor derived from this signature is associated with VI and recurrence in an independent validation cohort and is robust to intra-tumor heterogeneity in multi-regional sampling data. In a cohort of pre-surgical biopsies, predictor scores correlate with matched resections and show promising discrimination of VI. These findings indicate that VI-associated transcriptional changes extend across the tumor and are detectable in limited biopsy material, supporting further validation for preoperative risk stratification in stage I LUAD.

Introduction
Lung adenocarcinoma (LUAD) is the most common lung cancer subtype, and invasive LUAD represents 70–90% of surgically resected lung cancers1. Improved early detection through computed tomography (CT)-based screening programs is expected to increase the proportion of LUAD diagnosed at an early stage. Microscopic vascular invasion (VI), defined as tumor invasion within the lumen of veins or arteries, is a well-described route to metastatic dissemination and is consistently associated with higher rates of tumor recurrence among stage I LUAD2–7. VI is not included in the current World Health Organization (WHO) grading system for lung cancer, but may be a better predictor of recurrence than the most severe WHO-2021 grade, leading to our proposal for VI+ LUAD to be reclassified as distinct angioinvasive LUAD8.

Patients with stage I VI+ LUAD may benefit from adjuvant therapy9–12, but it is difficult to assess in resected tumor specimens. Elastics stains can be used to improve visualization of invaded blood vessels over hematoxylin and eosin (H&E), but comprehensive tumor histopathology review for small (<1 mm) VI foci is difficult and prone to false negatives. A biomarker of VI would improve the existing routine histopathologic review of high-risk cases at resection. Moreover, biopsy specimens do not provide enough tissue material for pathologists to identify VI+ LUAD prior to surgery, preventing the use of VI status to guide surgical or neoadjuvant treatment approaches. Retrospective analysis shows that VI+ patients receiving sublobar resection have increased rates of recurrence13. While recent evidence from large clinical trials, including JCOG0802/WJOG4607L and CALGB140503, suggests that lobectomy provides no clinical benefit over sublobar resection for patients with stage IA disease14,15, post-hoc analyses of pathologic factors specifically associated with recurrence in patients who underwent sublobar resection have yet to be published16,17. Therefore, there is an emerging need to identify if stage I patients will benefit from precision surgery. In parallel, biomarker approaches could also improve patient selection for neoadjuvant or adjuvant therapy in early-stage LUAD18,19.

Molecular profiling technologies are common in clinical pathology, but prior RNA-sequencing (RNA-seq) studies of LUAD have primarily focused on identifying signatures of outcome, which lack reproducibility20,21. Others have identified signatures associated with high-grade growth patterns or subtypes22,23, such as solid LUAD, but they have not been clinically deployed. To our knowledge, VI-associated gene expression has not been studied in lung cancer. Pathologists infrequently document VI as a separate entity from lymphatic invasion (LI), preferring to designate the presence of either type as lymphovascular invasion (LVI). This obfuscates the independent prognostic value of VI24 and complicates downstream molecular approaches seeking to disentangle these modes of tumor spread, the differences of which are still incompletely understood25. Spatial transcriptomics (stRNA-seq) allows for resolving transcriptomic changes associated with tumor invasion within the spatial context of the tumor and tumor microenvironment26. This represents an opportunity to elucidate VI+ LUAD biology with the aim to improve detection of these tumors.

In this study, we profile 163 tumors with detailed histopathologic annotations by RNA-seq, including 15 by stRNA-seq, from a diverse multi-institutional cohort of stage I LUAD patients to comprehensively define the transcriptional landscape of angioinvasive LUAD. We derive a molecular signature associated with VI+ LUAD using bulk RNA-seq, describe its association with histopathology in situ using stRNA-seq, develop and validate a machine-learning predictor of VI, evaluate its stability across intra-tumor heterogeneity (ITH), and demonstrate its proof-of-concept feasibility in a pilot cohort of pre-surgical tumor biopsies.

Results
VI is the invasion type most associated with stage I LUAD recurrence
We assembled a retrospective discovery cohort of stage I LUAD tumors (n=192) from patients across two institutions, Boston Medical Center (BMC), a large urban safety-net hospital, and Lahey Hospital and Medical Center (LHMC), a suburban hospital, a subset of a previously published cohort with detailed clinicopathology annotations (Table S1)8. Tumors were graded using our novel grading system, whose reproducibility has been validated by us and others8,13,27,28. This novel grading system classifies tumors into low malignant potential (LMP), no special type (NST), or VI, and is more associated with outcome than the WHO grading system (Supplementary Fig. 1a, b). VI was the strongest predictor (hazard ratio (HR)=7.62, 95% CI 3.70–15.8, Cox regression p=5.17×10−8) of 7-year recurrence-free survival (RFS) among the four observed modes of tumor spread (VI, LI, visceral pleural invasion (VPI), and spread through air spaces (STAS)), in agreement with previous reports of its prognostic value in stage I LUAD8,29–31 (Supplementary Fig. 1c–g). VI was present in 30% of tumors, and 84% of VI+ tumors contained multiple types of invasion (Supplementary Fig. 1h). Regardless, VI remained a significant predictor of recurrence when including clinical covariates and other invasion types (HR=11.04, 95% CI 3.80–32.1, multivariate Cox regression p=1.01×10−5; Supplementary Fig. 1i). Of the different invasion types, only LI was also an independent predictor of recurrence (HR=2.73, 95% CI 1.15–6.44, multivariate Cox regression p=0.02). When scoring tumors based on the proportion of LUAD growth patterns, VI was associated with solid (Bonferroni adjusted P value (p.adj)=0.004, Wilcoxon test), but not micropapillary (p.adj=0.08) or cribriform proportion (p.adj=0.35) (Supplementary Fig. 1j). In contrast, STAS was most strongly positively associated with increased micropapillary proportion (p.adj=5.3×10−7), as previously reported32. VI+ tumors were more strongly associated with distant recurrence (sub-distribution HR=19.1, 95% CI 4.27–85.0, multivariable Fine-Gray regression p=1.1×10−4) than loco-regional recurrence only (sub-distribution HR=3.86, 95% CI 1.22–12.2, multivariable Fine-Gray regression p=0.02), consistent with tumor spread through the vascular system (Supplementary Fig. 1k–m).

Gene expression changes in stage I LUAD with VI
To identify gene expression changes in tumors with VI, we profiled a subset of surgically resected stage I LUAD tumors from the discovery cohort (n=108, n=103, post-QC) by bulk RNA-sequencing (RNA-seq) (Fig. 1a, Table S2). Among tumors that passed quality control, 78 were VI- and 25 were VI+. We identified 474 genes differentially expressed between VI+ tumors and tumors of LMP using the three-level factor of tumor grade (LMP, NST, VI) as the independent variable at a significance threshold of FDR<0.01, which we clustered into four gene expression clusters (Fig. 1b). We selected the optimal number of clusters as four via consensus clustering (Supplementary Fig. 2a). Genes in three of the clusters had increased expression in LUAD with VI (clusters 1, 2, and 3) and genes in the other cluster had decreased expression (cluster 4). Pathway enrichment analysis using EnrichR (FDR<0.05) revealed that genes involved in different biological processes were enriched in each cluster (Fig. 1c). Cluster 1 (n=115 genes) was enriched for genes involved in the cell cycle. In contrast, cluster 2 (n=37 genes) was enriched for genes involved in tissue remodeling and vasculogenesis, including the epithelial to mesenchymal transition (EMT), extracellular matrix organization, and angiogenesis. The genes in cluster 3 (n=182 genes) were enriched for some pathways that overlapped with cluster 1, like mTORC1 signaling and E2F targets, but was distinguished by enrichment for pathways related to response to reactive oxygen species and hypoxia. Finally, the genes with decreased expression in VI tumors (cluster 4; n=140 genes) were enriched for pathways indicating reduced malignant progression (p53 pathway, regulation of cell growth), increased cell-cell adhesion via TGF-β signaling (regulation of pathway-restricted SMAD protein phosphorylation), and increased immune surveillance (IL-2/STAT5 signaling) in VI− tumors. In the discovery cohort, the mean z-score of the genes in each of the four clusters was a significant predictor of VI (Supplementary Fig. 2b) and was also able to distinguish between NST and VI+ tumors (Supplementary Fig. 2c), supporting that the differentially expressed genes reflect VI-specific biology. We observed that two VI+ tumors clustered with the LMP classified tumors. These both represented rare edge cases, with one being the only lepidic predominant VI+ sample and the other being the only VI+ sample with >50% micropapillary growth pattern. Larger cohorts containing more of these types of samples will be required to determine if the current signature might be less accurate in detecting VI in these conditions.

Given that VI and LI are routinely reported together as LVI8,24, we compared the VI signature with gene expression differences associated with LI in stage I LUAD. We added LI as a covariate to our model to determine whether the gene expression differences seen with VI are due to the covariance between LI and VI. We identified 133 genes associated with VI (VI+ vs. LMP, FDR<0.01), which included 130 genes from the original 474-gene VI signature (Supplementary Fig. 3a). At the same significance threshold (FDR<0.01), we recovered only 6 genes associated with LI (LI+ vs. LI-) using the same model, with 2/6 belonging to the VI signature (Supplementary Fig. 3b), suggesting the biological signal for LI is much weaker. To further interrogate the relationship between LI and VI-associated gene expression, we explored whether the various VI gene clusters from the 474-gene signature are enriched among the genes most differentially expressed in LI+ tumors using gene-set enrichment analysis (GSEA). Interestingly, we found that the genes in VI cluster 4 are significantly enriched among the genes decreased in LI+ tumors (Supplementary Fig. 3f). But unlike the genes in VI cluster 1 (cell-cycle) and VI cluster 3 (hypoxia), the genes in VI cluster 2 (EMT/angiogenesis), were not significantly enriched among the genes most increased in LI+ tumors (Supplementary Fig. 3c–e). Finally, genes previously found to be increased in VI+ endometrial cancer33 or LVI+ breast cancer34 were each significantly enriched among genes most increased in VI+ stage I LUAD tumors using GSEA (Supplementary Fig. 3g, h). This hints at common molecular changes associated with VI across epithelial tumors that should be explored in future work.

VI gene clusters are associated with specific LUAD histopathology features in stRNA-seq data
Although the genes differentially expressed in tumors with VI were enriched for biological processes implicated in tumor intravasation, the RNA for sequencing was isolated from histologic sections selected to be representative of the predominant histologic pattern and were not selected to include invaded blood vessel(s). Given that VI foci represent such a small percentage of tumor volume (diameter ≤1 mm), we hypothesized that the VI signature we identified reflected molecular changes extending beyond the site of intravasation. To assess the spatial architecture of the expression of the VI signature and its association with LUAD histopathology features, we profiled 15 (post-QC) resected stage I and stage II LUAD samples (n=8 VI−, n=7 VI+) from 13 patient tumors by high-resolution stRNA-seq using the 10× Genomics Visium platform, including 3 sections containing invaded vessels directly in the 6.5mm2 Visium capture area (Table S3). Of the VI− tumors, 7/8 were NST. All six non-mucinous LUAD growth patterns (solid, cribriform, micropapillary, acinar, papillary, lepidic) were represented across the stRNA-seq capture areas, along with three invasion types (VI, VPI, and STAS) (Fig. 2a). Adjacent normal appearing lung was present in 12/16 capture areas. The histopathology associated with VI+ tumors in the stRNA-seq data was generally representative of trends observed in our larger clinical cohort, although the micropapillary pattern tended to be underrepresented in VI+ tumors (p.adj=8.62×10−14) and the papillary pattern was overrepresented (p.adj=1.01×10−11) (Supplementary Fig. 4a). A solid growth pattern and the presence of desmoplastic stroma were both associated with angioinvasion (p.adj=6.38×10−31; p.adj=2.26×10−18). While VI co-occurred with other invasion types frequently in our clinical cohorts, we did not find more than one invasion type per capture area in the stRNA-seq dataset. In summary, our stRNA-seq dataset captures many of the pathology features that are know to occur in stage I LUAD and co-occur with VI.

After stRNA-seq quality control and filtering, we obtained expression estimates from 43,421 spots (50μm diameter) with a median depth of 8548 unique molecular identifiers (UMIs)/spot (Supplementary Fig. 4b) that detected a median 3996 genes/spot (Supplementary Fig. 4c). The low depth of sample 6 was likely due to tissue detachment that occurred during the stRNA-seq workflow, and this sample was removed prior to downstream analysis. Strong concordance was observed between pseudo-bulked stRNA-seq gene expression counts and bulk RNA-seq counts of matched tumor tissues across all genes (spearman R=0.85) (Supplementary Fig. 4d) and across VI signature genes (spearman R=0.78) (Supplementary Fig. 4e).

We next scored individual spots to create an enrichment score representing the activity of the genes from each of the four clusters in our bulk VI signature. Across all capture areas, spatially weighted correlation analysis of enrichment scores revealed that individual VI clusters were spatially distinct from one another (Supplementary Fig. 4f). These data suggest the distinct patterns of gene co-expression represented by the four clusters are also spatially distinct, indicating that the biological processes that these gene expression patterns represent are spatially variable and are unlikely to represent redundant biology.

To investigate if the VI gene clusters are more strongly expressed in regions with different histopathologic feature annotations in the stRNA-seq data, we modeled the enrichment score of each cluster as a function of histopathologic annotation using a linear mixed model with the section as a random effect. Unexpectedly, despite the VI gene clusters being derived in bulk RNA-seq, we observed a strong increase in expression directly in VI foci of both cluster 1 (p.adj=9.02×10−28) and cluster 3 (p.adj=8.92×10−73). Cluster 1 and cluster 3 were also increased in regions annotated as high-grade patterns (solid (p.adj=8.87×10−15; 7.90×10−75) micropapillary (p.adj=3.94×10−07; 1.40×10−23) and cribriform (p.adj=4.23×10−08; 2.39×10−58)) (Fig. 2b). Cluster 2 was most noticeably expressed in desmoplastic stroma (p.adj=1.22−122) with a greater than 2.95-fold higher mean enrichment score than in stroma (p.adj=1.87×10−25). Conversely, expression of cluster 4 was strongly decreased in VI foci (p.adj=5.79×10−79), followed by desmoplastic stroma (p.adj=6.53×10−50) and showed strong positive enrichment in adjacent normal-appearing lung (p.adj=2.57×10−04), stroma (p.adj=8.90×10−12), and papillary (1.68×10−13). When enrichment scores were visualized directly on tissue sections, the independent spatial patterning of each cluster was apparent (Fig. 2c). Given that clusters 1 and 3 were significantly enriched directly in invasive foci in addition to high-grade patterns, we were interested to see whether the expression of these genes was also higher in VI+ tumors independently of aggressive pattern. When we analyzed only spots annotated with high-grade LUAD patterns, we found that cluster 1 (p=0.047), 2 (p=0.004) and 3 (p=0.003) were all higher and cluster 4 (p=6.9×10−5) was lower in spots annotated with high-grade patterns belonging to VI+ tumors, suggesting that our VI signature is capturing an angioinvasive phenotype that is independent of aggressive pattern and the invasive focus (Fig. 2d). In total, these findings suggest that gene expression clusters of angioinvasive LUAD identified by bulk RNA-seq are differentially expressed between distinct histopathologic features in situ, show increased expression in the same tumor patterns if it is a VI+ tumor, and are not solely expressed in invaded blood vessels.

The VI signature is composed of both tumor-specific and tumor-microenvironment changes reflective of angioinvasion
Our supervised analysis above relied upon detailed pathologic annotation of the tissue containing Visium spots, but spots were only labeled if there was a clear consensus on morphology, especially for the histologic pattern. In some cases, tumor areas may not have been annotated if they did not clearly contain canonical features of one of the six histologic patterns35. Additionally, for poorly differentiated patterns such as solid, the annotation may include admixed stromal or immune cells. With the aim to further delineate the spatial biology of the VI clusters, we estimated the cell type composition of each stRNA-seq spot with CytoSPACE by using a published single-cell RNA seq atlas of 295,813 annotated cells from 124 stage I LUAD samples (primary tumor and normal adjacent lung) and segmenting nuclei per spot with StarDist36–39. LUAD pathology annotations showed distinct cell type composition (e.g., B cells with lymphoid aggregates) (Supplementary Fig. 5a). When we modeled predicted cell type proportion per spot as a function of histopathologic annotation using a linear mixed model with the section as a random effect, we observed a significant enrichment of predicted B cells, Plasma cells, and other immune cell types with lymphoid aggregates (Fig. 3a). Normal bronchus had a significantly higher proportion of predicted ciliated cells, as expected. Spots annotated as VI foci had significantly higher proportions of predicted tumor cells, while desmoplastic stromal regions had significantly higher proportions of peribronchial fibroblasts. We found a trend toward a higher proportion of peribronchial fibroblasts belonging to VI+ stRNA-seq samples (Fig. 3b), consistent with the tissue remodeling biology of VI cluster 2 described above. Peribronchial, but not adventitial or alveolar fibroblasts, showed a significantly higher enrichment for a published NSCLC myofibroblast-cancer-associated fibroblast (CAF) signature from Hanley et al.40 (Fig. 3c). Hanley et al. reported key markers of LUAD myofibroblasts as MMP11, POSTN, CTHRC1, COL1A1, and COL3A1, which are all VI cluster 2 genes. Consistently, VI cluster 2 was significantly more highly enriched in this myofibroblast cluster than with alveolar and adventitial clusters in LUAD scRNAseq samples from Hanley et al (n=39) (Fig. 3d). Deconvoluted VI+ stRNA-seq samples also had a trend toward a higher proportion of plasma cells (Fig. 3b). To validate these findings, we deconvoluted the bulk RNA-seq discovery cohort data using the same atlas from Salcher et al. Although the predicted plasma cell proportions were associated with pathologist-annotated plasma cell grade on case-matched H&E images (p=1.37×10−8) (Supplementary Fig. 5b), we did not observe an association with VI+ tumors (p.adj=1), suggesting that plasma cell infiltration may be overrepresented in the stRNA-seq data. In contrast, in the same bulk RNA-seq data we found a significantly higher predicted proportion of peribronchial fibroblasts (myofibroblasts) (p.adj=0.005) in VI+ tumors (Fig. 3e). To examine the associations between the expression of the VI gene clusters and predicted cell types in the stRNA-seq data, we first generated cell type specific gene signatures from the Salcher atlas (Supplementary Fig. 5c). We then used spatially weighted regression to identify correlation between cell type signatures and our VI gene clusters across the stRNA-seq data after removing low abundance cell types. Unsupervised clustering of the average spatially weighted correlation across all tumors revealed co-localization of cluster expression with distinct cell types (Fig. 3f). As expected from our previous observation of enrichment in desmoplastic stroma, VI cluster 2 was most spatially correlated with peribronchial fibroblasts, as well as other stromal cell types, including endothelial and smooth muscle cells. In contrast, VI clusters 1, 3, and 4 were all tightly spatially correlated with epithelial cells, suggesting the contribution of these gene clusters to the VI signature is dependent more on cell state than cell type. These associations were not as apparent when examining enrichment scores in the scRNA-seq atlas (Supplementary Fig. 5d), highlighting the added value of spatial deconvolution. Collectively, these results show that the VI gene clusters capture biological changes associated with angioinvasion from both tumor cells and the tumor microenvironment.

To further validate our stRNA-seq deconvolution findings about the relationship between VI cluster 2, desmoplastic stroma, and peribronchial fibroblasts, we performed COL1A1 RNAscope ISH and COL1A1 IHC labeling on serial sections close to stRNA-seq sample 4, a VI+ tumor with prominent desmoplastic stroma. We selected COL1A1 as it is a representative gene from VI cluster 2 that also showed significantly higher expression in myofibroblasts compared with adventitial or alveolar fibroblasts in Hanley et al (Supplementary Fig. 6a). After serial image registration to the H&E from the section used for stRNA-seq, we observe a significant correlation across spots between the fraction of ISH RNAscope COL1A1+ expressing cells per spot and the average VI cluster 2 gene expression per spot (r=0.53, p<2.2×10−16; Supplementary Fig. 6b). We also observe a similar significant correlation between the average VI cluster 2 gene expression per spot and the fraction of COL1A1+ cells per spot labeled by IHC (r=0.47, p<2.2×10−16; Supplementary Fig. 6c). We also performed IHC for THY-1, a VI cluster 2 gene that was not reported by Hanley et al as a key myofibroblast marker. The fraction of THY-1+ expressing cells per spot by IHC was significantly correlated with the fraction of COL1A1+ expressing cells measured by IHC, which further supports the co-expression of VI gene cluster 2 at the protein level (R=0.52, p<2.2×10−16; Supplementary Fig. 6d). When we examined the correlation between COL1A1+ cell fractions measured by either IHC or ISH with the fraction of myofibroblasts predicted by the CytoSPACE stRNA-seq deconvolution, we found that both were significantly more correlated with the fraction of peribronchial/myofibroblast CAFs than either adventitial or alveolar fibroblasts (Supplementary Fig. 6e). This further demonstrates that cells expressing COL1A1, a key marker of VI cluster 2, show a myofibroblast-like phenotype. Finally, we also observed a significant association between the fraction of COL1A1+ IHC fraction and spots annotated as desmoplastic stroma (Supplementary Fig. 6f), further suggesting that VI cluster 2 is associated with tissue remodeling.

To further generalize these findings beyond the stRNA-seq sample, and to directly demonstrate the expression of COL1A1 in myofibroblasts, we analyzed an independent cohort (Table S4) of 20 resected stage I LUAD samples (11 VI−, 9 VI+) using duplex IHC for COL1A1 and POSTN. POSTN is a canonical myofibroblast marker40,41 and a VI cluster 2 gene. In pathologist-annotated stroma regions (i.e., non-desmoplastic stroma), nearly all POSTN+ cells were also COL1A1+ (Supplementary Fig. 6g, h). The fraction of COL1A1+POSTN+ cells in these regions was significantly higher in VI+ than VI- tumors (median 68% vs 32%, p=0.006), consistent with VI-associated myofibroblast expansion (Fig. 3g). In analyzing desmoplastic stroma, we restricted our analysis to desmoplastic stroma from VI+ tumors, given the paucity of desmoplastic stroma in VI− tumors. Within desmoplastic stroma, POSTN+ cells are also nearly universally COL1A1+ (~100%, Supplementary Fig. 6i), confirming persistent co-expression within myofibroblasts. The fraction of all cells in desmoplastic stroma regions that are COL1A1+ was significantly higher than outside the desmoplastic stroma regions (median 87% vs 57%, p=7.09×10−7) (Fig. 3h), without a corresponding rise in POSTN+COL1A1+ cells (median 26.3% vs. 21.0%, p=0.263) (Supplementary Fig. 6j). This is a result of an increase in COL1A1+POSTN- stromal cells (median 61% vs. 39%, p=0.02) (Supplementary Fig. 6k). These findings support and extend the stRNA-seq results, demonstrating that the elevated expression of VI cluster 2 in VI+ tumors reflects both the VI-specific expansion of POSTN+COL1A1+ myofibroblasts in non-desmoplastic stroma, and the VI-associated expansion of desmoplastic stroma in which the majority of COL1A1+ cells are POSTN-.

A VI predictor derived from the signature validates in an independent stage I LUAD cohort
Our finding of biologically and spatially distinct gene expression clusters associated with VI+ LUAD suggests that all four clusters provide orthogonal information and might be combined to form a predictor of angioinvasive stage I LUAD. To test this hypothesis, we proceeded with a standard machine learning cross validation pipeline using only the samples from the discovery cohort (n=103; Fig. 4a). Training was performed in the discovery cohort (n=103) where the differentially expressed genes were first identified, using a cross-validation approach (Fig. 4a). Predictors utilizing a binomial logit generalized linear model (GLM) with ridge regression performed the best on the internal cross-validation test sets within our discovery cohort and we observed no dependence of gene set size on model performance (Supplementary Fig. 7a, b). When we trained a GLM with ridge regression on the full discovery cohort, a 48-gene predictor achieved an area under the receiver operating characteristic curve (AUROC) of 0.85 (95% CI 0.78–0.93, p=1.1×10−7) to separate VI+ and VI- LUAD (Fig. 4b). The predictor was independently associated with 7-year recurrence free survival (HR=3.99, 95% CI 1.77–8.99, Cox regression p=8.29×10−4) when including common clinicopathology variables, including TNM stage (Supplementary Fig. 7c). When ranked by feature importance in the model, the top 3 genes in the predictor were from cluster 3 and included MUC16, NOTCH3, and H19, but only cluster 2 genes were significantly enriched among the top ranked genes overall (two-sided Kolmogorov-Smirnov test p=0.02; Fig. 4C). The top features for clusters 1, 2, and 4 were MKI67, COL1A1 and CHD2, respectively.

The predictor retained a similar performance in predicting VI+ LUAD when applied to our independent validation cohort (Table S5; n=42 VI−, n=17 VI+) with an AUROC of 0.86 (95% CI 0.76–0.96, p=6.16×10−6) (Fig. 4d). It retained equivalent performance even after excluding LMP tumors (NST vs. VI only, AUROC 0.86, 95% CI 0.75–0.96, p=1.22×10−5) (Fig. 4e). It was associated with histologic grade (Supplementary Fig. 7d) and 7-year recurrence free survival (HR=1.98, 95% CI 1.09–3.60, Cox regression p=0.025) (Fig. 4f). Additionally, it was also a significant predictor of RFS among VI- LUAD tumors (HR=2.76, 95% CI 1.41–5.41, Cox regression p=0.003), and in tumors without VI, necrosis was associated with higher VI predictor scores in a univariate linear model (p.adj=0.003) (Fig. 4g). In our validation set, the VI predictor retained good performance across tumor subsets that varied by histological growth patterns, with the best performance in samples that contained no lepidic pattern (0% lepidic, n=22, AUROC 0.96, 95% CI 0.89–1.0, p=1.18×10−4) and the worst performance in samples that contained no solid pattern (0% solid, n=34, AUROC 0.70, 95% CI 0.48-0.93, p=0.198) (Fig. 4h). It also showed variable performance in predicting other invasion types, and notably was not significantly better than a random classifier at detecting LI+ tumors (AUROC 0.53, 95% CI 0.38–0.69, p=0.677) or LI+ VI- tumors (AUROC 0.50, 95% CI 0.31–0.69, p=1.0) (Fig. 4i), supporting our hypothesis that vascular invasion is a different molecular process than lymphatic invasion. Although we did not have separate VI annotations for published LUAD RNA-seq datasets, in TRACERx, a longitudinal study of NSCLC, we observed a difference in predictor scores between LVI+ and LVI− LUAD (p=0.0014) (Supplementary Fig. 7e). Interestingly, VI predictor scores were higher in tumors from patients with detected preoperative ctDNA (p=5.9×10−4; Supplementary Fig. 7f). Higher VI predictor scores were significantly associated with decreased 5-year RFS in TRACERx stage I LUAD (Fig. 4f; HR=1.69, 95% CI 1.04–2.75, Cox regression p=0.036) but were not significantly independent of TNM stage IA and IB (HR=1.60, 95% CI 0.97–2.65, Cox regression p=0.066). However, the predictor was significantly associated with 5-year RFS independently of TNM stage I and II (Supplementary Fig. 7g; HR=1.46, 95% CI 1.07–2.0, Cox regression p=0.018). Additionally, elevated VI predictor scores remained significantly associated with decreased 5-year RFS when including surgery type (lobectomy or sublobar resection) and adjuvant therapy as covariates in stage I (HR=1.66, 95% CI 1.02–2.71, multivariate Cox regression p=0.041; Supplementary Fig. 7h). Finally, higher VI predictor scores associated with decreased overall survival in stage IA LUAD tumors from both The Cancer Genome Atlas (TCGA) (n=127) (HR=1.72, 95% CI 1.05–2.81, Cox regression p=0.032) and Uppsala (n=44) (HR=2.11, 95% CI 1.05–4.27, Cox regression p=0.039) cohorts42,43 (Fig. 4f).

The VI predictor is robust to intra-tumor heterogeneity and predicts VI when measured in pre-surgical biopsies from stage I LUAD patients
Our stRNA-seq analysis of the expression of VI-associated genes implies that each VI-associated cluster may be spatially distinct. This suggests that a VI predictor combining genes from each cluster might overcome ITH. Overcoming ITH is crucial for molecular biomarkers that sample only a portion of the tumor volume, such as those measured on tissue available from biopsies44. This is particularly important for detecting VI due to the small size of invaded vessels. When we divided stRNA-seq spots from VI+ samples into distal VI+ (spots >1 mm beyond the edge of invaded vessels) and proximal VI+ (spots <1 mm from and including invaded vessels) (Fig. 5a), we found a significant difference in the enrichment score between VI− and distal VI+ tumors for the predictor genes (p=0.028) (Fig. 5b). This supports the hypothesis that the predictor may be able to detect VI away from the invasive focus. Next, we used RNA-seq from multi-region sampling of stage I LUAD (n=63) in TRACERx to evaluate the variability of the VI predictor to random tissue sampling in a larger cohort32. Although a label for VI was not available for the TRACERx tumors, we observed strong correlation between predictor scores from two randomly selected paired regions (n=126) across all tumors (R=0.87, p<2.2×10−16) (Fig. 5c). When we filtered the dataset to remove lowly expressed genes and ranked the remaining 15,800 genes in the dataset by correlation between paired regions, the 48 genes present in our predictor were significantly enriched among those with the least ITH (p=2.89×10−6) (Fig. 5d). Furthermore, inter-tumor heterogeneity of the VI predictor score (as measured by the absolute difference of scores from regions of two randomly selected tumors) was significantly higher than ITH (as measured by the absolute difference of scores from two randomly selected regions within the same tumor) (p=1.12×10−6) (Fig. 5e). Taken together, these results support our observations made from the stRNA-seq data that our VI signature represents a global shift in molecular aggressiveness toward an angioinvasive phenotype and that these changes can be leveraged to predict VI+ LUAD from any region of the tumor.

Finally, to demonstrate a proof of concept for how our VI predictor may be applied in routine clinical workflows, we performed RNA-seq on an independent cohort of formalin-fixed paraffin-embedded (FFPE) pre-surgical biopsies (Table S6; n=12 VI−, n=12 VI+), collected from stage I LUAD patients treated at a separate institution from our discovery or validation cohorts. The VI predictor measured in these pre-surgical biopsies predicted the detection of VI in the resected tumor with an AUROC of 0.90 (95% CI 0.75–1.0, p=4.96×10−4) (Fig. 5f). The sensitivity and specificity when using the Youden point from the discovery cohort as a cutoff in the biopsies was 0.83 (95% CI 0.52–0.98) and 0.92 (95% CI 0.62–1.0), respectively. Following standard methods for estimating diagnostic accuracy45, we derived posterior distributions for sensitivity and specificity from the biopsy cohort’s classification counts. We then used Monte Carlo simulations and applied Bayes’ theorem to project positive predictive value (PPV) and negative predictive value (NPV) across a realistic prevalence range (15–35%). The VI predictor’s prevalence-adjusted NPV estimate from the biopsies was high and precise (median 91%, CI 81–98% at 30% prevalence) (Fig. 5g). Because our balanced biopsy cohort design under-samples true negatives relative to general clinical prevalence, these CIs are conservative and would gain further precision in an unselected cohort. Conversely, the prevalence-adjusted PPV estimates from this biopsy cohort are less precise (median 73%, CI 47–95% at 30% prevalence), reflecting both the small number of VI+ biopsy cases and the moderate prevalence of VI+ stage I LUAD. In addition, we observed a strong correlation between VI predictor scores determined from pre-surgical biopsies (n=12) and matched resected tumors (r=0.77, p=0.0053) (Fig. 5h).

Methods
Clinical cohorts
Discovery cohort. A discovery cohort consisting of 192 resected tumors from 192 patients with 8th edition TNM stage 0/1 LUAD and not treated with neoadjuvant or adjuvant therapy was included in this study. Cases were from Boston Medical Center (BMC) and Lahey Hospital & Medical Center (LHMC). The BMC and LHMC Institutional Review Boards approved this study (BU/BMC IRB H-37859; Lahey Clinic IRB-518308). Relevant clinical information on all patients included can be found in Table S1. For all clinical cohorts, race was obtained from the electronic medical record and reflects patient self-reported categories. The 192 tumors represented a wide range of sample ages, so a small pilot batch of 12 samples was selected to determine which tumors might yield useful data (4/5 of the samples removed post-QC belonged to this pilot batch). We then selected an additional 96 samples that yielded at least 83 ng of library. RNA-seq was performed on these 108 tumors, with 103 passing quality control. Relevant clinical information on all patients included can be found in Table S2.

Validation cohort. An independent validation cohort consisting of 61 resected tumors from 60 patients with 8th edition TNM stage I/II LUAD and not treated with neoadjuvant or adjuvant therapy was included in this study. All samples were selected to be stage IA/IB at the time of collection, except for one tumor that was upstaged to stage IIA under the 8th TNM edition. Cases were from LHMC. RNA-seq was performed on 61 tumors, with all passing quality control. StRNA-seq was performed on a subset of these 61 tumors and included 16 samples taken from 14 tumors (13 stage I tumors and 1 stage II tumor). Relevant clinical information on all stRNA-seq samples included can be found in Table S3. Only the 8th edition TNM stage I LUAD from the validation cohort was included in the bulk RNA-seq (n=60) analysis. Relevant clinical information on all patients included can be found in Table S5.

IHC cohort. An independent cohort consisting of 20 resected tumors from 20 patients with 8th edition TNM stage I LUAD was included in this study. All samples were from patients not treated with neoadjuvant or adjuvant therapy and were receiving care at Inova Schar Cancer Institute. The Inova Institutional Review Board approved this study (IRB U23-06-5093). Relevant clinical information on all patients included in this cohort can be found in Table S4.

Pre-surgical biopsy cohort. RNA-seq was performed on an independent validation cohort consisting of 24 pre-surgical biopsies and 12 matched resected tumors from 24 patients with 8th edition TNM stage I LUAD. All samples were from patients not treated with neoadjuvant or adjuvant therapy and were receiving care at Inova Schar Cancer Institute. The Inova Institutional Review Board approved this study (IRB U23-06-5093). FFPE pre-surgical biopsy samples were derived from needle biopsy, forceps biopsy, or fine needle aspiration. The determination of novel grade was based upon a representative FFPE block of the resected tumor. Relevant clinical information on all patients included in this cohort can be found in Table S6.

Pathology review
An experienced thoracic pathologist (E.B.) reviewed all pathology cases. An inter-reader reproducibility assessment was previously performed on these pathology cases and found substantial agreement (κ=0.71) among 5 pathologists for the diagnosis of vascular invasion (VI)13. VI was defined as luminal invasion of a vein or muscular artery either within or adjacent to the tumor. Tumor proportions of lepidic, acinar, papillary, micropapillary, and solid patterns were assessed in 5% increments with distinction of simple tubular acinar from complex and cribriform acinar patterns. Adenocarcinoma in situ (AIS) was assigned to purely lepidic tumors ≤3cm, whereas minimally invasive adenocarcinoma (MIA) was diagnosed when non-lepidic foci measured ≤0.5cm as per WHO criteria. WHO-2021 grade was defined as G1, lepidic predominant with <20% high-grade patterns; G2, acinar or papillary predominant with <20% high-grade patterns; and G3, ≥20% high-grade patterns (solid, micropapillary, and/or complex glands)73,74. LMP adenocarcinoma was assigned as previously described27. LMP tumors were non-mucinous adenocarcinoma measuring ≤3cm in total size, with ≥15% lepidic growth, and without nonpredominant high-grade patterns (≥10% cribriform, ≥5% micropapillary, ≥5% solid), >1 mitosis per 2mm2, vascular, lymphatic, or visceral pleural invasion, STAS, or necrosis. AIS/MIA and LMP were analyzed together due to their identical outcome (100% 10-year disease-specific survival). One LMP recurred after wedge resection with a positive surgical margin. The tumor recurred at the staple line and was treated with SBRT, resulting in prolonged survival (over 10 years) without further recurrence or metastasis. NST designation was given for all other tumors not classified as VI or LMP. Stage assignments were retrospectively made using the 8th edition of the American Joint Committee on Cancer (AJCC).

Statistical analysis
All statistical analysis was performed with R version 4.2.1. Tables were created with the tableone package. P values were corrected for multiple hypothesis testing where appropriate using either Holm-Bonferroni adjustment (abbreviated as p.adj) or false-discovery rate correction (abbreviated as FDR), depending on the default method of the R package used, unless otherwise specified.

RNA-Seq library preparation, sequencing, and data processing
Total RNA was extracted from FFPE tissue with the Qiagen AllPrep DNA/RNA Universal Kit (discovery, validation cohorts) or Qiagen miRNeasy FFPE kit (biopsy cohort), and exome-targeted sequencing libraries were prepared with the Illumina TruSeq RNA Exome Library Prep Kit (discovery, validation cohorts) or Illumina RNA Prep with Enrichment Kit (biopsy cohort). Sample sequencing was performed on the Illumina NextSeq 500 to generate paired-end 50-nucleotide reads (discovery cohort), Illumina HiSeq 2500 to generate paired-end 75-nucleotide reads (validation cohort), or Illumina NextSeq 2000 to generate paired-end 100-nucleotide reads (biopsy cohort). Basespace (Illumina) was used for demultiplexing and FASTQ file generation. A Nextflow v21.10.6 RNA-Seq pipeline aligned reads using hg38 and STAR75 v2.6.0c. Counts were calculated using RSEM76 v1.3.1 using Ensembl v108 annotation. Quality metrics were calculated with STAR and RSeQC77. EdgeR was used to compute normalized data (library sizes normalized using the trimmed mean of M-values). Genes with count per million (CPM) >1 in at least 10% of samples were retained for further analysis. Samples in the discovery and validation cohorts were excluded if the transcript integrity number (TIN)78 calculated by RSeQC was >2 standard deviations below the mean TIN. Combat-seq79 was used to correct a batch effect between collection sites (BMC and LHMC) in the discovery cohort. It was also used to correct a batch effect between sequencing batches (batch 1 and 2) within the validation cohort.

Derivation of VI signature and biological annotation
Genes whose expression is associated with the presence of VI were derived from 103 tumors in the discovery cohort using a negative binomial generalized linear model via edgeR80. Our previously published novel grading system (LMP, NST, and VI)8,13 was the independent variable, and we compared VI and LMP groups using glmLRT(). Differentially expressed genes were identified using a false-discovery rate (FDR) threshold of 0.01. Four gene expression clusters were identified using the Ward 2 method of hierarchical clustering of genes and samples based on Euclidean distance81. When we explored other values for k (k=2–10 clusters) using consensus clustering, we found that clustering stability showed substantial improvement (large change in delta area under the CDF) up to but not beyond k=4. Additionally, k>4 did not return distinct significant biological enrichment terms for all clusters, suggesting a lack of biologically distinct gene co-expression with k>4. The top ten biological processes and pathways enriched in each of the four VI clusters were identified using Enrichr82 with queries to the following databases: MSigDB Hallmark 2020 and GO Biological Process 2021. To compare VI and LI-associated gene expression, a VI signature was first rederived after adding LI as a covariate to the edgeR model described above, and differentially expressed genes were again identified using an FDR threshold of 0.01 when comparing VI and LMP groups. To identify genes associated with LI, contrasts were set to LI using makeContrasts() and differentially expressed genes identifying using an FDR threshold of 0.01. Genes associated with LI were ranked by -log10(pvalue)*the sign of the logFC. The VI gene cluster enrichment against this ranked list was calculated using the GSEA() function from the clusterProfiler83 R package, and p-values were adjusted with p.adjust() for Bonferroni correction.

Spatial transcriptomics library preparation, sequencing, and data processing
Sixteen samples from 14 patient tumors within the validation cohort were selected for 10× Genomics Visium for FFPE spatial whole-transcriptome profiling. All samples were selected to be stage IA/IB at the time of collection, except for one that was upstaged to stage IIA under the 8th TNM edition (this sample was excluded from the bulk RNA-seq validation cohort when evaluating the performance of VI predictor scores). Samples were chosen based on (1) the presence of pathological features of interest (e.g., VI foci, representative LUAD histologic patterns) and (2) more than 50% of RNA fragments being greater than 200 nucleotides (DV200) after extraction with the AllPrep DNA/RNA Universal Kit (Qiagen) during bulk RNA-seq library preparation. New 5μm tissue sections were cut by a microtome, placed in a 42°C RNase-free water bath, and manually transferred onto 6.5×6.5mm tissue capture areas on Visium Spatial Gene Expression Slides (PN-1000185, 10× Genomics). Tissues were deparaffinized, H&E stained, imaged with a Leica Aperio AT2, and de-crosslinked according to the manufacturer’s recommended protocol, with the following user modifications made to deparaffinization to prevent tissue detachment: (1) the 15 min incubation step during deparaffinization was removed, (2) after immersion in 96% ethanol, one 85%, one 70% and one 50% ethanol immersion for 3 min each was added. Probe hybridization was performed with the Visium Human Transcriptome Probe Kit (PN-1000364), followed by probe ligation, extension, and elution for downstream qPCR cycle determination and 10× library construction with the Visium FFPE Reagent Kit (PN-1000362). Sequencing was done on an Illumina NextSeq 2000. The 10× Space Ranger pipeline was used to demultiplex and generate FASTQ files. FASTQ files were input into the Space Ranger FFPE count algorithm along with matching H&E.tiff images for alignment to the reference probe set and generation of count matrices for each Visium capture area.

Spatial transcriptomics pathology annotations
H&E images of sections captured for stRNA-seq analysis were annotated for pathological features of interest by an experienced thoracic pathologist (EB) using Loupe Browser v6 (10x Genomics).

Spatial transcriptomics data analysis
stRNA-seq data were processed using the Seurat R package, unless otherwise specified84. Sample 6 was removed due to tissue loss that occurred during wash steps. Low-quality spots and spots not covered by tissue were filtered if <250 genes were detected per spot. For the global sample analysis, samples were merged and normalized using SCTransform and scaled85. For individual stRNA-seq samples, spots were log-normalized. Spots were scored with gene signatures using the AddModuleScore function from the Seurat R package. Spatial differences in VI cluster expression between pathology regions were assessed with a generalized binomial linear mixed effect model of VI cluster enrichment score as a function of spot pathology annotation with sample as a random effect. Pathology regions with <200 spots annotated were removed prior to analysis, and all regions were downsampled to 200 spots each. Normal lung was used as the reference factor level. In order to bin spots into distal VI+ and proximal VI+, image annotation screening on the invaded vessel pathology annotations was performed with the plotSurfaceIAS() function from the SPATA286 R package with a distance and bin width of 1 mm.

Spot deconvolution with scRNA-seq data
Nuclear segmentation of high-resolution H&E images of tissue used for stRNA-seq analysis was performed with the StarDist38,39 algorithm using default settings with a prob_thresh set to 0.3. Per-spot cell number estimates were generated within Squidpy87. Spot deconvolution was performed with the CytoSPACE36 algorithm using the per-spot cell numbers as input and the Salcher stage I LUAD atlas37 as the reference scRNA-seq dataset. The Salcher atlas was downloaded from https://zenodo.org/records/6411868 and was filtered to stage I LUAD samples only. Spatial differences in cell type predictions between pathology regions were assessed with a generalized binomial linear mixed effect model of cell type proportion as a function of pathology annotation with sample as a random effect, after averaging per-spot cell type proportions within a given region annotated as the same pathology. Pathology regions with <200 spots annotated were removed prior to analysis, and all regions were downsampled to 200 spots each. Normal lung was used as the reference factor level.

Deconvolution of bulk RNA-seq data
Deconvolution of the bulk RNA-seq discovery cohort was performed with CIBERSORTx88 using the Salcher stage I LUAD atlas as input to the signature genes matrix generation step. The signature genes matrix was created by taking the average expression for each cell type per sample in the atlas.

Spatially weighted correlation analysis
Spatial transcriptomic spots cannot be assumed to be spatially independent; therefore, spatially weighted regression analysis is used to calculate the correlation between gene signatures. Unlike simple linear regression, which produces a single coefficient per variable, spatially weighted regression allows for coefficients to change locally with tissue coordinates. Data points closer to the coordinates at a given location within a set window (defined by the bandwidth) are given more weight within the model. This window can then slide over the whole tissue to calculate local correlation coefficients per-spot between variables. First, we used Seurat’s AddModuleScore to generate per-spot expression scores for our VI gene cluster signatures and cell type signatures from the Salcher atlas. Cell types that did not have at least one predicted cell by deconvolution across at least 20% of stRNA-seq capture areas were considered too sparse and were filtered out. Cell type expression signatures were generated for each of the remaining cell types annotated in the atlas using Seurat’s FindAllMarkers function with default settings. We selected the top 50 most significantly differentially expressed genes for each cell type. Spatially weighted regression was performed for each sample using the gwss function from the GWmodel89 package 2.3.1 with default settings and a bandwidth set to five. The mean Spearman’s rho of each correlation pair across all tissue spots was then taken for each sample and finally averaged across all samples.

In situ hybridization and immunohistochemistry
Chromogenic immunohistochemistry (IHC) and RNA in situ hybridization (ISH) were performed on a Ventana Discovery ULTRA (Roche Diagnostics). IHC staining for Thy1/CD90 was performed with a rabbit monoclonal antibody (D3V8A, #13801, Cell Signaling) at a 1:100 dilution after heat-induced epitope retrieval in alkaline buffer. IHC for COL1A1 was performed with a rabbit monoclonal antibody (E8F4L, #72026, Cell Signaling) at 1:200 dilution after heat-induced epitope retrieval in alkaline buffer. Primary antibody detection was performed using Discovery Hapten Detection with Discovery Anti-Rabbit HQ (#760-4815, Ventana) and visualized with the DAB chromogen (#760-124, Ventana). RNA ISH for COL1A1 was performed using RNAscope using 2.5 VS probe-HS-COL1A1 (#401899, Advanced Cell Diagnostics) using both protease and heat-induced epitope retrieval as per the manufacturer’s guidelines. Multiplex IHC was performed using the Ventana Discovery ULTRA platform using a rabbit polyclonal antibody against POSTN (HPA012306, Sigma Aldrich) at a 1:100 dilution and the antibody against COL1A1 at 1:200 dilution after heat-induced epitope retrieval for 24 min at 95°C using alkaline buffer (CC1, #950-124, Ventana). Primary antibody detection was performed using Discovery Hapten Detection with Discovery Anti-Rabbit HQ (#760-4815, Ventana), and visualized with purple (#760-229, Ventana) and yellow (#760-239, Ventana) chromogens for COL1A1 and POSTN, respectively. Primary antibody inactivation for the multiplex reactions was carried out using temperature-induced denaturation for 8 min at 95°C with CC2 buffer (#950-123, Ventana).

Spatial alignment of ISH/IHC data to stRNA-seq data
Individual ISH or IHC images were registered separately to the matching adjacent original stRNA-seq H&E image via affine transformation with annotated landmarks using the STalign90 python tool. The same affine transformation was used to align the stRNA-seq spot coordinate grid to the ISH or IHC image. Affine-aligned and filtered stRNA-seq spots were then imported into QuPath v0.5.191. Only spots that were present over tissue in both images and passed QC in the stRNA-seq pipeline were included. QuPath positive cell detection was used to quantify the number of ISH and IHC positive cells at defined staining thresholds to maximize signal-to-noise ratio. For correlation analyses between ISH/IHC and stRNA-seq, we calculated the mean IHC or ISH positive cell fraction, smoothed by k-nearest neighbors, with k=7 (representing a single spot and its immediate neighbor spots). Because the images were not exact serial sections, this approach was applied to smooth minor differences in spot alignment.

Dual IHC of POSTN/COL1A1 data analysis
QuPath v0.5.1 positive cell detection was used to quantify the single positive and double positive cells across all samples. Color convolution vectors were defined from representative single-stained IHC slides for POSTN and COL1A1, and identical intensity thresholds were applied to all images across the dataset. Pathologist-annotated regions of non-desmoplastic stroma were used for the comparison of VI+ and VI− tumors. For both non-desmoplastic stroma and desmoplastic stroma analyses, tumor nests and lymphoid aggregates within these regions were manually excluded from total cell counts.

VI predictor derivation
A nested cross-validation pipeline (Fig. 4A) was used to evaluate the performance of different machine learning models using only the discovery cohort samples and to prevent overfitting model hyperparameters. To accomplish this, samples from the discovery cohort were first divided into outer train and test sets using a 70/30 split 100 times.

Feature selection. Genes with VI-associated expression were derived within each train split as described above for the analysis of the entire discovery cohort: this included both selecting genes with FDR<0.01 and dividing the genes into four clusters using hierarchical clustering. The size of the gene set used for generating a VI predictor in each fold of the cross-validation was set to 48. To automate feature selection and to ensure that diverse VI-associated gene expression patterns were evaluated, genes were selected within each fold of the cross-validation such that the proportion of genes from each of the four clusters was the same as the proportion of genes in that cluster among all the genes with FDR<0.01 in that cross–validation fold. The gene set size of 48 was arbitrarily selected to allow for a 10-fold reduction of the size of the original signature in the full discovery cohort. We observed no dependence of gene set size on model performance from 2.5% up to 40% of the signature, which included median gene set sizes from 12 genes to 192 genes.

Prediction model selection. The outer train split was then divided into inner train and inner test folds using 5-fold cross-validation. This was to allow for evaluating predefined model hyperparameters (hyperparameter search) without learning information about the test data (the outer test fold) that would be used to inform final model selection. Model selection was based on comparing methods using this 5-fold cross-validation within the inner fold of the train split. This was done using the AutoML interface from the h2o.ai92 package v.3.40.0.4 and compared models including generalized linear model (GLM), gradient boosting machine (GBM), extreme gradient boosting (XGBoost), and distributed random forest (DRF). The best performing model and hyperparameter combination in the inner test fold was selected to develop a predictor on the outer train split and then applied to the held out outer test set. This was repeated over 100 cross–validation iterations using the discovery cohort. The final prediction model was selected as the model with the highest mean AUROC across all cross-validation iterations. This prediction model method and its associated hyperparameters, together with the feature selection step, were used with the gene expression data from the entire discovery cohort to develop a final 48-gene VI predictor using a binomial logit GLM utilizing ridge regression.

VI predictor validation
Genes in the validation cohort were subsetted to those in the filtered discovery cohort. The mean and variance of each gene’s log-transformed CPM in the validation cohort were adjusted to match the discovery cohort using reference ComBat93 with the discovery cohort as the reference batch. VI predictor scores were generated using the final 48-gene predictor, and performance was evaluated by AUROC. For evaluating the VI predictor scores by histologic subgroup, samples in the validation cohort were filtered based upon the percentage of annotated histopathologic growth patterns.

Survival analysis
Univariate and multivariate Cox proportional hazards regression was performed using the survival R package.

Additional datasets for survival analysis
Pre-processed RNA-seq and matching clinical data were downloaded using a Genomic Data Commons (GDC) query for TCGA LUAD using the TCGAbiolinks94 R package and filtered to 127 stage IA samples with data on overall survival. Pre-processed RNA-seq and matching clinical data from the Uppsala NSCLC cohort43 were downloaded using the Gene Expression Omnibus (GSE81089) and filtered to 44 stage IA LUAD samples with overall survival data. Because the novel histopathology grading system is based upon the 8th edition of the TNM stage I samples, and our predictor was derived in TNM 8th edition stage I samples, we focused our analysis on this subset. Therefore, since tumor size was not available for either cohort, we excluded stage IB samples because these included tumors up to 5 cm in previous TNM editions and may have been recommended for adjuvant therapy. Genes in both datasets were subsetted to those in the filtered discovery cohort, and the mean and variance of log-transformed CPM were adjusted to match the discovery cohort using reference ComBat with the discovery cohort as the reference batch.

TRACERx intra-tumor heterogeneity analysis
Processed RNA-seq data and clinical annotations from the TRACERx NSCLC cohort of multi-region tumor sampling data were downloaded from https://doi.org/10.5281/zenodo.7683605 and https://doi.org/10.5281/zenodo.760338632. Samples were filtered to TNM 7th edition stage I LUAD tumors ≤4cm with data from at least two regions. Genes were subsetted to those in the filtered discovery cohort, and the mean and variance of log-transformed CPM were adjusted to match the discovery data using reference ComBat with the discovery cohort as the reference batch. VI predictions were generated as described above. To determine intratumor heterogeneity in predictor scores, we calculated the absolute difference in VI predictor score between regions of the same tumor, when all tumors (n=63) were randomly downsampled to two regions each (n=126 regions total). For inter-tumor heterogeneity, the absolute difference was calculated between two regions that were randomly sampled from two different tumors. This was repeated 63 times to make the number of observations equivalent to the intratumor heterogeneity analysis. Spearman correlation between predictor scores from regions of the same tumor was also calculated. Predictor gene enrichment among all genes ranked by region-region correlation was determined using GSEA.

Bayesian estimation of predictive values from biopsy cohort
To estimate the diagnostic performance of the VI predictor in the biopsy cohort, we applied the decision threshold (Youden point) derived from the discovery cohort as a pre-specified cutoff to the pre-surgical biopsy samples. Based on this threshold, we tabulated the number of true positives (TP), false negatives (FN), true negatives (TN), and false positives (FP) among 12 VI+ and 12 VI− biopsy specimens. To model in diagnostic accuracy, we adopted a Bayesian framework in which sensitivity and specificity were treated as random variables with uniform Beta(1,1) priors. Posterior distributions for sensitivity and specificity were therefore Beta(1+TP, 1+FN) and Beta(1+TN, 1+FP), respectively. We generated 10,000 random draws from each posterior to obtain paired samples of sensitivity and specificity that reflect estimation uncertainty. For each simulated pair, positive and negative predictive values (PPV and NPV) are calculated at varying assumed prevalences of VI+ stage I LUAD (15, 20, 25, 30, and 35%) using Bayes’ theorem. For each prevalence level, the median and 2.5th–97.5th percentiles of the simulated distributions were reported as the point estimate and 95% credible interval, respectively.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Steiner et al report a new set of genes that correlate specifically with vascular invasion in their cohort of patients with stage 1 LUAD. They have found four clusters of genes associated with vascular invasion. They have then explored some other characteristics of these genes, such as their site of expression within the tumors in relation to cellular content and their robustness in regards to intra-tumor heterogeneity, using stRNAseq on 15 samples. They have then used a second 60 patient cohort to validate their findings.\nThe main issues with this study are as follows:\n1-The clinical utility of this set of genes as biomarkers for specific management of patients with stage 1 LUAD: As mentioned by the authors, the main management decision for this group of patients is the decision regarding lobectomy vs segmentectomy and two recent trails have failed to show inferiority of segmentectomy at this decision fork. This decision is made before surgical resection and possibly only based on a biopsy and imaging rather than sampling from the resected tumor. Therefore, a critical missing part from this study is validation of this gene signature in biopsies, rather than resected tumors.\nThe other clinically relevant point is if and how this gene expression biomarker adds to the universally used TNM staging, rather than the WHO grading of pathology. For instance, one expects to find a statistical analysis showing that adding vascular invasion to TNM staging has had an impact on the patients’ survival and that this factor has shown an independent predictive value in a multivariable analysis, as shown in Ma et al, 2023. It is of note that even in the latter study, which is one of the main reference points for this paper, vascular invasion was only one of the factors that showed predictive value and other factors such as tumor size, clean margins, and suboptimal nodal sampling had similar or lower p values for predicting survival\n2-Clinical feasibility in a “real world” setting: The proposed 48-gene set is derived from bulk-seq in 100 samples from the firsts cohort and then validated (presumably) using bulk-seq data on 60 tumor samples from the second cohort. Regarding the limitations of bulk seq performance and analysis in clinical settings, how would these 48 gene perform, if measured by more widely available gene probes/qPCR?\n3-Pathophysiological significance: Even though the authors go through extensive description of gene expression patterns based on stRNA-seq, the paper does not contain simple validations such as staining for a group of the high-ranking genes in clinical samples to place them within known cell types and possibly imply a function. The best evidence supporting a biomarker is biological/functional validity that will come out of mechanistic experiments, and this kind of evidence is an implied prerequisite for a high impact publication.\nThe following is a critique of the specific sections:\nThe main finding in the first part of the study is that among the three pathological patterns defined by the authors, vascular invasion has the highest hazard ratio for RFS. While this analysis is valid internally and a multi-variable analysis shows its independent predictive value, the TNM stage of the tumors is missing from this table.\nThe study on the correlation with the histopathological feature is entirely based on stRNA-seq on 15 tumor samples that are mixed stage 1 and 2. Authors have tried to study various issues such as intratumor heterogeneity and structural correlates based on these samples. While interesting, the value and general applicability of statistical reasoning based on findings within these samples are questionable.\nThe same critique applies to the findings on association of gene clusters with specific cell types. As an acceptable validation step for these correlative hypotheses, one would expect a simple staining showing a specific gene with a known function within a specific cell type, possibly together with a mechanistic experiment that captures at least one aspect of that gene function. In the absence of such data, statements such as “ We found a trend toward a higher proportion of peribronchial fibroblasts belonging to VI+ samples (Fig. 3b), consistent with the tissue remodeling biology of VI cluster 2 described above. VI+ samples also had a trend toward a higher proportion of plasma cells and B cells.” remain enigmatic and of questionable value.\nIn the derivation study: was a multivariable analysis, similar to the one performed for VI, done for the 48-genre set? Wouldn’t such analysis be required as the main message of this paper, ie showing if altered expression (levels defined for each gene) of a 48 gene set has an independent predictive value for tumor recurrence or survival.\nThe point about the distinction between vascular and lymphatic invasion is interesting and gene expression patterns are hinting at possible biological differences. However, to the reviewer, none of these findings prove or disprove whether lymphatic invasion is a stage in the continuum, or a biologically different and unique process with clinical implications.\nIn the validation set, the 48 gene signature is applied to 60 samples, only 17 of which had VI. Would a power calculation for this analysis be appropriate? I cannot find the data for the sentence “It was associated with histologic grade and 7-year recurrence free survival (HR=1.98, 95% CI 1.09–3.60, Cox regression p=0.024) (Extended Data Fig. 6b and Fig. 4f).”\nIn discussion, some of the pathophysiological explanations are not convincing or need more explanation. For instance, why would a group of “cell cycle” or “tumor hypoxia” genes specifically contribute to vascular invasion. Why does cluster 2 (EMT/Angiogenesis) correlate with the presence of CAFs? What is the value of necrosis? How would a gene signature of vascular invasion have prognostic value in VI- tumors?\nOn stylistic issues, the results section is filled with technical details, and this overrepresentation of details detracts from the main message of the paper and its proposed value as a biomarker of invasiveness. Also, some of the premises of the study are not solid. For instance, authors mention that in the Hattori et al study (2024) “…post-hoc analyses are beginning to show subsets of patients for which lobectomy might decrease the rate of recurrence.” The mentioned study showed that lobar resection was associated with higher locoregional recurrence among the full cohort, but survival was age and sex- dependent. Authors mention that gene signature for vascular invasion has not been studied. A number of studies in the literature including https://doi.org/10.1016/j.ajpath.2010.10.040, https://doi.org/10.1371/journal.pone.0216847, and https://doi.org/10.1038/s41416-019-0486-6, had visited this entity in cancer.", "reviewer_2": "This manuscript reports a lung adenocarcinoma transcriptional signature that is associated with the pathological detection of vascular invasion in resected stage 1 tumors.\nThe association was confirmed in an independent cohort of 60 patients with AUROC of 0.86. The signature was associated with recurrence free survival and was detectable in regions without invasion in tumors that were annotated as harboring vascular invasion.\nThe study design is sound and the observations are robust, although not conceptually novel.\n1. The determination of vascular invasion was by a single thoracic pathologist. Were these cases previously read by other clinical pathologists and what was the agreement of determination of vascular invasion in cases read by more than one pathologist? The authors should address inter-reader variability concerns with determination of this pathological finding.\n2. The introduction emphasizes the potential clinical significance of detection of vascular invasion and suggests that the finding may support use of adjuvant therapy or lobectomy. This is confusing and should be clarified. Further, it is possible that the signature impact on recurrence free survival may be independent of performance of lobectomy or administration of adjuvant therapy. Addressing the potential clinical utility aspects of the proposed classifier and how it can be assessed will be helpful.\n3. Figure 3 shows imputed cell type associations with the vascular invasion signature using spatial transcriptomic and bulk sequencing data. It will be important to validate key gene expression signatures in specific cell types using immunohistochemistry.", "reviewer_3": "In their manuscript “Identification of a gene expression signature of vascular invasion and recurrence in stage I lung adenocarcinoma via bulk and spatial transcriptomics”, Steiner et al. developed a transcriptomic signature of microscopic vascular invasion (VI) in stage I lung adenocarcinoma based on a cohort of 163 tumors. The VI-associated expression patterns clustered into four subclasses linked to unique biological processes. The authors utilized spatial transcriptomics to describe the spatial organization of the four gene clusters. The authors additionally developed a predictor of VI in an independent cohort and evaluated its performance in multiple regions of the same tumor to show the predictor is largely unaffected by intra-tumor heterogeneity. The novel transcriptional signature for vascular invasion in early-stage lung adenocarcinoma presented in this work may enable more precise diagnosis and treatment of patients.\nMy specific comments which may help to improve the manuscript are (as line numbers are not available, I am referencing the figure panels to identify the relevant section of the manuscript):\nExt. Data Fig. 1a-f: It would be helpful to add outcomes that include both VI and LI classification given that LI is used as a covariate for some of the differential expression analysis.\nExt. Data Fig. 1j: In the legend, “only” is cut-off.\nExt. Data Fig. 1j-k: Is the site of recurrence related to time to recurrence?\nFig. 1b: It may be helpful to have VI/NST/LMP as a single column annotation as they are mutually exclusive and based on the same grading system.\nFig. 1b: The differential expression results and cluster assignments should be provided as supplementary tables or in a data repository.\nFig. 1b: There are two VI samples that are clustering with LMP samples. Are there any biological or technical factors that may explain this grouping?\nFig. 1b: The NST samples are clustering with both VI and LST samples. Can you comment on how to best interpret the two VI classifications?\nFig. 1b: Cluster 1 and, to a lesser degree, cluster 3 appear to be separating the VI+ samples into two subsets. Is there any merit to treating those as independent groups for downstream analysis?\nExt. Data Fig. 3a: It may be helpful to add a brief explanation of the biological significance of adding LI as a covariate.\nExt. Data Fig. 3b: In the main text, it is not clear which specific comparison was performed to get these genes.\nExt. Data Fig. 4d: Are all genes in all samples shown? Is the correlation similar for just the VI gene signature? Is the correlation consistent for all samples?\nFig. 3a: It would be helpful to have the cell types order be consistent with Fig. 3b.\nFig. 3b: Are all the cell types listed in the legend present? For this visualization, it may be beneficial to collapse all the cell types with negligible proportions.\nFig. 3b: It would be useful to additionally show the cell type frequencies across the pathology annotations.\nFig. 3b: Are there other cell types in addition to peribronchial fibroblasts with significantly altered abundances between the two groups?\nFig. 3c: It would be interesting to see the stRNA-seq proportions in addition to the bulk RNA-seq ones.\nFig. 3d: If you check the enrichment of the VI cluster signatures across all the cell types (analogous to Ext. Data Fig. 5b), does that yield different associations?\nFig. 4a: The gene selection process includes dividing the genes into four clusters. Can you clarify the reasoning for the step?\nFig. 4b: Does increasing or decreasing the gene set size improve the performance of the predictor?\nFig. 5a: The figure legend states “VI+ tumors that did not contain VI foci in the capture area were considered as proximal VI+”. Would the tumors without VI foci be considered distal?\nThe URL provided in the Reviewer Assessment form (https://github.com/dst53/VI_manuscript) is different from the URL provided in the Code and Software Submission Checklist and the Reporting Summary (https://github.com/dst53/Steiner_NatComm). Neither one is accessible."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nDear reviewers,\nThank you for your insightful comments. We have made many changes to the manuscript\nunder your guidance and now report an important advance in our study. To demonstrate\nreal-world feasibility and clarify the clinical application of our method, we profiled 36 new\nstage I LUAD FFPE pre-surgical biopsy specimens (24 biopsies, 12 matched resections)\nfrom an independent institution. Our VI predictor measured in biopsies achieved an\nAUROC of 0.90 for predicting VI at resection and showed strong concordance with\nmatched resected tumors. These results underscore the potential of our predictor to help\ntailor pre- and peri-operative care and extend our biomarker’s applicability beyond\nresected tumors. Collectively, this additional data brings our manuscript from the initial\ndiscovery of a robust gene expression signature linked to VI to the threshold of clinical\nimplementation. In response to the reviewers’ request for orthogonal validation of key\nfindings from the spatial transcriptomics studies, we have now performed validation of a\nkey VI cluster 2 gene, COL1A1, by RNAscope ISH and IHC on serial sections adjacent to a\nVI+ stRNA-seq sample. We find that COL1A1 signal (both RNA and protein) correlates\ntightly with our cluster 2 enrichment scores and localizes to the expected desmoplastic\nstroma and myofibroblastic CAFs. These data validate our previous observations from the\nstRNA-seq data and clarify the link between VI cluster 2, CAFs, and tissue remodeling,\nsuggesting that there is a strong tumor microenvironment component to our VI gene\nexpression signature. Finally, we addressed all other reviewer comments in detail. We\nclarified our clustering methodology and included new data to support our discovery of\nfour main gene expression clusters linked to vascular invasion in response to Reviewer #3.\nWe have refined our discussion in response to Reviewer #1 to better contextualize the\nroles of cell-cycle, hypoxia, EMT/angiogenesis, CAFs, and necrosis in vascular invasion.\nWe have also included new multivariate survival analyses to show that our VI predictor is\nindependent of TNM stage, surgical type, and adjuvant therapy. We also note a minor\nchange to the original reported n of cases that passed QC in the validation cohort (now\nn=59, previously n=60). This was reported in error and has now been fixed, with no changes\nto the reported performance results. We provide all changes to the manuscript as\nhighlighted text in the updated file and elaborate on each reviewer comment below."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In response to my comments authors have performed a study on 24 selected FFPE samples and IHC and ISH for col1a1 on several slides from one sample. The following are my comments on these and other modifications:\n1-clinical utility:\nThe addition of the FFPE biopsies to the end of the paper increases the clinical applicability of the paper. However, unfortunately it falls short of the main purpose of my critique which was to show the actual clinical utility of this approach. Testing the results of the bulk seq on 12 VI+ and 12VI- preselected samples does not mean that this approach will work on FFPE biopsies. This groups of samples do not have the variation seen in a regular sample pool and the AUROC derived from this set is of questionable applicability for clinical practice. Also, the placement of this section as a proof of concept defies the main purpose of this critique; the fact that the main utility of this finding, if validated, would be iy diagnostic utility for FFPE samples.\n-The addition to the TNM staging:\nThe statistical analyses showing that VI remains an independent predictor of RFS with TNM stage added is satisfactory. (fig1i extended) However,\n-I cannot find the clinical characteristics of this 183 tumors that are “a subset” of previously published data (Yambayev et al). I can not find 1823 samples in that paper.\n-Also: what percentage of the data is missing for each variable?\n-About the age and pack year in these samples. Since these are continuous variables why are there no CIs?\n-Extrinsic factors: Even though author’s response about focusing on VI as an intrinsic factor (vs extrinsic factors such as clear margins and suboptimal sampling) is valid, the extrinsic features will still have importance in determining the predictive power of a biomarker in real world and should at least be mentioned in the discussion.\n2-Clinical feasibility:\nThe authors in their response mention that in a real-world scenario the remainder of an FFPE after pathologist review will probably be used for bulk-seq, stressing the importance of the FFPE sample validation. Yet perhaps the weakest point of this manuscript at present is the analysis of the FFPE samples (24 preselected VI+ and VI- samples). Increasing the number of the FFPE samples and testing this approach on a true representative sample pool (rather than as an add-on proof of concept piece) will show the basic utility of this study.\nThe comment about the use of RNA-seq in current clinical practice: the sentence “whole transcriptome RNA-seq test….fine-needle aspiration FFPE material used … is now standard of care for preoperative diagnosis…” should be moderated. The use of molecular probes for indeterminate thyroid nodules refers to a small group of samples assessed by FNA, and mRNA sequencing is only one of the three approaches used for these samples: mutational analysis, mRNA sequencing, and miRNA expression analysis. Authors are presenting Bulk-seq as a method for general classification of all early stage NSCLC samples, which is a very different premise.\nAlso, on the disparities between RNA-seq and qRT-PCR please refer to:\n-Int. J. Mol. Sci. 2021, 22(17), 9349; https://doi.org/10.3390/ijms22179349\nDespite the extensive discussion on the advantages of sequencing by some authors and commercial entities, I think we all know that the access to RNAseq is much more limited and technical caveats are much more prevalent, compared to hybridization probes or qPCR. These methods are currently used for biomarker analysis in many clinical settings and for various cancers, including breast and colon cancer. Therefore, validation by qPCR or RNA probes for some genes, which is routinely done for RNAseq studies is reasonable and adds to the possible clinical utility.\n3-pthophysiological significance:\nThis is one of the main critiques regarding the appropriateness of this paper for Nature communications. Despite the extensive work and the multitude of statistical analyses, the conclusions in the manuscript are based on correlative data and deficient in mechanistic insights, which is a prerequisite for publication in journals like Nature Comm. The involvement of CAFs in vascular invasion is known. However, the main mechanisms mentioned for this contribution are increased angiogenesis, expression of ECM remodeling enzymes, immune response, and metabolic programming of tumor cells.\n(https://doi.org/10.3389/fimmu.2025.1582532)\nAuthors have picked Col1a1 for validation. This gene is expressed as part of a fibrotic tissue response to injury/invasion. Mechanistic questions remain. For instance: Does col1a1 contribute to the process of invasion? If so, how? Or is it just a consequence? Does it identify “activated” CAFs? Is there any mechanism downstream of col1a1 that can be shown to enhance any aspect of vascular invasion? In a mechanistic paper these questions need to be answered by experiments, not just correlation and citing literature.\n-Critique on association with a specific cell type: staining for one of the genes.\nAuthors have performed IHC and ISH staining for col1a1, on multiple sections from one sample. Question arises, why could this not be done on more than one sample, as is regularly done for the histopathological validation of a gene involved in a pathological process? Staining even one or two sections from 10-20 samples for col1A1 and presenting the histopathological images that show the expression of this gene in myofibroblasts within the conventional “vascular invasion” sites would be appropriate for this response. The multiple statistical analyses show that within sample 4 the expression of col1A1 coincides with the expression of the other genes in cluster 2 in that sample. This kind of analysis adds to the rigor but does not provide an adequate validation and link between this new finding and the currently used histopathological studies that are the basis of the statistical analysis in extended figure 1i.\n-Multivariable analysis on the derivation study: Extended Figure 7c is adequate and a critical piece of data. The same questions about the CI for age and pack year apply and need to be answered. If the latter questions are adequately answered, I suggest moving this data to the main figures.\n-On the distinction between LI and VI : the changes are adequate. I would add that LI was also an independent predictor in extended fig.1i and would show the gene signature associated with it as a variable in extended figure 7c.\n-The explanations for tumor hypoxia and cell cycle genes are adequate.\n-There are still a multitude of technical details in result sections that take away from the main message in each paragraph.\nThe literature references are adequate.\nAdditionally:\nThe manuscript still has several stylistic and language issues and would benefit from cleaning up.\nSome examples:\n-What is the meaning of the sentence “ tumors were graded…. Or VI, is associated with outcome than the World…”\n-The text still has redundant and hard to understand details that do not clarify the significance: For instance:\n“ based on our analysis from consensus clustering the delta area under the cumulative….” This sentence and sentences alike should be modified/simplified or deleted.\n-Where can we find the list of the genes for the four clusters described in the second results paragraph?\n-In discussion paragraph 7:\n“High VI predictor scores were associated not only with… VI- tumors were shedding cells into the circulation…”\nThis sentence is confusing and needs clarification.", "reviewer_2": "The authors have responded to the main comments that are listed below.\n1. The determination of vascular invasion was by a single thoracic pathologist. Were these cases previously read by other clinical pathologists and what was the agreement of determination of vascular invasion in cases read by more than one pathologist? The authors should address inter-reader variability concerns with determination of this pathological finding.\nThe authors’ response is complete.\n2. The introduction emphasizes the potential clinical significance of detection of vascular invasion and suggests that the finding may support use of adjuvant therapy or lobectomy. This is confusing and should be clarified. Further, it is possible that the signature impact on recurrence free survival may be independent of performance of lobectomy or administration of adjuvant therapy. Addressing the potential clinical utility aspects of the proposed classifier and how it can be assessed will be helpful.\nThe discussion addresses these points. To further address potential clinical utility, the authors added analysis of a cohort of presurgical FFPE lung biopsy specimens acquired from 24 subjects in Figure 5f. The AUROC is shown. To better understand the performance of this proposed biomarker, it will be important to show the sensitivity, specificity, PPV and NPV of the predictor. In addition, it is important to carefully describe the specimens from which the predictor was evaluated. Were these consecutively acquired specimens or were they selected? Were the specimens derived from cell block or forceps biopsy? What was the evidence for vascular invasion on the biopsy specimens used for the analysis, if from forceps biopsy?\n3. Figure 3 shows imputed cell type associations with the vascular invasion signature using spatial transcriptomic and bulk sequencing data. It will be important to validate key gene expression signatures in specific cell types using immunohistochemistry.\nIn extended figure 6, The authors show analysis of COL1A1IHC and ISH in one sample. This does not address validation of key V1 biomarkers cell specific expression in lung adenocarcinoma specimens with and without vascular invasion. These data will be important to demonstrate the generalizability of the reported observations.", "reviewer_3": "The authors addressed my comments and introduced satisfactory corrections."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author): expertise in lung cancer TME\nSteiner et al report a new set of genes that correlate specifically with vascular invasion in\ntheir cohort of patients with stage 1 LUAD. They have found four clusters of genes\nassociated with vascular invasion. They have then explored some other characteristics of\nthese genes, such as their site of expression within the tumors in relation to cellular\ncontent and their robustness in regards to intra-tumor heterogeneity, using stRNAseq on\n15 samples. They have then used a second 60 patient cohort to validate their findings.\nThe main issues with this study are as follows:\n1-The clinical utility of this set of genes as biomarkers for specific management of patients\nwith stage 1 LUAD: As mentioned by the authors, the main management decision for this\ngroup of patients is the decision regarding lobectomy vs segmentectomy and two recent\ntrails have failed to show inferiority of segmentectomy at this decision fork. This decision is\nmade before surgical resection and possibly only based on a biopsy and imaging rather\nthan sampling from the resected tumor. Therefore, a critical missing part from this study is\nvalidation of this gene signature in biopsies, rather than resected tumors.\nWe thank the reviewer for calling attention to the potential clinical applicability of our gene\nsignature, and we agree that the ability to identify vascular invasive tumors prior to the\ndecision of lobectomy vs segmentectomy is a major potential application of our gene\nsignature. To demonstrate the potential utility of the VI gene signature in this setting, we\nhave now validated the gene signature by RNA-sequencing of pre-surgical biopsies and\nmatched resected tumors obtained from an independent cohort of stage I LUAD patients.\nIn total, we performed RNA-seq on 36 new samples, which included 24 pre-surgical tumor\nbiopsies and 12 from matched subsequently resected tumors. In this cohort, our VI\npredictor measured in presurgical tumor biopsies showed strong performance in\npredicting the detection of vascular invasion in the resected tumor (n=12 VI-, n=12 VI+)\nwith an AUROC of 0.90 (p=4.96x10-4). Furthermore, we observed strong correlation in VI\npredictor score between the pre-surgical tumor biopsies and the matched resected\ntumors (r=0.77, p=0.0053). These samples were evaluated retrospectively from residual\nclinical samples from patients who received their care at Inova Schar Cancer Institute\n(extending our biomarker validation to another institution) and contain no overlap with our\npreviously reported discovery or validation cohorts. They are also formalin-fixed paraffin\nembedded (FFPE) pre-surgical tumor biopsies, which is important for the clinical feasibility\nof our biomarker and ensures its compatibility with existing clinical workflows. After a\npathologist has examined H&E slides of the pre-surgical biopsy, remaining FFPE biopsy\nmaterial can be sent for molecular testing to generate a prediction of VI risk score,\nenabling timely pre-surgical decision making. Collectively, these results suggest that our VI\npredictor may have strong clinical utility in guiding pre-surgical management decisions for\nstage I LUAD patients. We have now updated the manuscript to reflect these important\nnew results, including the figures (see Fig. 5f and 5g), tables (see Table S4), abstract,\nmethods, results, and discussion sections.\nThe other clinically relevant point is if and how this gene expression biomarker adds to the\nuniversally used TNM staging, rather than the WHO grading of pathology.\nWe thank the reviewer for raising this relevant point. We agree that it is important to\nevaluate the gene expression biomarker in the context of TNM staging. While our validation\ncohorts are presently underpowered to stratify patient RFS by both TNM stage IA/IB + VI\npredictor, we have added additional analysis in the TRACERx cohort showing that our VI\npredictor is an independent predictor of RFS across TNM stage I and stage II LUADs. We\nhave now included this result as a new panel in Extended Data Fig. 7 (previously Extended\nData Fig. 6). We acknowledge the limitation of being underpowered to show that the VI\npredictor is useful for identifying poor prognosis clinical stage IA or good prognosis clinical\nstage IB and have added this to the discussion. We are currently running an ongoing\nprospective clinical trial that will further explore the added benefit of VI predictor staging in\nlarger cohorts. We have now added the following text to the results section of the\nmanuscript: “Higher VI predictor scores were significantly associated with decreased 5-\nyear RFS in TRACERx (Fig. 4f; HR=1.69, 95% CI 1.04–2.75, Cox regression p=0.036). The\npredictor was also significantly associated with 5-year RFS independently of TNM stage I\nand II (Extended Data Fig. 7g; HR=1.46, 95% CI 1.07-2.0, Cox regression p=0.018).”\nFor instance, one expects to find a statistical analysis showing that adding vascular\ninvasion to TNM staging has had an impact on the patients’ survival and that this factor has\nshown an independent predictive value in a multivariable analysis, as shown in Ma et al,\n2023.\nWe thank the reviewer for pointing out this omission. We have now included TNM stage as\nan additional covariate in our multivariate Cox regression analysis and updated Extended\nData Fig. 1i and the results accordingly.\nIt is of note that even in the latter study, which is one of the main reference points for this\npaper, vascular invasion was only one of the factors that showed predictive value and\nother factors such as tumor size, clean margins, and suboptimal nodal sampling had\nsimilar or lower p values for predicting survival.\nWe thank the reviewer for raising this important point. We agree that factors such as tumor\nsize, clean margins, and nodal sampling are all important determinants of prognosis. For\ndecision making about pre-surgical management, we focus on intrinsic features of the\ntumor (such as size – which is already an integral component of TNM staging, and vascular\ninvasion). In contrast, resection margins and the extent of nodal sampling reflect aspects\nof the surgical resection that we hope will be influenced by the preoperatively determined\nintrinsic features of the tumor.\n2-Clinical feasibility in a “real world” setting: The proposed 48-gene set is derived from\nbulk-seq in 100 samples from the firsts cohort and then validated (presumably) using bulk-\nseq data on 60 tumor samples from the second cohort. Regarding the limitations of bulk\nseq performance and analysis in clinical settings, how would these 48 gene perform, if\nmeasured by more widely available gene probes/qPCR?\nWe thank the reviewer for this question about clinical feasibility. Our new results show\nvalidation of our 48-gene set VI predictor in FFPE pre-surgical biopsy material by RNA-seq.\nWe anticipate that after a pathologist reviews H&E slides of the pre-surgical biopsy, the\nremaining FFPE material can be used directly for RNA extraction and sequencing, enabling\ncompatibility with existing clinical workflows. To that end, we are not aware of any\nlimitations of bulk RNA-seq performance and analysis in clinical settings. In fact, RNA-seq\nhas proven clinical use even on small FFPE cytology specimens. For example, Veracyte’s\nAfirma Genomic Sequencing Classifier is a CLIA-certified, whole transcriptome RNA-seq\ntest that is performed on the same fine-needle aspiration FFPE material used for cytology\nreview and is now standard of care for preoperative diagnosis of indeterminate thyroid\nnodules in thousands of patients (PMID: 29799911). This demonstrates that RNA-seq on\nFFPE clinical specimens can integrate into routine practice. While we cannot answer how\nour 48-gene predictor would perform when measured by qPCR, large multicenter efforts\n(SEQC/MAQC-III) have shown that RNA-seq expression measurements correlate tightly\nwith TaqMan qPCR (Pearson r > 0.9 for most genes) (PMID: 25150838). Furthermore, RNA-\nseq has additional advantages over qPCR that are desirable for clinical implementation of\nour VI predictor. RNA-seq based normalization is more robust and means that predictor\nperformance will be less sensitive to variation in any single reference gene used by qPCR.\nRNA-seq also has built-in sample QC metrics, such as mapping rates, and the ability to\nperform orthogonal checks on sample identity. We are nevertheless eager to determine the\ncost effectiveness of measuring the VI predictor using other laboratory methods in future\nstudies.\n3-Pathophysiological significance: Even though the authors go through extensive\ndescription of gene expression patterns based on stRNA-seq, the paper does not contain\nsimple validations such as staining for a group of the high-ranking genes in clinical\nsamples to place them within known cell types and possibly imply a function. The best\nevidence supporting a biomarker is biological/functional validity that will come out of\nmechanistic experiments, and this kind of evidence is an implied prerequisite for a high\nimpact publication.\nWe thank the reviewer for raising this important point. We have addressed the request for\nvalidation via staining in the reviewer’s related comment below about validation by\nstaining.\nThe following is a critique of the specific sections:\nThe main finding in the first part of the study is that among the three pathological patterns\ndefined by the authors, vascular invasion has the highest hazard ratio for RFS. While this\nanalysis is valid internally and a multi-variable analysis shows its independent predictive\nvalue, the TNM stage of the tumors is missing from this table.\nWe thank the reviewer for pointing out this omission. We have now added TNM stage 8th\nedition to both Table S1 and Table S3 and included TNM stage as a covariate in the\nmultivariable analysis presented in Extended Data Fig. 1i, as described above.\nThe study on the correlation with the histopathological feature is entirely based on stRNA-\nseq on 15 tumor samples that are mixed stage 1 and 2. Authors have tried to study various\nissues such as intratumor heterogeneity and structural correlates based on these\nsamples. While interesting, the value and general applicability of statistical reasoning\nbased on findings within these samples are questionable.\nWe thank the reviewer for raising the concern about sample size in our stRNA-seq analysis.\nAlthough or initial dataset comprised 15 sections (8 VI-, 7 VI+), we used a linear mixed-\neffects model with sample ID as a random effect to control for inter-sample variability,\nyielding robust associations (e.g. Cluster 2 with desmoplastic stroma – a finding that we\nhave since validated by IHC and ISH – see response below). To further demonstrate that\nthese findings are not cohort-specific, we validated our results using 11 additional Visium\nsections (6 VI-, 5 VI+; no overlap with our discovery or validation sets). In these new\nsamples, we similarly found that VI gene cluster 2 is strongly associated with desmoplastic\nstroma (p.adj=1.14x10-63). In addition, as before, we observed strong expression of Cluster\n1 and Cluster 3 genes directly in VI foci (p = 9.70x10-7 and 2.63x10-56 respectively). We also\nvalidated that cluster 1 and cluster 3 have increased expression in regions annotated as\nhigh-grade patterns: solid (p.adj= 6.84x10-7; 2.26x10-68), micropapillary (p.adj=1.32x10-3;\n3.94x10-101), and cribriform (p.adj=1.56x10-48; 9.50x10-100). Finally, expression of cluster 4\nwas strongly decreased in VI foci (p.adj=6.60x10-38) and showed positive enrichment in\nlepidic pattern (p.adj=8.68x10-8), normal bronchus (p.adj=1.10x10-5) and adjacent normal-\nappearing lung (p.adj=1.57x10-3). Validating these results in a second independent stRNA-\nseq cohort suggests that our statistical approach studying the association of VI cluster\nexpression by pathology annotation is valid and yields highly reproducible findings.\nBecause the new data used for validation is part of separate ongoing study and cannot yet\nbe made publicly available, we have opted not to include these validation results in the\npresent manuscript. We include them in the figure below in our response to the reviewer\nonly.\nReviewer Response Figure 1. Reproducing our stRNA-seq histopathology findings in an\nadditional cohort a. Pathology annotations present across new spatial transcriptomics\n(stRNA-seq) capture areas that passed QC (n=11). b. Association between expression of\nVI-associated clusters and pathology across all 11 samples. Only results with Bonferroni\nadjusted p values (p.adj) < 0.01 are shown. P.adj values were derived from a generalized\nbinomial linear mixed effect model of VI cluster enrichment score as a function of spot\npathology annotation with sample as a random effect. For visualization purposes the plot\nshows the mean of per-spot enrichment scores within a given region annotated as the\nsame pathology.\nThe same critique applies to the findings on association of gene clusters with specific cell\ntypes. As an acceptable validation step for these correlative hypotheses, one would expect\na simple staining showing a specific gene with a known function within a specific cell type,\npossibly together with a mechanistic experiment that captures at least one aspect of that\ngene function. In the absence of such data, statements such as “ We found a trend toward\na higher proportion of peribronchial fibroblasts belonging to VI+ samples (Fig. 3b),\nconsistent with the tissue remodeling biology of VI cluster 2 described above. VI+ samples\nalso had a trend toward a higher proportion of plasma cells and B cells.” remain enigmatic\nand of questionable value.\nWe thank the reviewer for providing this valuable feedback around validating our cell type\nassociation findings. We have now generated and included the requested data and added\nadditional analyses. Because our main finding about cell types (Figure 3) was the\nassociation between VI gene cluster 2 and peribronchial fibroblasts (or myofibroblasts),\ndesmoplastic stroma, and VI+ tumors, we have now performed in situ hybridization (ISH)\nand immunohistochemistry (IHC) for COL1A1, a key Cluster 2 gene. Below is a summary of\nthe new data and how it clarifies the link between VI cluster 2, peribronchial fibroblasts,\nand tissue remodeling.\nFirst, we previously demonstrated that the peribronchial fibroblasts from the Salcher LUAD\natlas we used for our stRNAseq deconvolution analysis are referred to as myofibroblasts\nby the Hanley et al 2023 LUAD fibroblast scRNAseq dataset. This is important because we\nused the Salcher et al data for deconvolution of our stRNAseq data. The original\nsubmission of the manuscript demonstrates that these differently named cell populations\nare very likely equivalent. As part of the revised manuscript, we have further scored the\nHanley et al data (an independent dataset consisting of 46 independent scRNAseq\nsamples) for expression of the genes in VI gene cluster 2. We observed significantly higher\nexpression of VI gene cluster 2 in Hanley’s myofibroblast cluster compared with alveolar\nand adventitial clusters. Hanley et al reported that the key markers of LUAD myofibroblasts\nare MMP11, POSTN, CTHRC1, COL1A1, and COL3A1, which are all VI cluster 2 genes.\nNext, to validate our stRNA-seq deconvolution results directly, we performed COL1A1\nRNAscope ISH and COL1A1 protein IHC labeling on serial sections close to stRNA-seq\nsample 4, a VI+ tumor with prominent desmoplastic stroma. We selected sample 4 due to\nthe difficulty in obtaining serial sections close enough to the original stRNA-seq section in\nthe other samples in our cohort. After computationally registering the adjacent sections to\nthe original Visium image and creating a ‘virtual’ Visium spot grid on the adjacent sections,\nwe observed a significant correlation across spots between the fraction of ISH RNAscope\nCOL1A1+ expressing cells per spot (quantified by QuPath) and the average VI cluster 2\ngene expression per spot (from spatial transcriptomics) (r=0.53, p < 2.2 x10-16). We also\nobserved a similar significant correlation across spots between average VI cluster 2 gene\nexpression and COL1A1+ protein expressing cell fraction measured by IHC (r=0.47, p < 2.2\nx10-16). We also did IHC staining for another VI cluster 2 gene, THY-1, that was not reported\nby Hanley et al as a key myofibroblast marker. The fraction of THY-1+ protein expressing\ncells measured by IHC was significantly correlated with the fraction of COL1A1+ protein\nexpressing cells measured by IHC (R=0.52, p < 2.2 x10^-16), which further supports the\ncoexpression of VI gene cluster 2 at the protein level. When we examined the correlation\nbetween COL1A1+ expressing cells measured by both IHC or ISH, and the proportion of\nmyofibroblasts predicted by our stRNA-seq deconvolution algorithm, we found that both\nCOL1A1 RNA and protein were significantly more correlated with the fraction of\nperibronchial (ISH r=0.43, p=7.95x10-159; IHC r = 0.36, p = 1.66x10-107 ) than adventitial (ISH\nr=-0.10, p=4.1x10-9; IHC r=0.11, p=2.47x10-10) or alveolar fibroblasts (ISH r=0.14,\np=1.05x10-17; IHC r=-0.05, p=0.001) (Steiger’s p value for difference in COL1A1 correlation\nbetween peribronchial vs adventitial = 7.85x10-17 for IHC and 4.18x10-78 for ISH; Steiger’s p\nvalue for difference in COL1A1 correlation between peribronchial and alveolar = 3.12x10-24\nfor IHC and 5.61x10-101 for ISH). This further demonstrates that cells expressing a key\nmarker of our VI gene cluster 2 show a myofibroblast-like phenotype. Finally, we also\nobserved a significant association between the fraction of COL1A1+ protein expressing\ncells measured by IHC and spots annotated by the pathologist as belonging to\ndesmoplastic stroma (p = 4.17 x 10-16), further suggesting that VI cluster 2 is associated\nwith desmoplastic stroma and tissue remodeling. We have now updated the figures (see\nFig. 3 and Extended Data Fig. 6) results section, methods section detailing these additional\nfindings and our analysis approach.\nWe have added the following text to the results section:\nTo further validate our stRNA-seq deconvolution findings about the relationship between VI\ncluster 2, desmoplastic stroma, and peribronchial fibroblasts, we performed COL1A1\nRNAscope ISH and COL1A1 IHC labeling on serial sections close to stRNA-seq sample 4, a\nVI+ tumor with prominent desmoplastic stroma. We selected COL1A1 as it is a\nrepresentative gene from VI cluster 2 that also showed significantly higher expression in\nmyofibroblasts compared with adventitial or alveolar fibroblasts in Hanley et al (Extended\nData Fig. 6a) After serial image registration to the H&E from the section used for stRNA-\nseq, we observed a significant correlation across spots between the fraction of ISH\nRNAscope COL1A1+ expressing cells per spot and the average VI cluster 2 gene expression\nper spot (r=0.53, p < 2.2x10-16; Extended Data Fig. 6b). We also observed a similar\nsignificant correlation between the average VI cluster 2 gene expression per spot and the\nfraction of COL1A1+ cells per spot labeled by IHC (r=0.47, p < 2.2x10-16; Extended Data Fig.\n6c). We also performed IHC for THY-1, a VI cluster 2 gene that was not reported by Hanley\net al as a key myofibroblast marker. The fraction of THY-1+ expressing cells per spot by IHC\nwas significantly correlated with the fraction of COL1A1+ expressing cells measured by\nIHC, which further supports the coexpression of VI gene cluster 2 at the protein level\n(R=0.52, p < 2.2x10-16; Extended Data Fig. 6d). When we examined the correlation between\nCOL1A1+ cell fractions measured by either IHC or ISH with the fraction of myofibroblasts\npredicted by the CytoSPACE stRNA-seq deconvolution, we found that both were\nsignificantly more correlated with the fraction of peribronchial/myofibroblastic CAFs than\neither adventitial or alveolar fibroblasts (Extended Data Fig. 6e). This further demonstrates\nthat cells expressing COL1A1, a key marker of VI cluster 2, show a myofibroblast-like\nphenotype. Finally, we also observed a significant association between the fraction of\nCOL1A1+ IHC fraction and spots annotated as desmoplastic stroma (Extended Data Fig.\n6f), further suggesting that VI cluster 2 is associated with tissue remodeling.\nWe have added the following text to the methods section:\nIn situ hybridization and immunohistochemistry. Chromogenic immunohistochemistry and\nRNA in situ hybridization were performed on a Ventana Discovery ULTRA (Roche\nDiagnostics). Immunohistochemical staining for Thy1/CD90 was performed with a rabbit\nmonoclonal antibody (D3V8A, #13801, Cell Signaling) at 1:100 dilution after heat-induced\nepitope retrieval in alkaline buffer. IHC for COL1A1 was performed with a rabbit\nmonoclonal antibody (E8F4L, #72026, Cell Signaling) at 1:200 dilution after heat-induced\nepitope retrieval in alkaline buffer. RNA ISH for COL1A1 was performed using RNAscope\nusing 2.5 VS probe-HS-COL1A1 (#401899, Advanced Cell Diagnostics) using both protease\nand heat-induced epitope retrieval as per the manufacturer’s guidelines.\nSpatial alignment of ISH/IHC data to stRNA-seq data. Individual ISH or IHC images were\nregistered separately to the matching adjacent original stRNA-seq H&E image via affine\ntransformation with annotated landmarks using the STalign87 python tool. The same affine\ntransformation was used to align the stRNA-seq spot coordinate grid to the ISH or IHC\nimage. Affine-aligned and filtered stRNA-seq spots were then imported into QuPath\nv0.5.188. Only spots that were present over tissue in both images and passed QC in the\nstRNA-seq pipeline were included. QuPath positive cell detection was used to quantify the\nnumber of ISH and IHC positive cells at defined staining thresholds to maximize signal to\nnoise ratio. For correlation analyses between ISH/IHC and stRNA-seq, we calculated the\nmean IHC or ISH positive cell fraction, smoothed by k-nearest neighbors, with k=7\n(representing a single spot and its immediate neighbor spots). Because the images were\nnot exact serial sections this approach was applied to smooth minor differences in spot\nalignment.\nIn the derivation study: was a multivariable analysis, similar to the one performed for VI,\ndone for the 48-genre set? Wouldn’t such analysis be required as the main message of this\npaper, ie showing if altered expression (levels defined for each gene) of a 48 gene set has\nan independent predictive value for tumor recurrence or survival.\nWe thank the reviewer for pointing out this omission. We have now included a multivariate\nsurvival analysis with the VI predictor in the discovery cohort as an additional figure panel\nin Extended Data Fig. 7 (previously Extended Data Fig. 6). We show here that the VI\npredictor is independently predictive of RFS when including common clinicopathology\nvariables, including TNM stage. We have added the following text to the results section of\nthe manuscript: “The predictor was independently associated with 7-year recurrence free\nsurvival (HR= 3.31, 95% CI 1.38 – 7.96, Cox regression p=1.74x10-3) when including\ncommon clinicopathology variables, including TNM stage (Extended Data Fig. 7c).”\nThe point about the distinction between vascular and lymphatic invasion is interesting and\ngene expression patterns are hinting at possible biological differences. However, to the\nreviewer, none of these findings prove or disprove whether lymphatic invasion is a stage in\nthe continuum, or a biologically different and unique process with clinical implications.\nWe thank the reviewer for this question about whether LI and VI represent a continuum or\ndistinct process in LUAD. We have further analyzed our data and clarified in the\nmanuscript that, although LI and VI are often combined as lymphovascular invasion (LVI) in\nclinical reports, our findings indicate that LI and VI exhibit distinct gene expression profiles\nand biological characteristics. For example, genes involved in epithelial–mesenchymal\ntransition (EMT) and tissue remodeling (Cluster 2) were strongly upregulated in VI+ tumors\nbut not in LI+ tumors, indicating a VI-specific invasion phenotype. Consistently, the VI\npredictor developed in our study did not identify LI+ cases in the absence of VI (AUC=0.51)\n– supporting that it predicts VI specifically rather than a general LVI state. Finally, the\ntranscriptional signature associated with VI was far more robust than that of LI. After\nadjusting for LI status in our model, we identified 133 genes significantly associated with VI\n(FDR < 0.01), whereas only 6 genes met this threshold for LI. Taken together, these\nobservations suggest that while LI and VI often co-occur, they are associated with distinct\ngene expression patterns and tumor biology in our cohort.\nHowever, we do acknowledge the limitations in conclusively determining whether LI and VI\nlie along a biological continuum or represent entirely separate pathways. Thus, while our\nevidence supports a biological distinction (with VI conferring unique molecular and clinical\nfeatures), we cannot exclude the possibility that LI and VI are related steps in tumor\nprogression. We have revised our conclusions in the discussion to reflect this uncertainty.\nFuture investigations, such as spatial transcriptomic analysis focused on lymphatic\ninvasion sites or longitudinal studies (e.g., tracking tumors over time) could determine\nwhether tumors progress from LI to VI or if they arise via distinct routes. Such studies will\nbe needed to establish whether LI and VI are part of a continuum of metastatic spread or\nare fundamentally distinct biological processes.\nIn the validation set, the 48 gene signature is applied to 60 samples, only 17 of which had\nVI. Would a power calculation for this analysis be appropriate?\nIn the validation set the 48 gene signature achieved an AUROC of 0.86 with a 95%\nconfidence interval of 0.76–0.96 and a p value of 6.16x10-6. This robust performance\ndemonstrates that, even with only 17 VI+ samples, the probability that the true AUROC is\nabove 0.50 is effectively 100%. Given these results, a power calculation would be of\nlimited additional value. While we acknowledge that a larger sample size would narrow the\nuncertainty around the effect size, our current analysis is adequately powered to support\nthe conclusion that the 48-gene signature is predictive of VI status. Future studies with\nlarger cohorts can help to refine these estimates further. For additional clarity, we have\nnow added the 95% confidence interval and p values in the text of the results section when\nreporting AUROC values.\nI cannot find the data for the sentence “It was associated with histologic grade and 7-year\nrecurrence free survival (HR=1.98, 95% CI 1.09–3.60, Cox regression p=0.024) (Extended\nData Fig. 6b and Fig. 4f).”\nWe apologize for the confusion. We have now separated the figure references for these\nstatements for clarity. The text now reads: “It was associated with histologic grade\n(Extended Data Fig. 7d) and 7-year recurrence free survival (HR=1.98, 95% CI 1.09–3.60,\nCox regression p=0.024) (Fig. 4f).”\nIn discussion, some of the pathophysiological explanations are not convincing or need\nmore explanation. For instance, why would a group of “cell cycle” or “tumor hypoxia”\ngenes specifically contribute to vascular invasion. Why does cluster 2 (EMT/Angiogenesis)\ncorrelate with the presence of CAFs? What is the value of necrosis? How would a gene\nsignature of vascular invasion have prognostic value in VI- tumors?\nWe thank the reviewer for these questions. We have now substantially revised the\ndiscussion section to clarify the biological interpretation our findings, drawing on evidence\nfrom our data and the literature.\nOn stylistic issues, the results section is filled with technical details, and this\noverrepresentation of details detracts from the main message of the paper and its\nproposed value as a biomarker of invasiveness.\nWe thank the reviewer for providing stylistic feedback on the results section. We have now\nmade substantial changes to the results section, which includes reporting the results of\nmany new analyses that were done to address reviewer comments. During this process,\nwe were careful to avoid including unnecessary technical details and have moved relevant\ndetails about existing and new results to the methods section.\nAlso, some of the premises of the study are not solid. For instance, authors mention that in\nthe Hattori et al study (2024) “…post-hoc analyses are beginning to show subsets of\npatients for which lobectomy might decrease the rate of recurrence.” The mentioned study\nshowed that lobar resection was associated with higher locoregional recurrence among\nthe full cohort, but survival was age and sex-dependent.\nWe appreciate the reviewer pointing out our lack of clarity here. We cited the Hattori\n(2024) study to highlight that post-hoc analysis are just beginning to be published on the\nmore recent randomized control trials but recognize this may be confusing as stated. We\nhave simplified the language of paragraph 2 of the introduction to make this more clear\nand further improved clarity in response to Reviewer #2’s feedback. We have modified the\nfollowing text in Paragraph 2 of the introduction:\nOriginal text:\nPatients with stage I VI+ LUAD may benefit from adjuvant therapy9–12, but it is difficult to\nassess in resected tumor specimens. Elastic stains can be used to improve visualization of\ninvaded blood vessels over hematoxylin and eosin (H&E), but comprehensive tumor\nhistopathology review for small (<1mm) VI foci is difficult and prone to false negatives.\nSimilarly, biopsy specimens do not provide enough tissue material to identify VI+ LUAD\nprior to surgery, preventing the use of VI status to guide surgical or neoadjuvant treatment\napproaches. Recent evidence from large clinical trials including JCOG0802/WJOG4607L\nand CALGB140503 suggests that lobectomy provides no clinical benefit over sublobar\nresection for patients with stage IA disease13,14, although post-hoc analyses are beginning\nto show subsets of patients for which lobectomy might decrease the rate of recurrence15.\nMoreover, retrospective analysis shows that VI+ patients receiving sublobar resection have\nincreased rates of recurrence, highlighting an emerging need to identify if these patients\nwill benefit from precision surgical approaches16.\nRevised text:\nPatients with stage I VI+ LUAD may benefit from adjuvant therapy9–12, but it is difficult to\nassess in resected tumor specimens. Elastic stains can be used to improve visualization of\ninvaded blood vessels over hematoxylin and eosin (H&E), but comprehensive tumor\nhistopathology review for small (<1mm) VI foci is difficult and prone to false negatives. A\nbiomarker of VI would improve existing routine histopathologic review of high-risk cases at\nresection. Moreover, biopsy specimens do not provide enough tissue material for\npathologists to identify VI+ LUAD prior to surgery, preventing the use of VI status to guide\nsurgical or neoadjuvant treatment approaches. Retrospective analysis shows that VI+\npatients receiving sublobar resection have increased rates of recurrence13. While recent\nevidence from large clinical trials including JCOG0802/WJOG4607L and CALGB140503\nsuggests that lobectomy provides no clinical benefit over sublobar resection for patients\nwith stage IA disease14,15, post-hoc analyses of pathologic factors have yet to be\npublished16. Therefore, there is an emerging need to identify if stage I patients will benefit\nfrom precision surgery. In parallel, biomarker approaches that are informed by the tumor\nmicroenvironment could also improve patient selection for neoadjuvant or periadjuvant\ntherapy in early-stage LUAD17,18.\nAuthors mention that gene signature for vascular invasion has not been studied. A number\nof studies in the literature\nincluding https://doi.org/10.1016/j.ajpath.2010.10.040, https://doi.org/10.1371/journal.po\nne.0216847, and https://doi.org/10.1038/s41416-019-0486-6, had visited this entity in\ncancer.\nWe apologize for the lack of clarity. We have re-written the introduction text to more clearly\nstate that “To our knowledge, VI-associated gene expression has not been studied in lung\ncancer.” To further highlight the novelty and robustness of our gene expression signature\nof VI, we have compared our signature with two of the three papers the reviewer referenced\nthat examined RNA differences (the third investigated miRNA differences). The first paper\nwas from Mannelqvist et al. 2011 (PMID: 21281818), who compared VI+ and VI- primary\nendometrial cancers by DNA microarray and identified a 35-gene signature VI signature (26\nupregulated, 9 downregulated). The second paper came from Kurozumi et al. 2019 (PMID:\n31114020), who derived a 99-gene signature (42 upregulated, 57 downregulated) between\nLVI+ and LVI- breast cancers by RNA-seq. We found limited overlap between the\nupregulated genes in these signatures and the upregulated genes in our signature. But\nwhen we performed gene-set enrichment analysis (GSEA) on each signature’s upregulated\ngenes against a ranked list of genes ordered by association with VI (18,285 genes\nmeasured), we found that both the 26-gene Mannelqvist 2011 signature upregulated with\nVI+ endometrial cancer and the 42-gene Kurozumi 2019 signature upregulated with LVI+\nbreast cancer were significantly enriched. This suggests that some of the molecular\naspects of VI may be shared among different epithelial cancers and that this is an area that\nshould be further explored. We have added the following text to the results section of the\nmanuscript titled ‘Gene expression changes in stage I LUAD with VI’: “Finally, genes\npreviously found to be increased in VI+ endometrial cancer30 or LVI+ breast cancer31 were\neach significantly enriched among genes most increased in VI+ stage I LUAD tumors using\nGSEA (Extended Data Fig. 3g-h). This hints at common molecular changes associated with\nVI across epithelial tumors that should be explored in future work.”\nReviewer #2 (Remarks to the Author): clinical expertise in lung cancer development and\nprogression\nThis manuscript reports a lung adenocarcinoma transcriptional signature that is\nassociated with the pathological detection of vascular invasion in resected stage 1 tumors.\nThe association was confirmed in an independent cohort of 60 patients with AUROC of\n0.86. The signature was associated with recurrence free survival and was detectable in\nregions without invasion in tumors that were annotated as harboring vascular invasion.\nThe study design is sound and the observations are robust, although not conceptually\nnovel.\n1. The determination of vascular invasion was by a single thoracic pathologist. Were\nthese cases previously read by other clinical pathologists and what was the agreement\nof determination of vascular invasion in cases read by more than one pathologist? The\nauthors should address inter-reader variability concerns with determination of this\npathological finding.\nWe thank the reviewer for raising the important point of inter-reader robustness for VI\namong pathologists. We previously published on inter-reader reproducibility for VI (PMID\n38204657 and reported substantial agreement (κ = 0.71) among 5 pathologists for the\ndiagnosis of vascular invasion. The cases that were used for assessing inter-reader\nreproducibility for VI included histopathology samples that were profiled as part of our\ndiscovery cohort. We have now added the following statement to the methods to reference\nthis prior inter-reader variability assessment: “An inter-reader reproducibility assessment\nwas previously performed on these pathology cases and found substantial agreement (κ =\n0.71) among 5 pathologists for the diagnosis of VI (PMID 38204657).”\n2. The introduction emphasizes the potential clinical significance of detection of vascular\ninvasion and suggests that the finding may support use of adjuvant therapy or lobectomy.\nThis is confusing and should be clarified.\nWe thank the reviewer for suggesting clarification of the clinical significance of VI\ndetection. We have re-written paragraph #2 of the introduction for clarity to distinguish\nbetween two use cases for a VI biomarker: 1) predicting VI using resection tissue to guide\nadjuvant therapy and 2) predicting VI pre-surgically using biopsies to guide precision\nsurgery or neoadjuvant therapy. To support the use of the VI biomarker in presurgical\nmanagement decisions, and as detailed in our response to Reviewer #1 (see above), we\nhave now validated the gene signature by RNA-sequencing of pre-surgical biopsies and\nmatched resected tumors obtained from an independent cohort of stage I LUAD patients.\nWe have modified the following text in paragraph 2 of the introduction of the manuscript:\nPatients with stage I VI+ LUAD may benefit from adjuvant therapy9–12, but it is difficult to\nassess in resected tumor specimens. Elastic stains can be used to improve visualization of\ninvaded blood vessels over hematoxylin and eosin (H&E), but comprehensive tumor\nhistopathology review for small (<1mm) VI foci is difficult and prone to false negatives. A\nbiomarker of VI would improve existing routine histopathologic review of high-risk cases at\nresection. Moreover, biopsy specimens do not provide enough tissue material for\npathologists to identify VI+ LUAD prior to surgery, preventing the use of VI status to guide\nsurgical or neoadjuvant treatment approaches. Retrospective analysis shows that VI+\npatients receiving sublobar resection have increased rates of recurrence13. While recent\nevidence from large clinical trials including JCOG0802/WJOG4607L and CALGB140503\nsuggests that lobectomy provides no clinical benefit over sublobar resection for patients\nwith stage IA disease14,15, post-hoc analyses of pathologic factors specifically associated\nwith recurrence in patients who underwent sublobar resection have yet to be published16.\nTherefore, there is an emerging need to identify if stage I patients will benefit from\nprecision surgery. In parallel, biomarker approaches could also improve patient selection\nfor neoadjuvant or adjuvant therapy in early-stage LUAD17,18.\nFurther, it is possible that the signature impact on recurrence free survival may be\nindependent of performance of lobectomy or administration of adjuvant therapy.\nAddressing the potential clinical utility aspects of the proposed classifier and how it can be\nassessed will be helpful.\nWe thank the reviewer for this insightful comment. Yes, in the TRACERx LUAD dataset we\nfind that when we add surgery type (lobectomy or sublobar resection) and administration\nof adjuvant therapy as covariates in our Cox regression model, the VI predictor scores\nremain significantly associated with decreased 5-year RFS. We have now included this\nresult in Extended Data Fig. 7h and reported this in the results section. Our hypothesis is\nthat any difference in RFS between VI+ and VI- patients in the group receiving sublobar\nresection would be attenuated in the group receiving lobar resection. These groups are not\nlarge enough in TRACERx to conduct this analysis, but we are planning a well-powered\nclinical study to address this question in the future.\nIt is important to note at the time of dataset collection adjuvant therapy was not routinely\ngiven to stage I LUAD patients, and only 4/82 patients in this analysis received adjuvant\ntherapy. Furthermore, adjuvant therapy was not administered to any patients in our\ndiscovery or validation cohorts. We have now clarified this by changing the following text in\nour methods section from “A discovery cohort consisting of 192 resected tumors from 192\npatients with 8th edition TNM stage 0/1 LUAD and not treated with neoadjuvant therapy\nwere included in this study” to “A discovery cohort consisting of 192 resected tumors from\n192 patients with 8th edition TNM stage 0/1 LUAD and not treated with neoadjuvant or\nadjuvant therapy were included in this study.” We also changed the text, “An independent\nvalidation cohort consisting of 61 resected tumors from 60 patients with 8th edition TNM\nstage I/II LUAD and not treated with neoadjuvant therapy were included in this study” to\n“An independent validation cohort consisting of 61 resected tumors from 60 patients with\n8th edition TNM stage I/II LUAD and not treated with neoadjuvant or adjuvant therapy were\nincluded in this study.” We have added the following text to the results section of the\nmanuscript: “Additionally, VI predictor scores remained significantly associated with\ndecreased 5-year RFS when including surgery type (lobectomy or sublobar resection) and\nadjuvant therapy as covariates (HR = 1.66, 95% CI 1.02-2.71, multivariate Cox regression\np=0.041; Extended Data Fig. 7h).”\n3. Figure 3 shows imputed cell type associations with the vascular invasion signature\nusing spatial transcriptomic and bulk sequencing data. It will be important to validate\nkey gene expression signatures in specific cell types using immunohistochemistry.\nWe thank the reviewer for raising the important point about validation via IHC. We have\nperformed the requested analysis as described above in response to Reviewer #1’s similar\ncomment.\nReviewer #3 (Remarks to the Author): expertise in lung ST and scRNA-seq\nIn their manuscript “Identification of a gene expression signature of vascular invasion and\nrecurrence in stage I lung adenocarcinoma via bulk and spatial transcriptomics”, Steiner et\nal. developed a transcriptomic signature of microscopic vascular invasion (VI) in stage I\nlung adenocarcinoma based on a cohort of 163 tumors. The VI-associated expression\npatterns clustered into four subclasses linked to unique biological processes. The authors\nutilized spatial transcriptomics to describe the spatial organization of the four gene\nclusters. The authors additionally developed a predictor of VI in an independent cohort and\nevaluated its performance in multiple regions of the same tumor to show the predictor is\nlargely unaffected by intra-tumor heterogeneity. The novel transcriptional signature for\nvascular invasion in early-stage lung adenocarcinoma presented in this work may enable\nmore precise diagnosis and treatment of patients.\nMy specific comments which may help to improve the manuscript are (as line numbers are\nnot available, I am referencing the figure panels to identify the relevant section of the\nmanuscript):\nExt. Data Fig. 1a-f: It would be helpful to add outcomes that include both VI and LI\nclassification given that LI is used as a covariate for some of the differential expression\nanalysis.\nWe thank the reviewer for this helpful feedback. We have now added an additional panel to\nExtended Data Fig. 1 that shows patients stratified by lymphovascular invasion (VI+ and/or\nLI+) in the same Kaplan-Meier plot and updated the corresponding figure reference\n(Extended Data Fig. 1g) in the results section.\nExt. Data Fig. 1j: In the legend, “only” is cut-off.\nWe thank the reviewer for pointing this out. We have now fixed the plot legend for the\noriginal Extended Data Fig. 1j (now Fig. 1k).\nExt. Data Fig. 1j-k: Is the site of recurrence related to time to recurrence?\nThis is an interesting question. In our discovery cohort there is no clear association\nbetween the site of recurrence and the time to recurrence. When we plot the cumulative\nincidence functions, we observe substantial overlap in the 95% confidence intervals for\nloco-regional only and distant recurrences, indicating that any differences in timing are\nminimal. We have now included this plot as an additional panel in Extended Data Fig. 1 and\nupdated the figure legend and reference to the figure (Extended Data Fig. 1m) in the results\nsection accordingly.\nFig. 1b: It may be helpful to have VI/NST/LMP as a single column annotation as they are\nmutually exclusive and based on the same grading system.\nWe thank the reviewer for this suggestion. We have updated the VI/NST/LMP column\nannotation in the Fig. 1b and Extended Data Fig. 3a-b heatmaps to better reflect their\nmutual exclusivity.\nFig. 1b: The differential expression results and cluster assignments should be provided as\nsupplementary tables or in a data repository.\nThank you for the helpful suggestion. We have now included this information in the source\ndata file.\nFig. 1b: There are two VI samples that are Ving with LMP samples. Are there any biological\nor technical factors that may explain this grouping?\nThis is an astute observation. Yes, we have provided additional details to the results\nsection of the manuscript about these two samples: “We observed that two VI+ tumors\nclustered with the LMP classified tumors. These both represented rare edge cases, with\none being the only lepidic predominant VI+ sample and the other being the only VI+ sample\nwith >50% micropapillary growth pattern. Larger cohorts containing more of these types of\nsamples will be required to determine if the current classifier might be less accurate in\ndetecting VI in these conditions.”\nFig. 1b: The NST samples are clustering with both VI and LST samples. Can you comment\non how to best interpret the two VI classifications?\nWe apologize for the confusion. The “No Special Type” (NST) histopathology grade refers\nto tumors that do not contain the required characteristics to be considered low malignant\npotential (LMP) and do not contain vascular invasion. NST is thus a ‘catch all’ for\n‘intermediate’ aggressiveness. When we derive our 48-gene VI predictor, we train it via\nsupervised learning to discriminate VI tumors from both NST and LMP tumors.\nFig. 1b: Cluster 1 and, to a lesser degree, cluster 3 appear to be separating the VI+ samples\ninto two subsets. Is there any merit to treating those as independent groups for\ndownstream analysis?\nWe thank the reviewer for this suggestion. We had performed additional analysis to\nsupport the observation that there may be two subsets of VI+ tumors present in the\ndiscovery cohort. However, when we divided the VI+ samples into the two subsets\nreferenced by the reviewer, we did not find a significant difference in RFS between them\n(log-rank p value = 0.16). The primary differences between these two VI+ samples were in\nhistopathologic pattern. The VI+ tumors with higher co-expression of cluster 1 genes had\nhigher mean % solid (39 vs. 8) and mean % cribriform (18 vs. 7). These findings further\nhighlight the importance of building our VI predictor to be independent of the LUAD high-\ngrade patterns, which we demonstrated in our stRNA-seq data.\nExt. Data Fig. 3a: It may be helpful to add a brief explanation of the biological significance\nof adding LI as a covariate.\nWe thank the reviewer for this suggestion. Adding LI as a covariate helps isolate gene\nexpression changes specific to VI by accounting for any overlap between VI and LI. We\nhave now added the following text to the results section of the manuscript (that also\naddresses the reviewer’s point below). Previously the text stated, “When we added LI as a\ncovariate to our model, we identified 133 genes associated with VI (FDR < 0.01), which\nincluded 130 genes from the original 474-gene VI signature (Extended Data Fig. 3a).”\nIt now reads: “We added LI as a covariate to our model to determine whether the gene\nexpression differences seen with VI are due to the covariance between LI and VI. We\nidentified 133 genes associated with VI (VI+ vs. LMP, FDR < 0.01), which included 130\ngenes from the original 474-gene VI signature (Extended Data Fig. 3a). At the same\nsignificance threshold (FDR < 0.01), we recovered only 6 genes associated with LI (LI+ vs.\nLI-) using the same model, with 2/6 belonging to the VI signature (Extended Data Fig. 3b),\nsuggesting the biological signal for LI is much weaker.”\nExt. Data Fig. 3b: In the main text, it is not clear which specific comparison was performed\nto get these genes.\nWe apologize for the lack of clarity. We hope the changes outlined in the preceding\nresponse also clarify this point.\nExt. Data Fig. 4d: Are all genes in all samples shown? Is the correlation similar for just the\nVI gene signature? Is the correlation consistent for all samples?\nYes, Extended Data Fig. 4d shows the expression of 17,544 genes measured by both the\nstRNA-seq and bulk RNAseq assays across all samples (n=15). The correlation is\nconsistent for all samples, with the range in spearman correlation coefficients for\nindividual samples between 0.80 and 0.87. We have added the sample-level correlation in\nthe source data file. When we consider only the correlation for the VI signature genes\npresent in both assays, we observe a similarly high correlation (0.78) across all samples.\nWe have added this correlation plot as Extended Data Fig. 4e. We have modified the\nfollowing sentence from the results section. Previously, it read: “Strong concordance was\nobserved between pseudo-bulked stRNA-seq gene expression counts and bulk RNA-Seq\ncounts of matched tumor tissues (spearman R=0.85) (Extended Daga Fig. 4d).” It now\nstates: “Strong concordance was observed between pseudo-bulked stRNA-seq gene\nexpression counts and bulk RNA-Seq counts of matched tumor tissues across all genes\n(spearman R=0.85) (Extended Data Fig. 4d) and across VI signature genes (spearman\nR=0.78) (Extended Data Fig. 4e).”\nFig. 3a: It would be helpful to have the cell types order be consistent with Fig. 3b.\nWe thank the reviewer for this suggestion. We have updated the cell type order in Fig.3a\naccordingly.\nFig. 3b: Are all the cell types listed in the legend present? For this visualization, it may be\nbeneficial to collapse all the cell types with negligible proportions.\nWe thank the reviewer for this suggestion. While all the cell types listed in the legend were\npresent, the reviewer is correct to point out that some had negligible proportions. We have\nnow revised Fig. 3b to only plot cell types that had proportions >1%.\nFig. 3b: It would be useful to additionally show the cell type frequencies across the\npathology annotations.\nWe thank the reviewer for this suggestion. We agree that this would be an informative plot.\nWe have now replaced the original Extended Data Fig. 5a panel and updated the figure\nlegend in favor of this new plot showing cell type proportions across LUAD pathology\nannotations. Previously the results section read: “Visualization of cell types within the\nstRNA-seq data showed estimated cell types localizing with expected pathology structures\n(e.g., B cells with lymphoid aggregates) (Extended Data Fig. 5a).” It now reads: “LUAD\npathology annotations showed distinct cell type composition (e.g., B cells with lymphoid\naggregates) (Extended Data Fig. 5a).”\nFig. 3b: Are there other cell types in addition to peribronchial fibroblasts with significantly\naltered abundances between the two groups?\nWe highlighted peribronchial fibroblasts and macrophages as cell types with different\nabundances between VI- and VI+ samples, as these showed statistically significant\ndifferences in the bulk RNAseq data after correction for multiple hypothesis testing. We\nhave highlighted that these are the only cell types that differ significantly between the VI-\nand VI+ samples in the Results.\nFig. 3c: It would be interesting to see the stRNA-seq proportions in addition to the bulk\nRNA-seq ones.\nAlthough we reported trends for higher proportion of peribronchial fibroblasts, B cells, and\nplasma cells with VI+ tumors in our stRNA-seq data shown in Fig 3B, these were not\nstatistically significant due to low N. Thus, we focused only on changes in proportions\nbetween VI+ and VI- tumors that we could validate by deconvolution in our bulk RNA-seq\ndiscovery cohort, after accounting for multiple hypothesis testing.\nFig. 3d: If you check the enrichment of the VI cluster signatures across all the cell types\n(analogous to Ext. Data Fig. 5b), does that yield different associations?\nWe thank the reviewer for this suggestion. We have now examined the enrichment of the VI\ncluster signatures across all the cell types in Fig 3d and included this as Extended Data Fig.\n5d for comparison. Importantly, the tumor cells present in the single-cell atlas are not the\nsame as those in our stRNA-seq data and may reflect different etiologies. We note\nsimilarities from this analysis with Fig. 3d; for instance, VI gene cluster 4 shows enrichment\nin club cells and shows some spatial correlation with club cells in our LUAD stRNAseq\ndata. VI gene cluster 2 shows higher enrichment in peribronchial fibroblasts than alveolar\nor adventitial fibroblasts, as expected. However, this analysis also highlights the limitation\nin not considering spatial location. For instance, spatial correlation in our stRNA-seq data\nis better able to detect the relationship between VI gene cluster 1 or 3 and tumor cells\nwithin the tissue. We have added the following text to the results section of the\nmanuscript: “These associations were not as apparent when examining enrichment scores\nin the scRNA-seq atlas (Extended Data Fig. 5d), highlighting the added value of spatial\ndeconvolution.”\nFig. 4a: The gene selection process includes dividing the genes into four clusters. Can you\nclarify the reasoning for the step?\nWe agree this is a useful point to clarify. We selected k=4 clusters in our original signature\n(Fig 1B.) based on several factors. First, we evaluated the biological enrichment terms\nreturned by EnrichR as reported in Fig 1C. When selecting k>4 clusters, we did not find\ndistinct and significant biological enrichment terms between all clusters, suggesting lack\nof biologically distinct gene co-expression with k>4. For a more quantitative metric, we\nexplored other values for k (k = 2-10) using consensus clustering and found that adding\nclusters greatly improved cluster stability up to but not beyond k=4. We have now included\nthe consensus clustering results as an additional panel in Extended Data Fig. 2. We also\nadded the following text to the results: “Based on our analysis of k=2-10 from consensus\nclustering, the delta area under the cumulative density function indicated substantial\nimprovement in clustering stability up to but not beyond k=4 (Extended Data Fig. 2a).”\nWe also added the following text to the methods in the “Derivation of VI signature and\nbiological annotation” section: “K=4 gene expression clusters were identified using the\nWard2 method of hierarchical clustering of genes and samples based on Euclidean\ndistance72. When we explored other values for k (k = 2-10) using consensus clustering, we\nfound that clustering stability showed substantial improvement (large change in delta area\nunder the CDF) up to but not beyond k=4. Additionally, k>4 did not return distinct\nsignificant biological enrichment terms between all clusters, suggesting a lack of\nbiologically distinct gene co-expression with k>4.” Finally, when we perform the gene\nselection process as part of our cross-validation pipeline in Fig. 4a, we preserve this\ndivision into four clusters in the outer train loop to ensure that we are consistently training\na VI predictor that contains all aspects of VI+ tumor biology.\nFig. 4b: Does increasing or decreasing the gene set size improve the performance of the\npredictor?\nWe initially picked a gene set size representing 10% of our 474 gene VI signature for ease of\nmodel evaluation, as opposed to trying to fine-tune the most optimal biomarker. However,\nwithin our cross-validation pipeline in our discovery cohort we varied our gene set size\nfrom 2.5% of the total 474 gene VI signature up to 40% of the signature and observed\nconsistent performance using our GLM model. We have now included an additional panel\nin Extended Data Fig. 7 (previously Extended Data Fig. 6) that demonstrates this lack of\ndependence of gene set size on VI predictor performance. We also added the following text\nto the results: “Predictors utilizing a binomial logit generalized linear model (GLM) with\nridge regression performed the best on the internal cross-validation test sets within our\ndiscovery cohort and we observed no dependence of gene set size on model performance\n(Extended Data Fig. 7a and 7b). We also added the following text to the methods in the “VI\npredictor derivation” section: “We observed no dependence of gene set size on model\nperformance up to 40% of the signature.”\nFig. 5a: The figure legend states “VI+ tumors that did not contain VI foci in the\ncapture area were considered as proximal VI+”. Would the tumors without VI foci be\nconsidered distal?\nWe thank the reviewer for catching this error. The Fig. 5a legend should read: “Examples of\nbinning VI+ tumors into distal VI+ (spots > 1mm outside invaded foci boundary) and\nproximal VI+ (spots < 1mm from and including the invaded foci). Any VI+ tumors that did\nnot contain VI foci in the capture area were considered as distal VI+”. We have updated the\ntext in figure legend.\nReviewer #3 (Remarks on code availability):\nThe URL provided in the Reviewer Assessment form\n(https://github.com/dst53/VI_manuscript) is different from the URL provided in the Code\nand Software Submission Checklist and the Reporting Summary\n(https://github.com/dst53/Steiner_NatComm). Neither one is accessible.\nWe thank the reviewer for bringing this to our attention. The correct URL is\nhttps://github.com/dst53/Steiner_NatComm and the code should now be accessible.\nRESPONSE TO REVIEW ERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nIn response to my comments authors have performed a study on 24 selected FFPE\nsamples and IHC and ISH for col1a1 on several slides from one sample. The following\nare my comments on these and other modifications:\n1-clinical utility:\nThe addition of the FFPE biopsies to the end of the paper increases the clinical\napplicability of the paper. However, unfortunately it falls short of the main purpose of my\ncritique which was to show the actual clinical utility of this approach. Testing the results\nof the bulk seq on 12 VI+ and 12VI- preselected samples does not mean that this\napproach will work on FFPE biopsies. This groups of samples do not have the variation\nseen in a regular sample pool and the AUROC derived from this set is of questionable\napplicability for clinical practice. Also, the placement of this section as a proof of\nconcept defies the main purpose of this critique; the fact that the main utility of this\nfinding, if validated, would be iy diagnostic utility for FFPE samples.\nR1C1. We thank the reviewer for these comments. We agree with the reviewer about the\nneed to demonstrate clinical utility prior to clinical implementation, while acknowledging\nthat true clinical utility – showing that the use of the test changes patient management and\nimproves outcomes – requires a prospective, randomized controlled trial with multi-year\nendpoints that is beyond the scope of this revision. With this in mind, we now introduce\nadditional findings related to the diagnostic performance of our classifier in our pre-\nsurgical biopsy cohort that further defines the potential clinical application of our\nbiomarker. Below is a summary of these findings:\nWe first specified a classification threshold for the biomarker score based on the Youden\npoint in the discovery cohort and used this threshold to classify the presurgical biopsy\nsamples. We now report the VI predictor has a sensitivity of 0.83 (exact binomial 95%\nconfidence interval (CI) 0.52 - 0.98) for predicting VI+ tumors at resection and a specificity\nof 0.92 (95% CI 0.62 to 1.0). As the reviewer points out, our biopsy cohort has a balanced\ndesign with equal numbers of VI+ and VI- tumors. We intentionally over-sampled VI+\ntumors relative to their real-world prevalence (~25-30%) to ensure equivalent precision for\nestimating both sensitivity and specificity in this pilot cohort.\nNext, in keeping with standard practice for estimating diagnostic accuracy\n(PMID: 35403239), we used independent and uniform Beta priors for sensitivity and\nspecificity, updated them to posterior distributions based on our biopsy cohort's true-\npositive/false-negative and true-negative/false-positive classification counts, simulated\n10,000 specificity/sensitivity pairs and then applied Bayes' theorem to estimate positive\npredictive value (PPV) and negative predictive value (NPV) over a realistic range of VI\nprevalence (15-35%).\nDespite our small, preselected biopsy cohort, the signature's prevalence-adjusted negative\npredictive value (NPV) estimate is high and precise (median 93%, credible interval (CI) 84-\n98% at 25% prevalence, dropping to median 91%, CI 81-98% at 30% prevalence),\nsuggesting that a negative VI predictor classification may be used to recommend sublobar\nresection given the low probability that vascular invasion will be detected on final\npathology of the resected specimen. Notably, because our balanced biopsy cohort design\nunder-samples true negatives relative to observed real-world prevalence, these CIs are\nconservative and would tighten in an unselected cohort.\nConversely, the prevalence-adjusted positive predictive value estimates from this biopsy\ncohort are less precise (median 73%, CI 47%-95% at 30% prevalence), reflecting both the\nsmall number of VI+ biopsy cases and the moderate prevalence of VI+ stage I LUAD in the\nreal-world population. As a result, it is unclear if a positive VI call will be subicient to\nrecommend lobectomy. Larger prospective studies and incorporation of orthogonal assays\n(e.g. radiomics) will be needed to narrow PPV confidence bounds and demonstrate the\nsignature's clinical utility for guiding more extensive surgery.\nWith the current n=24 FFPE biopsy specimens we have demonstrated the classifier's\npotential to be measured in routine biopsy tissue and its high NPV. Larger, prospective\ncohorts will be needed to determine its PPV and determine true clinical utility. We have\nnow included the projected NPV and PPV result as a new panel in Figure 5 and updated the\nmethods, results and discussion with this new analysis, including the limitations described\nabove.\n-The addition to the TNM staging:\nThe statistical analyses showing that VI remains an independent predictor of RFS with\nTNM stage added is satisfactory. (fig1i extended) However,\n-I cannot find the clinical characteristics of this 183 tumors that are “a subset” of\npreviously published data (Yambayev et al). I can not find 1823 samples in that paper.\n-Also: what percentage of the data is missing for each variable?\n-About the age and pack year in these samples. Since these are continuous variables\nwhy are there no CIs?\n-Extrinsic factors: Even though author’s response about focusing on VI as an intrinsic\nfactor (vs extrinsic factors such as clear margins and suboptimal sampling) is valid,\nthe extrinsic features will still have importance in determining the predictive power of a\nbiomarker in real world and should at least be mentioned in the discussion.\n2-Clinical feasibility:\nThe authors in their response mention that in a real-world scenario the remainder of an\nFFPE after pathologist review will probably be used for bulk-seq, stressing the\nimportance of the FFPE sample validation. Yet perhaps the weakest point of this\nmanuscript at present is the analysis of the FFPE samples (24 preselected VI+ and VI-\nsamples). Increasing the number of the FFPE samples and testing this approach on a\ntrue representative sample pool (rather than as an add-on proof of concept piece) will\nshow the basic utility of this study.\nThe comment about the use of RNA-seq in current clinical practice: the sentence “whole\ntranscriptome RNA-seq test….fine-needle aspiration FFPE material used … is now\nstandard of care for preoperative diagnosis…” should be moderated. The use of\nmolecular probes for indeterminate thyroid nodules refers to a small group of samples\nassessed by FNA, and mRNA sequencing is only one of the three approaches used for\nthese samples: mutational analysis, mRNA sequencing, and miRNA expression\nanalysis. Authors are presenting Bulk-seq as a method for general classification of all\nearly stage NSCLC samples, which is a very different premise.\nAlso, on the disparities between RNA-seq and qRT-PCR please refer to:\n-Int. J. Mol. Sci. 2021, 22(17), 9349; https://doi.org/10.3390/ijms22179349\nDespite the extensive discussion on the advantages of sequencing by some authors\nand commercial entities, I think we all know that the access to RNAseq is much more\nlimited and technical caveats are much more prevalent, compared to hybridization\nprobes or qPCR. These methods are currently used for biomarker analysis in many\nclinical settings and for various cancers, including breast and colon cancer. Therefore,\nvalidation by qPCR or RNA probes for some genes, which is routinely done for RNAseq\nstudies is reasonable and adds to the possible clinical utility.\n3-pthophysiological significance:\nThis is one of the main critiques regarding the appropriateness of this paper for Nature\ncommunications. Despite the extensive work and the multitude of statistical analyses,\nthe conclusions in the manuscript are based on correlative data and deficient in\nmechanistic insights, which is a prerequisite for publication in journals like Nature\nComm. The involvement of CAFs in vascular invasion is known. However, the main\nmechanisms mentioned for this contribution are increased angiogenesis, expression of\nECM remodeling enzymes, immune response, and metabolic programming of tumor\ncells.\n(https://doi.org/10.3389/fimmu.2025.1582532)\nAuthors have picked Col1a1 for validation. This gene is expressed as part of a fibrotic\ntissue response to injury/invasion. Mechanistic questions remain. For instance: Does\ncol1a1 contribute to the process of invasion? If so, how? Or is it just a consequence?\nDoes it identify “activated” CAFs? Is there any mechanism downstream of col1a1 that\ncan be shown to enhance any aspect of vascular invasion? In a mechanistic paper\nthese questions need to be answered by experiments, not just correlation and citing\nliterature.\n-Critique on association with a specific cell type: staining for one of the genes.\nAuthors have performed IHC and ISH staining for col1a1, on multiple sections from one\nsample. Question arises, why could this not be done on more than one sample, as is\nregularly done for the histopathological validation of a gene involved in a pathological\nprocess? Staining even one or two sections from 10-20 samples for col1A1 and\npresenting the histopathological images that show the expression of this gene in\nmyofibroblasts within the conventional “vascular invasion” sites would be appropriate for\nthis response. The multiple statistical analyses show that within sample 4 the\nexpression of col1A1 coincides with the expression of the other genes in cluster 2 in\nthat sample. This kind of analysis adds to the rigor but does not provide an adequate\nvalidation and link between this new finding and the currently used histopathological\nstudies that are the basis of the statistical analysis in extended figure 1i.\n-Multivariable analysis on the derivation study: Extended Figure 7c is adequate and a critical piece of data. The same\nquestions about the CI for age and pack year apply and need to be answered. If the latter questions are adequately\nanswered, I suggest moving this data to the main figures.\n-On the distinction between LI and VI : the changes are adequate. I would add that LI was also an independent predictor in extended fig.1i and would show the gene signature associated with it as a variable in extended figure 7c.\n-The explanations for tumor hypoxia and cell cycle genes are adequate.\n-There are still a multitude of technical details in result sections that take away from the main message in each paragraph.\nThe literature references are adequate.\nAdditionally:\nThe manuscript still has several stylistic and language issues and would benefit from cleaning up.\nSome examples:\n-What is the meaning of the sentence “ tumors were graded…. Or VI, is associated with outcome than the World…”\n-The text still has redundant and hard to understand details that do not clarify the significance: For instance:\n“ based on our analysis from consensus clustering the delta area under the cumulative….” This sentence and sentences alike should be modified/simplified or deleted.\n-Where can we find the list of the genes for the four clusters described in the second results paragraph?\n-In discussion paragraph 7:\n“High VI predictor scores were associated not only with… VI- tumors were shedding cells into the circulation…”\nThis sentence is confusing and needs clarification.\nReviewer #2 (Remarks to the Author):\nThe authors have responded to the main comments that are listed below.\n1. The determination of vascular invasion was by a single thoracic pathologist. Were\nthese cases previously read by other clinical pathologists and what was the agreement\nof determination of vascular invasion in cases read by more than one pathologist? The\nauthors should address inter-reader variability concerns with determination of this\npathological finding.\nThe authors’ response is complete.\n2. The introduction emphasizes the potential clinical significance of detection of vascular\ninvasion and suggests that the finding may support use of adjuvant therapy or\nlobectomy. This is confusing and should be clarified. Further, it is possible that the\nsignature impact on recurrence free survival may be independent of performance of\nlobectomy or administration of adjuvant therapy. Addressing the potential clinical utility\naspects of the proposed classifier and how it can be assessed will be helpful.\nThe discussion addresses these points. To further address potential clinical utility, the\nauthors added analysis of a cohort of presurgical FFPE lung biopsy specimens acquired\nfrom 24 subjects in Figure 5f. The AUROC is shown. To better understand the\nperformance of this proposed biomarker, it will be important to show the sensitivity,\nspecificity, PPV and NPV of the predictor. In addition, it is important to carefully describe the specimens from which the predictor was evaluated. Were these consecutively acquired specimens or were they selected? Were the specimens derived from cell block or forceps biopsy? What was the evidence for vascular invasion on the biopsy specimens used for the analysis, if from forceps biopsy?\n3. Figure 3 shows imputed cell type associations with the vascular invasion signature using spatial transcriptomic and bulk sequencing data. It will be important to validate key gene expression signatures in specific cell types using immunohistochemistry.\nIn extended figure 6, The authors show analysis of COL1A1IHC and ISH in one sample. This does not address validation of key V1 biomarkers cell specific expression in lung adenocarcinoma specimens with and without vascular invasion. These data will be important to demonstrate the generalizability of the reported observations.\nReviewer #3 (Remarks to the Author):\nThe authors addressed my comments and introduced satisfactory corrections.\nRESPONSE TO REVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nThe authors' responses are satisfactory.\nReviewer #1 rebuttal\nR1C1.\n1. To address comments regarding diagnostic accuracy, the authors present results\nfrom 24 selected biopsy specimens acquired from a collaborator. The results are\nsupportive but are insufficient to support the statement that the classifier should be\nmeasured in routine biopsy tissue. Validation of diagnostic accuracy in a multi-center\nclinical trial is required to support this claim.\nWe thank the reviewer for this comment. We agree with the reviewer’s concern\nregarding diagnostic accuracy and clinical validation of the classifier in routine biopsy\ntissue. As detailed in our responses below, we have revised the Abstract, Introduction,\nResults, and Discussion to clarify that the biopsy cohort represents proof-of-concept\nfeasibility rather than evidence of routine clinical applicability, and to explicitly state that\nmulti-center clinical trials will be required to establish diagnostic accuracy and clinical\nutility.\n2. Abstract. The results do not support the claim that “VI-associated gene expression\nchanges …..can be used to predict the presence of VI from routinely collected pre-\nsurgical tumor biopsy specimens” The results do support the need for further studies to\nvalidate the diagnostic utility and to demonstrate the clinical utility of this assay.\nWe thank the reviewer for this comment. We have now revised the previously second to\nlast sentence of the abstract (now the last sentence) to emphasize that VI-associated\ngene expression changes are detectable in pre-surgical biopsies and support further\nvalidation. The sentence now reads: “These findings indicate that VI-associated\ntranscriptional changes extend across the tumor and are detectable in limited biopsy\nmaterial, supporting further validation for preoperative risk stratification in stage I\nLUAD.”\n3. Introduction The final sentence needs to be revised for clarity and accuracy.\nSpecifically, the term “robustness” should be clarified and the phrase “validated the VI\npredictor in real-world pre-surgical tumor biopsies” needs to be rewritten for accuracy-\nthe data presented are not acquired from a “real-world” setting.\nWe thank the reviewer for these comments. We have now revised the final sentence of\nthe introduction to clarify what is meant by robustness (i.e., stability across intra-tumor\nheterogeneity and sampling variability) and to remove the term “real-world,” which could\nbe misinterpreted. The revised text now accurately reflects the controlled, pilot nature of\nthe biopsy cohort. The sentence now reads: “We derive a novel molecular signature\nassociated with VI+ LUAD using bulk RNA-seq, describe its association with\nhistopathology in situ using stRNA-seq, develop and validate a machine learning based\npredictor of VI, evaluate its stability across intra-tumor heterogeneity (ITH) and\ndemonstrate its proof-of-concept feasibility in a pilot cohort of pre-surgical tumor\nbiopsies.”\n4. Results. Similarly, the phrase “real-world population” should be deleted from the final\nparagraph.\nWe thank the reviewer for this comment. We have removed the phrase “real-world\npopulation” from the final paragraph of the Results section. The sentence now reads\n“Conversely, the prevalence-adjusted PPV estimates from this biopsy cohort are less\nprecise (median 73%, CI 47%-95% at 30% prevalence), reflecting both the small\nnumber of VI+ biopsy cases and the moderate prevalence of VI+ stage I LUAD.” We\nhave also removed the only other use of “real world” in the final paragraph of the\ndiscussion. The sentence “Our biopsy results support real world feasibility of predicting\nVI from presurgical tissue.” has been rewritten to “Our biopsy results support the\ntechnical feasibility of predicting VI from presurgical tissue.”\n5. Discussion. The first paragraph on page 19 needs to be revised to place the results in\nappropriate context in light of the limitations noted above.\nWe appreciate the reviewer’s suggestion to further contextualize the biopsy findings\nconsidering their limitations. While several revisions had already been made to this\nsection in prior rounds, we have now further revised the first paragraph on page 19 of\nthe discussion to more explicitly emphasize the pilot nature of the biopsy cohort, remove\nlanguage that could be interpreted as implying clinical decision-making, and clearly\ndistinguish technical feasibility from diagnostic or clinical utility.\nR1C3. Thank you for clarifying the confusion about samples listed in Extended Data\nFigure 1i and for specifying that data for 16 of the 192 tumors is missing. The authors\nshould examine and report the potential impact of the missing data on the multivariate\nanalysis. Simply excluding the missing data may introduce bias.\nWe thank the reviewer for raising this point. As noted in the figure legend for Extended\nData Fig. 1i, missing covariate data affected a small subset of tumors (16/192) primarily\ndue to missing race or smoking pack-year information, as well as mucinous tumors that\nwere excluded from WHO pathology grading. To assess the potential impact of\nexcluding these cases, we performed a sensitivity analysis by refitting the multivariate\nCox model in the full cohort (n=192 tumors) using a reduced covariate set that avoided\nexclusion due to missing data. The association between vascular invasion and\nrecurrence remained strong and statistically significant (HR = 7.07, 95% CI 2.97–16.8,\nmultivariate Cox regression p=9.82x10-6), indicating that exclusion of cases with missing\ncovariates does not affect the conclusions. We have now clarified this in the figure\nlegend."}]

Decision letter