
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39256467
72335
10.1038/s41598-024-72335-w
Article
Full-length target sequences of GeoMx digital spatial profiling probes reveal that gene-promiscuity predicts probe sensitivity to EDTA tissue decalcification
Oszwald André 1
Zisser Lucia 2
Schachner Helga 1
http://orcid.org/0000-0002-7795-2726
Kaltenecker Christopher 1
Wasinger Gabriel 1
Rohrbeck Johannes 1
Kozakowsky Nicolas 1
Tiefenbacher Andreas 1
Rees Andrew J. 1
http://orcid.org/0000-0002-2428-543X
Kain Renate renate.kain@meduniwien.ac.at

1
1 https://ror.org/05n3x4p02 grid.22937.3d 0000 0000 9259 8492 Department of Pathology, Medical University of Vienna, Währinger Gürtel 18-20, 1090 Vienna, Austria
2 https://ror.org/05n3x4p02 grid.22937.3d 0000 0000 9259 8492 Division of Nuclear Medicine, Department of Biomedical Imaging and Image-Guided Therapy, Medical University of Vienna, Währinger Gürtel 18-20, 1090 Vienna, Austria
10 9 2024
10 9 2024
2024
14 211568 6 2023
5 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
GeoMx Digital Spatial Profiling (Nanostring) is a commercial spatial transcriptomics method to selectively analyze regions of interest within intact tissue sections. We show that decalcification with ethylene-diamine-tetra-acetic (EDTA) variably attenuates probe counts, while probes that are more resistant to this effect consequently appear overexpressed after quantile normalization. By determining the undisclosed full-length target sequences of probes used in the human whole transcriptome panel, hereby updating target transcripts and genes, we find that the gene-promiscuity of probes is an important factor that determines sensitivity to EDTA incubation.

Subject terms

Gene expression analysis
Computational biology and bioinformatics
General departmental funds of the department of Pathologyissue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Spatial profiling provides unique opportunities to study biology in health and disease by enabling transcriptomic analysis of specific cell types or tissue compartments within their topohistological context1–3. GeoMX Digital Spatial Profiling (DSP) (Nanostring, Seattle) quantifies transcripts indirectly using hybridizing DNA probes conjugated to photo-cleavable DNA barcodes that are eluted by spatially controlled ultraviolet illumination4. Formalin fixed, paraffin-embedded (FFPE) tissues are compatible with this method and can be used to capture transcriptomic profiles of highly flexible microscopic regions of interest. Archived clinical FFPE specimens may stem from different times or sites of one individual, or from well-characterized patient cohorts, and are typically collected and stored under standardized conditions, minimizing pre-analytical variability. However, the comparison of samples from mineralized and non-mineralized compartments, such as primary tumors and their bone marrow metastases, is an important exception because mineralized tissues must first be decalcified, typically with ethylene-diamene-tetra-acetic acid (EDTA), before they can be sectioned5.

Currently, the essential data on whether the decalcification protocol affects DSP analysis is missing, and so we assessed the influence of a standard EDTA protocol on the DSP assay.

EDTA decalcification reduces probe counts while retaining the overall profile

