2.1 GBM explants cohort
2.1.1 GBM explants cohort description: In this study, we used the data presented by Morelli et al. (21), These data are related to 33 samples of glioblastoma multiforme (GBM) that were cultured as explants, an ex vivo pre-clinical model consisting of minimally handled fresh tumor tissues. Out of the 33 GBM samples, 18 were sequenced as fresh-frozen samples, while the other 15 as GBM explants cultured in vitro. For the sake of simplicity, we will refer to this dataset as the “GBM explants”.
All GBM explants were exposed to temozolomide (TMZ) drug. The protocol for drug treatment, as well as the culture methods of GBM explants, are described in Morelli et al. (21).
To determine which GBM explants were responsive to the benchmark drug (temozolomide, TMZ), Morelli et al. labelled the samples as responder (RES) or non-responder (NRES) based on the NAD(P)H-fluorescence lifetime imaging microscopy (FLIM) technique (21). In this context, FLIM exploits the intrinsic auto-fluorescence molecular properties of NAD(P)H, a metabolic enzymatic cofactor associated with the metabolic state of the cell/tissue.
2.1.2 RNA-seq data, alignment and annotation: We aligned the data using STAR (version 2.7.9a) against the GRCh38 version of the human genome. Read counts were generated during the alignment process (as described on the GDC guideline: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). We obtained the raw count matrix using an R script. All subsequent analyses were conducted using R version 4.2.0 within the RStudio 2022.07.02+576 suite or Python version 3.10.
2.1.3 Batch correction: Although recent works reported that fresh-frozen and cultured tissues have similar gene expression levels (23, 24), we investigated the presence of a batch effect in the RNA-seq data due to the different sample types. With this aim, we performed principal component analysis (PCA) using the PCATools R library. PCA revealed that the samples labelled by sample type clustered along the principal components 1 and 7, confirming the presence of a batch effect. To correct for this, we applied the ComBat-seq algorithm, which is specifically designed for RNA-seq data (25). In addition, to preserve the biological variability of the samples, we set the ComBat-seq argument "group" as the NAD(P)H-FLIM classification (RES/NRES). Finally, we applied PCA again to verify that the corrected data were no longer clustered by sample type.
2.1.4 Differential expression analysis: After the RNA-seq data of the 33 GBM samples were batch-corrected, we normalized them and conducted a differential expression (DE) analysis using the DESeq2 algorithm (26). The DE analysis compared the two groups of samples identified by the NAD(P)H-FLIM technique: RES (16 samples) and NRES (17 samples). To reduce the effect of genes with low counts, we corrected the estimated log-fold values using the normal shrinkage method (26). We identified differentially expressed genes (DEGs) as the ones with a fold change greater than 1.5 in absolute value (log2-fold change greater than 0.58) and p value adjusted for multiple comparisons less than 0.1. A custom R script was used to perform the analysis.
2.2 RNA-seq data of TCGA
2.2.1 TCGA RNA-seq cohort description: To validate the NAD(P)H-FLIM based classification, we conducted a comparable DE analysis on the GBM RNA-seq data sourced from The Cancer Genome Atlas (TCGA). We downloaded the RNA-seq dataset and the corresponding clinical data from the Genomic Data Commons (GDC) Data Portal. The raw counts table was obtained using the same alignment and annotation procedures as the GBM explants. From the whole set of patients, we only selected the ones treated for at least one cycle with TMZ (104 patients).
2.2.2 Differential expression analysis: Among the TMZ treated patients, we only selected the ones whose PFI value was non-censored. Then, we grouped them based on the progression-free interval (PFI): we defined two groups of patients with extreme PFI values, one group with PFI values within the 1st quartile (Short PFI, or S-PFI group, 22 patients) and the other with PFI values within the 4th quartile (Long PFI, or L-PFI group, 22 patients). Finally, the differential expression (DE) analysis was conducted on the selected RNA-seq data following the same approach and statistical thresholds used for the GBM explants dataset. Specifically, we used the DESeq2 algorithm to normalize raw counts data and then we conducted a DE analysis comparing the L-PFI with the S-PFI group (26). Differentially expressed genes (DEGs) were identified as the ones with a log2-fold change (corrected using normal shrinkage) higher than 0.58 (fold change higher than 1.5) and a p-value adjusted lower than 0.1. Data selection was performed using a custom Python script, while the DE analysis was conducted using a custom R script.
2.3 Functional-clinical differentially expressed genes
By only considering the DEGs whose expression pattern was consistent (preserved up- or down-regulation) among the two datasets, we found 19 DEGs shared by the two DE analyses. Subsequently, in order to select the most informative genes, we performed a supervised feature selection. With this aim we first standardized the RNA-seq expression values of the 19 shared DEGs of both GBM explants and TCGA. Next, we aggregated the two standardized datasets into a unique dataset (77 samples) and the Long PFI (L-PFI) label was converted to TMZ responder (RES), as well as the Short PFI (S-PFI) label was converted to TMZ non-responder (NRES). Finally, we applied a linear discriminant analysis (LDA, (27)) and selected the genes having a discriminant coefficient (i.e. the loading) higher than the mean of the coefficients, as shown in Figure 3B. This process led to the selection of 6 genes among the initial 19. The described analyses were conducted in Python using custom-made scripts.
2.4 Microarray data of TCGA Agilent
Here, we selected a different dataset: the TCGA AgilentG4502A_07_2 Microarray, downloaded from UCSC (University of California Santa Cruz) Xena (https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.GBM.sampleMap%2FAgilentG4502A_07_2.gz). We decided to use a new dataset to avoid any information leakage. Additionally, we excluded the 36 patients who were also present in the RNA-seq dataset of TCGA. To be consistent with the previous analyses, we only selected patients who were treated with TMZ, resulting in a final dataset of 247 patients.
2.5 Building a predictive model
To make the results as interpretable as possible, we implemented the Cox proportional hazards model using the lifelines library developed in Python (28, 29). An important advantage of using a Cox regression is that a hazard ratio (HR) is computed for each of the model’s covariates, which, in our case, are the 6 selected genes. Each HR gives insights on the effect of the corresponding covariate on the predicted clinical variable (progression-free in our case), i.e., whether it increases or decreases the event probability. Finally, the model can be used to divide subjects into two risk groups, i.e., high risk or low risk, by using the partial hazard (PH) assigned to each patient by the model.
2.6. Leave-One-Subject-Out (LOO) training and testing procedure
To train and test the model we adopted a Leave-One-Subject-Out (LOO) approach, rather than a classical training and testing procedure. Practically, the dataset was split into train and test set, with the latter consisting of only one sample at a time. A Cox model was then trained using all the samples except for the one left out, which was used to test the model. This procedure was repeated as many times as the total number of patients, with each patient left out of the training set in turn and treated as a new observation by the model. It is important to remark that the train set related to two different iterations differs from just one observation, i.e., the left-out-patient.
Finally, for each patient in the dataset we obtained a risk prediction, which is called partial hazard when using the Cox regression. We split the patients into high and low risk groups by using as threshold the median of the computed partial hazards. By conducting a univariate Cox regression that uses as covariate the risk group, we can compute the corresponding hazard ratio (HR), which indicates the model ability to correctly categorize patients into high and low risk groups.
2.7. Statistical analyses
Throughput this paper, we used the Kaplan-Meier plot to visualize time-to-event data. Log-rank test (https://lifelines.readthedocs.io/en/latest/lifelines.statistics.html) was applied to compare the progression-free interval (PFI) distributions of two sets of data. The Wilcoxon rank-sum test (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wilcoxon.html) was used to compare the risk scores distribution of MGMT-methylated and MGMT-unmethylated patients.