Overall, 31 tumor tissue samples were taken from 14 male and 11 female donors who were diagnosed between 2003 and 2018 with SIA (23 pilocytic astrocytomas, 4 glioblastomas, 2 anaplastic astrocytomas, and 2 astrocytomas with uncategorized histological subtype; Figure 1 ). The mean age was 32.73 years (range 18–69 years) and 30.00 years (18–56 years), respectively. The biosamples were formalin-fixed paraffin-embedded (FFPE) histologically characterized tumor samples with at least 70% cancer cells. Clinical annotations of tumor tissue specimens investigated in this study and their patient origin are summarized in Supplementary Table S1 . Kaplan–Meier plots for progression-free survival and overall survival are shown in Supplementary Figure S1 .

A total of 10,509.3 million reads were obtained for 31 independent libraries of SIA tissues (median 55.3 million reads per sample). Most reads reached Phred-like quality scores (Q-scores) at the Q30 level, indicating that the probability of an incorrect base call is 0.001%. The average coverage of sequencing depth reached approximately 53.45× of the human transcriptome. After alignment, 98.57% to 99.12% uniquely aligned reads were mapped to the reference human genome.

To further characterize SIA transcriptomic data, we compared using principal component analysis (PCA) distributions of RNA-seq profiles among the SIA samples (n = 31) and publicly available datasets of normal spinal cord (n = 159) and brain (n = 9) tissues from GTEx and ANTE databases, respectively. PCA was performed to investigate cross-dataset compatibility in order to select proper reference group(s) for SIA comparison. The profiles from ANTE database were chosen because they were obtained using the same reagents, equipment, and protocols as for the current experimental SIA sampling (24). The GTEx reference group of samples was selected because this is, to our knowledge, currently the biggest publicly available collection of healthy spinal cord RNA-seq profiles [30]. Performing two-step expression analysis allowed us to select the DEGs between SIA and normal neural tissue explored using the same RNA-seq platform (brain samples from ANTE) and using a different platform but for the same tissue type (GTEx spinal cord samples). Unfortunately, normal spinal cord samples were not available in the ANTE database. We hope that this approach allowed to establish differential gene expression profiles without the influence of platform-specific batch effect.

PCA was performed first in the space of log₁₀ transformed quantile normalized gene counts. We observed tissue-specific sample clustering corresponding to the biological nature of the datasets under analysis, where SIA samples formed a separate cluster ( Figure 1 A ) . In addition, we performed PCA for brain and spinal cord from GTEx, SIA, and ANTE normal brain data normalized using quantile normalization, DESeq2, and harmonized/batch corrected using XPN (26), CuBlock (27), and Shambhala (28). It appeared that GTEx brain profiles clustered with ANTE normal brain even in case of DESeq2 normalization, indicating that further batch correction was not necessary ( Supplementary Figure S2 ).

Then, we performed PCA based on pathway activation levels (PALs) of 1,611 molecular pathways (29) calculated using the same transcriptomic data for each sample under study ( Figure 1B ). In case of pathway upregulation or downregulation, PALs can take positive or negative values, respectively, thus quantitatively reflecting the extent of a pathway activation or inhibition relatively to the control group of samples. Zero PAL values suggest unaffected activity of a molecular pathway. Thus, PAL values can be used as the quantitative functional characteristic of the interactome under analysis (14). We calculated PALs according to Borisov et al. (12), with the 168 healthy CNS tissue samples taken as the controls.

On the PAL-based PCA plot, we observed similar clustering as for the gene expression-based PCA ( Figure 1B ), thus strongly suggesting that SIA samples should be independently compared to each of the above healthy CNS tissue datasets.

Subsequently, we performed paired differential gene expression analysis between SIA samples relatively to each dataset of healthy CNS tissues ( Supplementary Table S2 ). Overall, 1,949 genes, 1,766 (90.61%) upregulated and 183 (9.39%) downregulated, were found to be statistically significantly differentially expressed [|log₂ FC|>5, false discovery rate (FDR)-adjusted p-value <0.05] between SIA and GTEx spinal cord samples ( Figure 2A ). In turn, 382 DEGs, 102 (26.70%) upregulated and 280 (73.30%) downregulated, were found for the comparison between SIA and ANTE healthy brain samples ( Figure 2B ).

These DEG sets were then intersected with respect to log₂ FC sign ( Figures 3A, B ). In order to test whether an observed number of common differential genes can support random or non-random intersection hypothesis, we performed perturbation test for randomness according to Sorokin et al. (30) with 1,000 random gene sets. The percentile of the observed case precedent in the distribution of random intersections was considered as a measure of statistical significance.

In total, 32 genes, 23 (71.88%) upregulated and 9 (28.12%) downregulated, were commonly differentially expressed in SIA samples according to both comparisons, which supported the hypothesis that the intersections between the DEGs were non-random (p < 0.01).

To evaluate potential functional similarities of the above 32 SIA-specific differential genes and the underlying molecular and cellular processes, we then performed Gene Ontology (GO) terms enrichment analysis ( Figures 4A, B ). We identified significantly enriched 563 functional GO terms, where 340 terms (60.39%) were for upregulated and 223 (39.61%) were for downregulated DEGs. For statistical estimates, we used Benjamini–Hochberg method for FDR correction (31) and p-value threshold 0.05 (32).