We compared results of 24 sets of matched tissue samples from different organs of 20 patients, each consisting of four samples taken from adjacent tissue of clinical specimens. Within each set, we left one sample untreated and incubated the others with EDTA for one, three, or seven days, respectively (Fig. 1A, raw data and annotations in Table S1). Total target probe counts (i.e., the sum of all quantified target molecules within a region of interest) and their interquartile ranges decreased significantly after three (p < 0.05) and seven days (p < 0.0005) EDTA treatment compared to untreated tissues (Fig. 1B). By contrast, EDTA had no effect on the limit of quantification (LOQ, a threshold derived from non-target control probe counts, used to discriminate expressed genes from background noise), thus reducing the signal-to-noise ratio. Consequently, the number of probe species (of 18,676 biological target probes present in the WTA panel) detected above the LOQ decreased progressively during 7 days of treatment (from approximately 8000 to 5000 probes, p < 0.005) (Fig. 1B). Nevertheless, the overall expression profiles of samples were maintained (more clearly so after quantile normalization), as shown by continued tissue- and sample-specific clustering in tSNE plots (Fig. 1C), despite a progressively worsening correlation between the probe counts of matched EDTA-treated and untreated samples (Fig. 1D). This implies that EDTA does not have a uniform effect, but that some probes are more sensitive to it than others.Fig. 1 EDTA tissue treatment reduces signal-to-noise ratio and affects differential gene expression and pathway analysis. (A) Schematic outline of experiment. 24 different tissues were sampled from 20 different patients. Each tissue was divided into four similar parts and either left untreated as control, or incubated with EDTA solution for one, three, or seven days. Tissues were then assembled on a tissue microarray for the DSP run. In two tissues, two measurements were made (each from distinct tissue compartments). In all other tissues, one measurement per tissue was made. Two regions of interest had to be excluded due to sampling errors. Created with BioRender.com. (B) Top: sum and interquartile ranges of all probes per region of interest. Bottom: Limit of quantification (LOQ), calculated per region of interest as geometric mean of all negative probes * (geometric standard deviation of all negative probes ^ 2) according to default settings suggested by the manufacturer; number of probes per region of interest with counts higher than the regions’ respective LOQ. (C) t-SNE plots (perplexity 15) of probe counts, using non-normalized (left) or quantile-normalized data (right). Colors indicate tissue type, numbers indicate clinical specimens from which matched adjacent tissue samples were obtained (D) Correlation coefficients (Pearson, Fisher Z-transformed) between probe counts of each region of interest (untreated vs one, three, or seven days of EDTA incubation). (E) Differential gene expression analysis of probe counts between untreated control tissues and those incubated with EDTA for one (top), three (middle), or seven days (bottom), using non-normalized (left) or globally quantile-normalized data (right). (F) Gene Set Enrichment Analysis (GSEA) using the differential expression analysis data from quantile-normalized counts between untreated samples and those treated for seven days with EDTA was performed using https://www.webgestalt.org. Significance levels: ns not significant; *, p < 0.05; **, p < 0.005; ***, p < 0.0005.

To quantify this and to identify probes with particular EDTA sensitivity, we scaled the probe counts of each set to the respective untreated sample, in order to compensate for intrinsic differences and performed differential expression analysis. When using the non-normalized data, we observed a large number of significant differences with a clear bias towards genes with lower expression (Fig. 1E), indicating probes that are sensitive to EDTA treatment and represent technical variability. In this curated collection of samples, the observed differences in gene expression are a result of the decrease in probe counts after EDTA incubation and not true differences in gene expression. Consequently, in other situations, in order to deduce subtle biological differences between similar sets of decalcified and untreated samples (e.g., non-decalcified primary tumors vs. paired decalcified bone metastases), it would be necessary to remove this technical effect prior to performing differential gene expression analysis. We therefore sought to determine if quantile normalization, a common method to mitigate technical variation between samples6, would suitably remove the technical effect of EDTA, and to identify probes that appeared to be resistant or artificially altered by normalization. We found that quantile normalization (both global, i.e., across all samples, and within biological groups, i.e., tissue types) strongly reduced the number of significant differences induced by 7 days of EDTA treatment (5051 non-normalized vs. 571 and 625, respectively) and yielded a more balanced ratio of genes with significantly higher or lower counts (Fig. 1E). Thus, quantile normalization is effective in compensating the overall decrease in counts caused by EDTA incubation, but causes a subset of probes to appear overexpressed. To demonstrate the consequences of EDTA and quantile normalization on downstream analyses, we performed gene set enrichment analysis on the differential expression results between untreated samples and those treated for seven days, we identified several pathways that appeared significantly altered after normalization, including those related to antigen processing and DNA-damage (Fig. 1F). Thus, we concluded that the specific differences induced by EDTA are important because they may skew comparisons between decalcified and non-decalcified tissues, even when applying quantile normalization globally or within groups. To better understand this effect, we aimed to characterize the affected probes, which required knowledge of the full-length probe sequences.

In-silico characterization of probe target sequences

DSP probes are 35 to 50 nucleotides long, and are identified in the publicly available configuration file7 by a 35-nucleotide extract, the identifiers of mRNA transcripts they target, and the genomic coordinates that define the borders of the genomic sequences that give rise to the targeted transcript sequences. 98% of the 126,346 transcripts specified in the configuration file could be identified in a past RefSeq8 database version (GrCh38.p13 release 2020/05/22). However, because transcript databases are repeatedly updated with transcripts being added or removed, potential targets of DSP probes need to be repeatedly re-evaluated; this is critical because a more current version of RefSeq (GrCh38.p14, accessed 2023/04/08) allows identification of only 42% of the transcripts in the configuration file. Accordingly, we developed a strategy for inferring the undisclosed target sequences by combining the specifications provided by Nanostring with publicly available data (Fig. 2A).Fig. 2 Analysis of full length target sequences of DSP probes. (A) The DSP configuration file is publicly available (see reference 6) and lists targeted transcripts, the genomic coordinates that delimit the genomic sequence that gives rise to the targeted sequence of the target mRNA transcript, and a 35-nucleotide long extract of the probe target sequence. However, the corresponding transcript coordinates first need to be identified using RefSeq and the GenomicFeatures R package. Transcript coordinates allow the retrieval of potential target sequences, which are then filtered by those matching the truncated sequence provided by Nanostring. Sequences are then aligned to known human transcripts using BlastN. Created with BioRender.com. (B) Overview of results of the probe identification process. “Probe-gene association” is an identified pair of probe and target gene. One probe may have multiple probe-gene associations. Each probe-gene association may involve multiple transcripts. Percentages are rounded to two decimal points and are relative to the number of the “parent” element (C) Probes ranked by adjusted p value of differential gene expression analysis between untreated tissues and those incubated with EDTA for one, three, or seven days (vertical axis). On the horizontal axis are shown the number of gene targets per probe (length of bar) and the direction of the log2FC of the differential expression analysis (color and direction of bar). Blue bars indicate negative log2FC, whereas red bars indicate a positive log2FC. Note that in non-normalized data, probes with many target genes (long bars) are located at the bottom end of the graph after 7 days (i.e., show only small change after EDTA treatment). After quantile normalization, these probes appear to show the highest change and are located at the top of the graph. The dotted line indicates the portion of probes with significant differential expression.

For each probe, the configuration file specifies a variable number of target genes, transcripts and genomic coordinates, but not the correspondence between these data. For this reason, we first used RefSeq to determine matching pairs of target transcripts and their genomic coordinates. This allowed us to accurately convert the genomic coordinates to coordinates within targeted transcripts, which in turn enabled us to retrieve the precise transcript sequences that the probes are specified to hybridize with (summarized in Fig. 2B).

The Nanostring human whole transcriptome atlas panel7 comprises 18,815 probes (18,676 target probes and 139 negative control probes). We determined the unequivocal full-length target sequences of 18,300 of the target probes (~ 98%), whilst 376 remained ambiguous: 213 with missing genomic coordinates (including two erroneous coordinates); 91 matching multiple potential target sequences; 56 probes with up to three nucleotides differing from the given extract; 12 probes whose given extract lay within a specified transcript but outside of the given genomic coordinates; and 4 probes where the coordinates indicated a sequence longer than the limit of 50 nucleotides specified by Nanostring.

We further determined that the majority of unequivocally identified probes do not span introns (15,105, 83%), and that a minority of probes (26, 0.14%) span introns only in a subset of their target transcripts. The negative control probes (designed not to align against any human sequence) have no specified genomic coordinates or target transcripts, and therefore could not be deduced. We provide the entire dataset containing all probes, including unequivocal and ambiguous sequences in (Supplementary Table 2), but in the following summary only report results from unequivocally identified probes (e.g., excluding negative control probes), unless explicitly stated.