Interestingly, most of the enriched terms for upregulated DEGs were related to the regulation of an innate and adaptive immune response, thus supporting a concept that immune microenvironment may play a crucial role in the development of spinal astrocytomas (33). In contrast, for downregulated DEGs, the most strongly enriched terms were linked with complex neuronal processes, such as cognition, learning ability, and regulation of neurotransmitter secretion and transport. The latter supports specific functional impairments occurring in astrocytomas in comparison with healthy CNS tissues (34).

We then performed differential PAL analysis for SIA samples relative to healthy CNS tissues ( Figure 5 ). In total, 1,611 molecular pathways including 10 and more gene products were interrogated from Oncobox pathway databank (29).

When comparing SIA and GTEx healthy spinal cord samples, we identified 468 differentially regulated pathways, 272 (58.12%) of them were activated and 196 (41.88%) were suppressed ( Figure 5 ). In turn, for the comparison between SIA and ANTE normal brain samples, 464 differential pathways were identified; among them, 254 (54.74%) were activated and 210 (45.26%) were suppressed ( Figure 5 ).

The intersections between these pathway sets returned 397 common differential molecular pathways [222 (55.92%) activated and 175 (44.08%) inhibited; Supplementary Tables S3 , S4 ]. This supported non-random intersection between the differential pathways in both comparisons, p < 0.001 ( Figure 5 ).

In both comparisons, top activated pathways deal with intracellular signal transduction and with the immune response, whereas top downregulated pathways are responsible for translational regulation and neurotransmitter activities ( Figure 6 ). Except for the new feature of translational regulation, this trend was in line with the results obtained previously for the GO terms enrichment in the SIA differential genes ( Figure 4 ).

A chemotherapeutic treatment of SIA remains a challenging and poorly investigated field, and we performed computational simulation whether anticancer targeted drugs (ATDs) that are currently in use for other CNS tumors could be repurposed as second-line treatment options for SIA. To this end, we utilized Oncobox method for predicting efficiencies of ATDs based on gene expression and molecular pathway activation data (35). This returns for every drug a tumor sample-specific value of balanced drug efficiency score [drug score (DS)]. DS reflects an expected responsiveness of a tumor to a specific drug, where higher values mean higher expected efficacy of an ATD. Furthermore, drugs with positive DS are predicted to be potentially beneficial, and drugs with negative DS—potentially harmful (35). This method was shown to be clinically beneficial in a prospective clinical investigation on high-grade human solid tumors [16] and was effective for individual selection of experimental/off-label chemo- and targeted therapeutic settings, e.g., Buzdin et al. (36) and Moisseev et al. (37). In glioblastoma, Oncobox method could effectively predict tumor response on temozolomide, a DNA-alkylating agent whose activity is antagonized by MGMT gene products (38).

By using Oncobox algorithm, we identified 85 ATDs with positive DS in the SIA–GTEx spinal cord comparison, and 70 ATDs with positive DS in the SIA–ANTE healthy brain comparison ( Figures 7 , 8 ). In these lists, there were 66 common drugs, thus evidencing non-random intersection between the two comparison results ( Figure 8 and Supplementary Tables S5 , S6 ).

We then assessed available clinical trial reports for the top 20 DS-ranked drugs among these common 66 ATDs with the biggest DS values for different CNS tumors ( Table 1 and Figure 7 ). Interestingly, the top predicted drugs mostly represented the classes of tyrosine kinase inhibitors (i.e., regorafenib, lenvatinib, nintedanib, sorafenib, dovitinib, sunitinib, tivozanib, pazopanib, imatinib, foretinib, dasatinib, erdafitinib) and thalidomide analogs (thalidomide and pomalidomide). Many of these drugs were previously investigated for CNS tumors and related cancers like neuroblastoma ( Table 1 ).

On the other hand, drugs with the predicted negative drug scores that were, therefore, algorithmically not recommended belonged mainly to cyclin-dependent kinase inhibitors and androgenic and anabolic steroid (AAS) classes.

Chromosomal rearrangements resulting in fusion genes and abnormal transcripts in some cases may become clinically actionable targets of specific cancer therapeutics (71). Fusion transcripts combine exons of 2 or more genes and may serve as the oncogenic drivers in many cancers including CNS tumors (72, 73).

We used RNA-seq profiles for SIA patients to detect fusion transcripts presenting in spinal astrocytoma and focused on the fusions where at least on partner gene was a serine/threonine or tyrosine kinase. This allowed to select potentially druggable fusion genes. We found, in total, 16 different fusion transcripts identified by STAR-Fusion software (74) ( Figure 9 and Supplementary Figure S3 ). One of them, KIAA1549-BRAF, was found in two SIA patients and preserved BRAF kinase domain ( Figures 10A, B ). Interestingly, this fusion transcript was previously reported to confer a clinically less aggressive phenotype in pediatric low-grade astrocytoma (75). It was also found less abundant in the adult compared to pediatric patients with pilocytic astrocytoma (76). Other fusions, to our knowledge, were not reported previously and were not found in ChimerDB fusion database.