We next aligned the target sequences with known human transcripts using BlastN9. We identified 182,161 transcripts with complete alignment identity and 11,568 transcripts with incomplete (81% to < 100%) alignment identity. Unexpectedly, even among “official” target transcripts previously specified by Nanostring, we found a minority (391) with incomplete (89%-98%) alignment. Grouping the transcripts by probe and gene, we identified a total of 21,867 probe-gene associations. Similarly, even among those previously specified by Nanostring, we found a minority (321) with incomplete alignment (93%-98%). This suggests that complete sequence identity was not an absolute requirement for Nanostring to designate target transcripts or genes. We therefore also refrained from further filtering our BlastN hits by alignment identity.

We determined that the majority of the probe-gene associations (16,050) involved more than one transcript per gene, whereas a minority involved only a single one (5817; overall, with a median of 4 and mean of 8.8 transcript isoforms per probe). To identify novel transcript isoforms, we compared the base accession numbers of transcripts (without version suffixes) between those specified by Nanostring and those identified via BlastN, yielding 90,634 novel transcripts, and 30,641 that are either not targeted in their current version, or were entirely removed (and not updated) by RefSeq. Transcripts that were merely updated (i.e. showed an increment in the version suffix) were not considered novel and excluded from this summary.

Grouping transcripts by probe, we found that the majority of unequivocally identified probes (16,707) targeted only one gene, whereas a minority (1559) targeted more than one gene, and 34 probes appeared not to target any transcripts currently listed in RefSeq. 1255 probes had at least one target gene that differed from the Nanostring specifications.

Summarizing all aligned transcripts, we identified 19,937 target genes (19,277 with complete alignment), including 19,061 of the 19,505 genes specified by Nanostring. When including not only unequivocally identified probes, but also given extracts and ambiguous sequences, the number of target genes was 20,347 (19,707 with complete identity), including 19,403 of those specified by Nanostring. This indicates that some probes were designed to bind targets no longer thought to be expressed as mRNA, and the corresponding gene can no longer be considered to be assessed by the panel.

The given extracts of two negative control probes showed complete alignment with transcripts listed in RefSeq, with the caveat that their full-length sequence could not be identified.

In conclusion, this approach identified the full-length sequences and corresponding targets in the human whole transcriptome atlas panel, and enabled search for properties that contributed to specific EDTA sensitivity.

Gene promiscuity of probes predicts sensitivity to EDTA calcification

Inspection of the normalized scaled probe counts showed that sensitivity of probes to EDTA incubation was associated with the number of genes targeted by the probe. However, the directionality of this association was dependent on whether the data had been quantile normalized or not (Fig. 2C). Using non-normalized data, gene-promiscuous probes tended to show the least differential expression, suggesting they were most resistant to the reduction in probe counts by EDTA incubation. In contrast, after normalization, gene-promiscuous probes tended to show the greatest differential expression, suggesting that these probe counts are inappropriately augmented by quantile normalization. Using normalized data, multiple linear regression between the B statistic of the test for differential expression, the number of targeted genes, the number of targeted transcripts, the length of target sequences and the GC content of target sequences, confirmed a robust association between the number of target genes per probe and EDTA sensitivity at three and seven days (t = 14.6 and t = 15.7, p < 2e−16). There was an independent weaker but significant association with the number of target transcripts per probe (t = −6.4 and −4.8, p < 1.4e−6). The association with target sequence length was only significant when the GC content (t = 5.5 and t = 5.6, p < 1.7e−8) was left out of the model.

Discussion