This study was performed in agreement with the ethical principles of Declaration of Helsinki. Retrospective biosamples were obtained from patients diagnosed with SIAs who had undergone surgery at the spinal department of Burdenko Neurosurgical Center, Moscow. From all the patients involved or from their legal representatives, informed written consents to participate in this study were collected. The study design and consent collection procedure were approved by the local ethical committee of the Burdenko Neurosurgical Center. For all patients enrolled and for their biosamples, the consent was obtained for publication of age, sex, histological tumor type, diagnosis, and molecular data including RNA-seq profiles but not including whole-genome and/or whole-exome sequencing data.

Biosamples were FFPE tumor tissue blocks that were evaluated and confirmed by a pathologist who estimated a proportion of tumor cells and determined the histological type of a tumor. In this study, only FFPE blocks with at least 70% tumor cells were analyzed. In total, 31 samples for 25 SIAs meeting the above criteria were obtained for further molecular screenings ( Supplementary Table S5 ).

RNA libraries were generated and sequenced according to Suntsova et al. (24). RNA was extracted using RecoverAll™ Total Nucleic Acid Isolation Kit (Invitrogen). RNA concentrations were measured with Qubit RNA Assay Kit, and Agilent 2100 bioanalyzer was used to measure RNA Integrity Number (RIN). Depletion of ribosomal RNA was performed using RNA Hyper Kit (Roche), and then library concentrations and fragment length distributions were measured with Qubit (Life Technologies) and Agilent Tapestation (Agilent), respectively. The RNA-seq was performed using Illumina NextSeq 550 engine for 50-bp single-end reads and approximately 30 million raw reads per sample using standard protocol. Single-end sequencing was used because SIA samples were FFPE tissue blocks that typically have a strong degree of RNA degradation. Primary sequencing data quality control was performed with Illumina SAV, and demultiplexing was made according to Suntsova et al. (24) with Illumina Bcl2fastq2 v 2.17 software.

SIA profiles were processed according to Suntsova et al. (24). STAR aligner (91) was used to process FASTQ files from RNA-seq in “GeneCounts” mode for Ensembl human transcriptome annotation GRCh38.89. The gene names for 36,596 annotated genes were converted to HGNC (HUGO Gene Nomenclature Committee) gene symbols from Ensembl IDs according to Complete HGNC dataset, version of August 17, 2021 (https://www.genenames.org). Further quality control metrics for RNA-seq data were obtained with NCBI MAGIC software (92–94). RNA-seq profiles were preprocessed by quantile normalization method (95), and then differential expression analysis was performed using DESeq2 (96), and visualized with R package EnhancedVolcano (97). Genes that were considered significantly differentially expressed had to pass a threshold of FDR-adjusted p-values <0.05 (98). GO enrichment analysis was conducted using clusterProfile (v.4.2.1) and org.Hs.eg.db (v.3.8.2) R packages. Fusion transcripts were detected with STAR-Fusion tool (74), and PCA and visualization were done for log₁₀ transformed counts of all genes using pca2d R (v.3.6.0) and prcomp software. Code for data analysis is available at: https://github.com/raevskymichail/SIA_analysis.

We used healthy tissue transcriptomic profiles obtained for normal human spinal cord biosamples from GTEx project portal (25) and for normal brain samples from the ANTE database (24). In total, 9 ANTE and 159 GTEx normal CNS samples were analyzed. Raw count quantification in GENCODEv26 annotation was obtained from GTEx portal.

Algorithmically annotated molecular pathway graphs were taken from our previously published database (29). PALs were calculated with the Oncobox bioinformatic platform. It allows quantitative assessment of PALs using RNA-seq data and functionally annotated collection of molecular pathways (12, 29). We used a set of 1,611 pathways with 10 or more gene products included because of previously reported poor theoretically estimated data aggregation effect for smaller pathways (15).

This method of calculating PAL showed a strong potential to suppress batch effects (15, 16, 21) and to minimize the artifacts introduced by the methods of experimental transcriptome analysis (13, 99). An absolute value of PAL reflects the strength of the pathway up/downregulation, while a positive or negative sign indicates its activation or suppression, accordingly (12). To calculate PAL, each sample RNA-seq profile was normalized on mean geometric levels of gene expression in the relevant control dataset.

Drug score [Balanced Efficiency Score (BES)] for cancer-targeted drugs was calculated according to Tkachev et al. (35), whose method is based on the analysis of targeted molecular pathway activation and relative expression levels of drug target genes.

To test whether an observed number of overlapping differential genes or pathways between the two intersecting datasets is significant, for every comparison, we performed 1,000 random intersections according to Sorokin et al. (30). In every case, two random samples from the corresponding gene sets under comparison were taken. Then, these random samples were intersected for 1,000 iterations, and numbers of randomly obtained common genes were registered. Then, p-value of intersection significance was calculated as an expected fraction of random intersects that is equal to or higher than the experimentally observed number of overlapping genes.