A detailed assessment of the effect of decalcification ideally requires comparison of matched decalcified and non-decalcified samples of the same specimen, in order to compensate for biological and pre-analytical variables. When using DSP (or any tissue section-based method), this ideal set of samples cannot be realized using mineralized tissues as decalcification is a prerequisite for cutting them. In routine practice, only mineralized tissues are subjected to decalcification. In order to approximate an ideal set of samples, we assembled a representative set of specimen that could be sectioned without decalcification and subjected them to our decalcification protocol for increasing intervals of duration. Hence, although we did not use mineralized tissues for our study, in our opinion the data conclusively show that EDTA treatment for one day has a small effect; while decalcification for longer durations (three to seven days) is feasible, but certain effects on global and specific probe counts should be expected. Bone marrow trephines are typically decalcified only for several hours to 2 days, but larger resections containing cortical bone may require up to several weeks 5,10.

We show that EDTA reduces probe counts and may represent an important source of technical variation which needs to be controlled before performing differential expression analysis. Probes with a high number of target genes appear to be less sensitive to EDTA; we hypothesize that this may be due to their lower specificity and therefore higher relative abundance of target molecules after EDTA treatment, thus retaining higher counts than other probes. Subsequently, quantile normalization effectively mitigates the overall reduction in counts (as determined by fewer differentially expressed genes), but leads to artificial over-expression of gene-promiscuous probes.

Regarding the artefactual overexpression after normalization, quantile normalization may not be appropriate to control technical variation from prolonged EDTA incubation, which appears to have a non-uniform effect on probes. Alternatively, global normalization methods may be inappropriate due to substantial biological variation among samples (in this case, samples from different tissue types). Smooth quantile normalization was recently developed for this situation by introducing a weighted average of global and within-group normalization11. However, our observed effect is independent of whether normalization is performed globally or within biological groups. Other methods to control for systematic technical variation include batch effect adjustment12. However, in a hypothetical experiment (e.g. comparing primary tumors and bone metastases), ComBat12 requires internal controls (i.e., EDTA-treated and untreated samples of both primary and bone metastases), which are unlikely to be available from clinical tissue archives. Our expression data (which is available via GEO accession GSE272995) may be useful as an external reference dataset for researchers aiming to identify and remove technical effects of EDTA and subsequent normalization artefacts from similar experiments.

Probes with multiple targets are specified in the DSP configuration file, but are not clearly advertised. It is likely that many of these genes cannot be distinguished on basis of short sequences due to strong homology; accordingly, this ambiguity may not constitute a design flaw but an inevitable shortcoming of whole-transcriptome analysis via short sequences.

Our results could be specific to EDTA treatment or reflect the prolonged incubation in aqueous solution at room temperature required therefore, but are representative of decalcified clinical FFPE specimens routinely processed by pathology laboratories. Consequently, gene-promiscuous probes and those with a high GC content should be carefully validated in DSP studies using calcified and non-calcified samples when using quantile normalization. Nevertheless, we demonstrate that valuable DSP data can be acquired from archived pathology specimens subjected to EDTA decalcification protocols. Importantly, we also provide a method for users to acquire the full-length DSP probe target sequences, which allows re-alignment of the deduced probe sequences to obtain updated target transcripts and genes, and greatly enhances the assays’ documentation, transparency and utility.

Materials and methods

Probe target sequence analysis

Nanostring provides publicly available configuration files for DSP assays as JSON-formatted text files7 containing assay- and probe-level information including unique probe identifier codes, targeted genes and transcripts, the genomic coordinates that flank the exonic sequences comprising the hybridization sequence, and a truncated version of the hybridization sequence.

Using the R packages GenomicFeatures13, BSGenome14, and the NCBI RefSeq annotation sets of the GrCh38.p13 genome assembly, we retrieved the full hybridization sequences as outlined in brief as follows (full code including specification of other used packages is available via supplementary material):Although the DSP configuration file specifies the GrCh38.p13 genome build, it does not specify which version RefSeq version was used to identify the transcripts that are contained in the file. We determined that the GrCh38.p13 RefSeq release of 2020/08/15 contains almost all transcripts in the DSP configuration file, but that a few transcripts are only present in other versions. In order to maximise the number of target transcripts that we could use in our analysis, we combined transcript collections of several RefSeq GrCh38.p13 releases (see code for details). We then determined transcripts that are present in both the DSP configuration file and the merged annotation sets (including all but 8 transcripts of the DSP configuration file)

The configuration file specifies transcripts and genomic coordinates in an n:n ratio, which is why we needed to determine corresponding pairs before proceeding. We searched for these pairs for each probe individually by determining if a provided genomic coordinate lies within an exon sequence of one of the provided transcripts. For each identified pair of transcript and genomic coordinates, we then converted the genomic coordinates to transcript coordinates.

We retrieved the transcript sequences delimited by the transcript coordinates using the GRCh38.p13 genome (prepared using the BSgenome package, see code) and filtered the results, retaining only the sequence that contain the provided truncated sequence. This was necessary due to the unspecified n:n ratio of transcripts and coordinates in the configuration file, which can yield multiple (off-target) results. Probes with no retrieved sequences were investigated by manual comparison of the provided truncated sequence and the retrieved sequence(s) to identify the source of the error.

Probe target sequences that were unequivocally identified were aligned to known human transcripts using BlastN9 (using RefSeq and RefSeq select collections, word size = 16).

Tissue collection

For the purpose of this study, we prepared series of non-calcified tissues that were incubated with EDTA for increasing time periods, and matched control samples from the same source tissue, which were left untreated. Specimens submitted for routine tissue analysis were collected for research purposes after approval of the ethics commission of the Medical University of Vienna (2167/2021, registered at https://ekmeduniwien.at/core/catalog/2022/) as graphically shown in Fig. 1A). In the course of routine pathology workup, 24 clinical tissue specimens (Table S1) of different organs stemming from 20 patients were fixed in 7.5% buffered formalin for 2.0 to 3.5 (mean 3.3, median 3.0) days according to the laboratory’s routine procedure. At macroscopic dissection, an excess piece of tissue was removed and divided into 4 macroscopically similar samples. Tissues were then either processed immediately in a Tissue-Tek VIP® 6 AI (Sakura Finetek) or incubated at room temperature in EDTA solution (Tritiplex, Glatt-Koller, 403212270) for 1, 3, and 7 days before processing them analogous to the untreated tissue. Samples were embedded in paraffin blocks according to standard protocols using low-melting point Paraffin (Histo-Comp, ATS-200856).

Microarray construction

Blocks were sectioned, stained with hematoxylin/eosin and scanned (Pannoramic 250 Flash II, 3DHistech) to confirm the content of equivalent histological structures among samples sourced from the same specimen. Areas for tissue microarray (TMA) cores were selected using the digital slide scans. Tissue microarrays were then assembled using an automated TMA device (TMA Grandmaster, 3DHistech) using 1 mm punches and a recipient block of 5mm depth. TMA were tempered by alternating storage at 37 °C, 4 °C, and room temperature, for 5 times (2h each), following by a single incubation at 65 °C for 10 min.

Digital spatial profiling

Slides were sectioned using a rotating microtome at 5 µm thickness. One section was used for DSP analysis as outlined below, while an adjacent section was used for validation of the TMA cores using hematoxylin & eosin stain.

DSP was performed according to protocols provided by Nanostring using the following standard parameters. In brief, sections were deparaffinized and subjected to heat-induced protein epitope retrieval (using retrieval buffer at pH 9, Thermo Fisher Scientific, 00-4956-58) for 15 min at 100 °C using a steam cooker. Retrieval for mRNA was performed using proteinase K (Thermo Fisher Scientific, 25530049) at 1µg/ml for 15 min. Slides were then incubated with the human Whole Transcriptome Atlas (WTA) probe panel (Nanostring, NA-GMX-RNA-NGS-HuWTA-4, Lot HWTA21003) at 37 °C in a humidified environment overnight. Sections were washed according to standard protocols and direct immunofluorescence was performed with antibodies against CD34 conjugated to Alexa-647 (Novus Biologicals, clone QBEnd-10, NBP2-34713AF647), Pan-cytokeratin conjugated to to Alexa532 (Novus Biologicals, clone AE1/AE3, NBP2-33200), CD45 conjugated to Alexa594 (clone 2B11 + PD7/26, Novus Biologicals, NBP2-34528), and a nuclear counterstain (Syto-13, Nanostring).

DSP device run and library preparation were performed according to the manufacturers’ protocol. The library was sequenced using an Illumina NextSeq 550 sequencer with the specifications provided by Nanostring (paired-end reads at length 27, index length 8) at a library concentration of 1.6 pM with 5% PhiX spike-in using a NextSeq 550 High Output 75-cycle kit (Illumina, 20024906).

Nuclei counts of regions of interest were obtained via the default automatic nuclear segmentation settings of the DSP device software.

Data analysis

Read counts were processed in R statistics15 (4.2.0) using R Studio (2022.02.3). After TMA core assembly, two samples were excluded from analysis due to inadverted selection of non-representative tissue. For comparisons of total probe counts, interquartile ranges, level of quantification (LOQ) and number of genes above quantification, we used two-sided pairwise Wilcoxon tests with Holm’s correction applied. For the comparison of correlation coefficients, we used a two-sided pairwise t-test, also with Holm’s correction. In the plots, the box hinges correspond to the first and third quartile of the distributions, and the middle line corresponds to the median.

Differential gene expression analysis was performed using limma16. Graphs were generated using ggplot217, volcano plots using enhancedVolcano18. Prior to differential expression analysis, we removed the counts of probes where we identified target sequences but did not find any corresponding transcript sequences with 100% alignment in BlastN, following to our independent probe hybridisation sequence analysis. We restricted the analysis of differential probe counts to those with counts greater than the LOQ in at least 10% of samples in order to focus on those most likely to be biologically informative13. The absolute probe counts were first log2-transformed and then normalized using quantile normalisation to compensate for the non-specific EDTA-induced count reduction. The normalized probe counts were then scaled to the probe counts in the matched untreated samples to create a ratio suitable for comparing relative EDTA-susceptibility of probes, regardless of baseline expression level.

Supplementary Information

Supplementary Information 1.

Supplementary Information 2.

Supplementary Information 3.

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-72335-w.

Acknowledgements

This work was funded by departmental funds of the Department of Pathology of the Medical University of Vienna.

Author contributions

André Oszwald conceived the project, performed all experiments, data analysis, and drafted the manuscript. Lucia Zisser assisted in drafting the manuscript and critical manuscript review. Helga Schachner provided essential technical support in performing experiments. Christopher Kaltenecker provided essential support in data analysis. Gabriel Wasinger, Johannes Rohrbeck, Nicolas Kozakowsky, Andreas Tiefenbacher participated in sample collection and critical manuscript review. Andrew J Rees and Renate Kain reviewed and edited the manuscript.

Data availability

Raw data and processed data are available at the NCBI GEO repository under accession GSE272995. Processed data and code for R statistics software used to identify probe target sequences are available in the supplementary materials.

Competing interests

A.O. has presented at conferences for Nanostring, Inc., where travel and accommodation was provided by Nanostring, Inc. A.O. has also presented at conferences for Illumina, Inc. and received compensation by Illumina, Inc. All other authors declare no potential conflict of interest.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Marx V Method of the year: Spatially resolved transcriptomics Nat. Methods 2021 18 9 14 10.1038/s41592-020-01033-y 33408395
Marx, V. Method of the year: Spatially resolved transcriptomics. Nat. Methods 18, 9–14 (2021).33408395 10.1038/s41592-020-01033-y
2. Zhang Q The spatial transcriptomic landscape of non-small cell lung cancer brain metastasis Nat. Commun. 2022 13 5983 10.1038/s41467-022-33365-y 36216799
Zhang, Q. et al. The spatial transcriptomic landscape of non-small cell lung cancer brain metastasis. Nat. Commun. 13, 5983 (2022).36216799 10.1038/s41467-022-33365-y
3. Carter JM Distinct spatial immune microlandscapes are independently associated with outcomes in triple-negative breast cancer Nat. Commun. 2023 14 2215 10.1038/s41467-023-37806-0 37072398
Carter, J. M. et al. Distinct spatial immune microlandscapes are independently associated with outcomes in triple-negative breast cancer. Nat. Commun. 14, 2215 (2023).37072398 10.1038/s41467-023-37806-0
4. Merritt CR Multiplex digital spatial profiling of proteins and RNA in fixed tissue Nat. Biotechnol. 2020 38 586 599 10.1038/s41587-020-0472-9 32393914
Merritt, C. R. et al. Multiplex digital spatial profiling of proteins and RNA in fixed tissue. Nat. Biotechnol. 38, 586–599 (2020).32393914 10.1038/s41587-020-0472-9
5. Sterchi, D. L. 17—Bone. In Bancroft’s Theory and Practice of Histological Techniques (eds. Suvarna, S. K., Layton, C. & Bancroft, J. D.). 18th Ed. 280–305. 10.1016/B978-0-7020-6864-5.00017-7 (Elsevier, 2019).
6. Bolstad BM Irizarry RA Åstrand M Speed TP A comparison of normalization methods for high density oligonucleotide array data based on variance and bias Bioinformatics 2003 19 185 193 10.1093/bioinformatics/19.2.185 12538238
Bolstad, B. M., Irizarry, R. A., Åstrand, M. & Speed, T. P. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19, 185–193 (2003).12538238 10.1093/bioinformatics/19.2.185
7. GeoMx DSP Configuration Files. NanoString. https://nanostring.com/products/geomx-digital-spatial-profiler/geomx-dsp-configuration-files/.
8. O’Leary NA Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation Nucleic Acids Res. 2016 44 D733 745 10.1093/nar/gkv1189 26553804
O’Leary, N. A. et al. Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 44, D733-745 (2016).26553804 10.1093/nar/gkv1189
9. Zhang, Z., Schwartz, S., Wagner, L. & Miller, W. A greedy algorithm for aligning DNA sequences. J. Comput. Biol. J. Comput. Mol. Cell Biol. 7, 203–214 (2000).
10. Ramsay, A., Pomplun, S. & Wilkins, B. Tissue Pathways for Lymph Node, Spleen and Bone Marrow Trephine Biopsy Specimens.
11. Hicks SC Smooth quantile normalization Biostatistics 2018 19 185 198 10.1093/biostatistics/kxx028 29036413
Hicks, S. C. et al. Smooth quantile normalization. Biostatistics 19, 185–198 (2018).29036413 10.1093/biostatistics/kxx028
12. Johnson WE Li C Rabinovic A Adjusting batch effects in microarray expression data using empirical Bayes methods Biostatistics 2007 8 118 127 10.1093/biostatistics/kxj037 16632515
Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8, 118–127 (2007).16632515 10.1093/biostatistics/kxj037
13. Lawrence, M. et al. Software for Computing and Annotating Genomic Ranges. PLoS Comput. Biol. 9 (2013).
14. Pagès, H. BSgenome: Software Infrastructure for Efficient Representation of Full Genomes and Their SNPs (2023).
15. R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).
16. Ritchie ME limma powers differential expression analyses for RNA-sequencing and microarray studies Nucleic Acids Res. 2015 43 e47 10.1093/nar/gkv007 25605792
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015).25605792 10.1093/nar/gkv007
17. Wickham H Ggplot2: Elegant Graphics for Data Analysis 2016 Springer
Wickham, H. Ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).
18. Blighe, K., Rana, S. & Lewis, M. EnhancedVolcano: Publication-Ready Volcano Plots with Enhanced Colouring and Labeling (2022).
