Joint Modeling of RNAseq and Radiomics Data for Glioma Molecular Characterization and Prediction

RNA sequencing (RNAseq) is a recent technology that profiles gene expression by measuring the relative frequency of the RNAseq reads. RNAseq read counts data is increasingly used in oncologic care and while radiology features (radiomics) have also been gaining utility in radiology practice such as disease diagnosis, monitoring, and treatment planning. However, contemporary literature lacks appropriate RNA-radiomics (henceforth, radiogenomics) joint modeling where RNAseq distribution is adaptive and also preserves the nature of RNAseq read counts data for glioma grading and prediction. The Negative Binomial (NB) distribution may be useful to model RNAseq read counts data that addresses potential shortcomings. In this study, we propose a novel radiogenomics-NB model for glioma grading and prediction. Our radiogenomics-NB model is developed based on differentially expressed RNAseq and selected radiomics/volumetric features which characterize tumor volume and sub-regions. The NB distribution is fitted to RNAseq counts data, and a log-linear regression model is assumed to link between the estimated NB mean and radiomics. Three radiogenomics-NB molecular mutation models (e.g., IDH mutation, 1p/19q codeletion, and ATRX mutation) are investigated. Additionally, we explore gender-specific effects on the radiogenomics-NB models. Finally, we compare the performance of the proposed three mutation prediction radiogenomics-NB models with different well-known methods in the literature: Negative Binomial Linear Discriminant Analysis (NBLDA), differentially expressed RNAseq with Random Forest (RF-genomics), radiomics and differentially expressed RNAseq with Random Forest (RF-radiogenomics), and Voom-based count transformation combined with the nearest shrinkage classifier (VoomNSC). Our analysis shows that the proposed radiogenomics-NB model significantly outperforms (ANOVA test, p < 0.05) for prediction of IDH and ATRX mutations and offers similar performance for prediction of 1p/19q codeletion, when compared to the competing models in the literature, respectively.

RNA sequencing (RNAseq) is a recent technology that profiles gene expression by measuring the relative frequency of the RNAseq reads. RNAseq read counts data is increasingly used in oncologic care and while radiology features (radiomics) have also been gaining utility in radiology practice such as disease diagnosis, monitoring, and treatment planning. However, contemporary literature lacks appropriate RNA-radiomics (henceforth, radiogenomics) joint modeling where RNAseq distribution is adaptive and also preserves the nature of RNAseq read counts data for glioma grading and prediction. The Negative Binomial (NB) distribution may be useful to model RNAseq read counts data that addresses potential shortcomings. In this study, we propose a novel radiogenomics-NB model for glioma grading and prediction. Our radiogenomics-NB model is developed based on differentially expressed RNAseq and selected radiomics/volumetric features which characterize tumor volume and sub-regions. The NB distribution is fitted to RNAseq counts data, and a log-linear regression model is assumed to link between the estimated NB mean and radiomics. Three radiogenomics-NB molecular mutation models (e.g., IDH mutation, 1p/19q codeletion, and ATRX mutation) are investigated. Additionally, we explore gender-specific effects on the radiogenomics-NB models. Finally, we compare the performance of the proposed three mutation prediction radiogenomics-NB models with different well-known methods in the literature: Negative Binomial Linear Discriminant Analysis (NBLDA), differentially expressed RNAseq with Random Forest (RF-genomics), radiomics and differentially expressed RNAseq with Random Forest (RF-radiogenomics), and Voom-based count transformation combined with the nearest shrinkage classifier (VoomNSC). Our analysis shows that the proposed radiogenomics-NB model significantly outperforms (ANOVA test, p < 0.05) for prediction of IDH and ATRX mutations and offers similar performance for prediction of 1p/19q codeletion, when compared to the competing models in the literature, respectively.
Keywords: RNA sequencing, radiomics, radiogenomics, negative binomial, molecular mutation INTRODUCTION Radiomics is increasingly being applied to radiology practice in disease diagnosis, grading, monitoring, and treatment planning (1,2). Radiomics is extracted from various radiological images of a targeted area of the disease. Fusing the important radiomics and genomics information in the proper computational machine learning (ML) model may helpto achieve a more comprehensive disease diagnosis, prognosis, and treatment planning scheme (3)(4)(5). Different studies have evaluated the association between glioma molecular subtypes and radiomics (e.g., tumor shape and size) (6)(7)(8), or between different form of genomics (e.g., RNA sequencing (RNAseq) gene expression, protein expression, copy number, molecular mutations, or DNA methylation) and glioma subtypes (9)(10)(11).
Conventional ML models do not adequately model the count-based nature of the RNA-sequence data as these models are usually designed to work with data that has a normal distribution. In order to alleviate the lack of appropriate ML models, researchers propose to transform the RNAseq read-count data to approximate a normal distribution. The transformation to normal distribution allows the use of existing methods such as the nearest shrinkage method (12,13) or Random Forest for classification. However, such transformation removes the count-based nature of the RNAseq read counts data, and hence, lacks the ability to fully preserve the strong mean-variance relationship that is otherwise useful for glioma classification and prediction (14,15). In order to appropriately model RNAseq read-count data, Negative Binomial (NB) and Poisson distributions are commonly used (16). The Poisson distribution is a single parameter distribution with its mean equals to its variance, which makes it rather restrictive. On the other hand, NB is similar to a Poisson distribution with an additional parameter called "dispersion" that allows the NB distribution to modify its variance without affecting the mean.
RNAseq uses high-throughput or next-generation sequencing technology (NGS) and has emerged as a novel alternative to microarray-based techniques for quantifying gene expression. The microarray technique is known to suffer from background noise. Gene expression level is measured as the relative frequency of the RNAseq reads that are mapped to one gene (17). RNAseq is a very sensitive technique that provides high resolution and a thorough understanding of the transcriptome and has revealed many novel gene structures.
RNAseq distribution requires an appropriate model that adapts and preserves the nature of RNAseq read counts data, and such classification models that preserve the nature of RNAseq are lacking in the traditional ML literature. The NB distribution is an appropriate choice to model such discrete reads counts data (16). Even though traditional ML tools that are developed based on NB are lacking, the choice of using NB distribution in differential gene expression and RNAseq analysis has been adapted by different studies in the literature such as in EdgeR (18)(19)(20), DESeq (21), and NBPSeq (22).
An example of a count-based classifier that fits a NB distribution is the Negative Binomial Linear Discriminant Analysis (NBLDA). NBLDA is a well-known classifier that is developed by fitting NB to RNAseq and the mean and dispersion parameter are estimated from the RNAseq data (23). A different type of classifier, known as VoomNSC, is developed based on the transformed count data. VoomNSC is a combination of Voom (an acronym for mean-variance modeling at the observational level) transformation (12) and the nearest shrunken centroids classifier (NSC) (24).
Consequently, the aim of this work is to implement a joint radiogenomics-NB model that predicts and classifies glioma molecular mutations following the 2016 World Health Organization's (WHO) updated guidelines for classification of tumors of the Central Nervous System (CNS) (including high grade and diffuse low-grade gliomas) (25). This work is critical especially when the RNAseq of some cases are unknown and a careful assessment is needed to avoid mischaracterization of lower grade gliomas. In this work, we utilize both volumetric features (radiomics) and RNAseq to implement and learn a radiogenomics-NB model. Then, the trained radiogenomics model is used to predict and classify the unknown RNAseq data. In the proposed model, a log-linear regression modeling is fitted to the estimated mean of the NB distribution and is linked with radiomics. We introduce this step to fuse the continuous radiomics data with the RNAseq count-based data without the need to transform RNAseq data into a normal distribution. Finally, we compare our radiogenomics-NB model performance with that of different genomics and radiogenomics state-of-theart methods in the literature.
The rest of the paper is organized as follows. A complete step-by-step mathematical derivation of the radiogenomics-NB model and parameters' estimations are presented in section Methodology. Section Experimental Results addresses the dataset used in this study, the data preparation, and the effect of using different numbers of differentially expressed genes in the radiogenomics-NB model. Furthermore, in section Experimental Results, a comparative analysis is discussed in which we compare the proposed radiogenomics-NB model's performance with different well-known methods in the literature. Moreover, in section Experimental Results, we investigate the effect of gender by developing a gender-specific radiogenomics-NB model for glioma molecular grading. Finally, the study's discussion is addressed in section Discussion.

METHODOLOGY
In this study, we propose a radiogenomics-NB method for glioma molecular grading and prediction. Figure 1 illustrates an overall flow diagram of the proposed radiogenomics-NB model. In Figure 1A, we fit the NB distribution to RNAseq read counts of the training dataset and estimate the model mean and dispersion parameter. Then, we use the estimated mean along with the predictor radiomics vector in a log-linear regression model to estimate the model regression coefficients. The dispersion parameter is estimated using the weighted likelihood empirical Bayes method (19). In Figure 1B, the estimated parameters of regression coefficients and the dispersion parameters along with the sample radiomics and its RNAseq read counts are utilized to predict the class label of a future test sample. A complete mathematical derivation of the radiogenomics-NB model is presented in the following subsection.

Prediction Using Negative Binomial Regression Model
To fuse radiomics with RNAseq read counts data in an NB model, the following parametrization is defined: Let C be the total number of classes, and I c ∈ (1, . . . , n c ) be the indices of samples in class c for c = 1, . . . , C. The examples of different classes include: IDH mutated vs. wildtype IDH (C = 2), 1p/19q codeletion: codeletion vs. non-codeletion (C = 2), Mutated ATRX vs. wildtype (C = 2). Let Y i = y i1 , y i2 , . . . , y iG be the RNAseq read counts training sample in the class label c and G is the total number of RNAseq. The purpose of this study is to predict the class label c of a future observation Y t using training samples associated with known class labels: p ( c|Y t ) ∝ p (Y t |c) p c , where p c is the probability of class c.
Using Bayes' rule, we have, where, p(Y i | c) is the pdf of the sample Y i in class c, and p c is the prior probability that one sample comes from class c. The pdf of class-specific c of RNAseq read counts of sample Y i and of RNAseq g is, ( In this parameterization, Y ig represents a count response of RNAseq, where µ igc represents the mean, φ g represents the dispersion parameter, E Y ig = µ igc , and Var Y ig = µ igc + µ igc φ 2 g . Note we assume that all RNAseq are independent of each other, so we have, Evaluating Equation (1) requires an estimation of p (Y i |c) and p c . The model in Equation (2) states that Y ig ∼ NB µ igc , φ g . We first estimate φ 1 , φ 2 , . . . , φ G , and µ i1c , µ i2c , . . . , µ iGc of all the training samples n c , and all RNAseq G. The mean is estimated as µ igc = s ic λ gc , where s ic is the size factor (26,27) which is used to scale RNAseq counts for the ith sample (in class c), λ gc is the total number of reads of RNAseq g across all samples in class c. For prior p c , we assume all classes are equally likely, p c = 1/C. Note that µ igc , s ic , and λ gc are estimated for each class c.
Next, plugging these estimates into Equation (2) and using the assumption of independent RNAseq, Equation (1) yields, The log-likelihood log p (Y i |c) is written as, Equation (5) can be written as, Rewriting Equation (6) yields, The proposed NB model of genomics relates to the radiomics (imaging features) X through the mean parameters µ igc (estimated mean of an ith sample and RNAseq g in class c).
We assume a log-linear regression model for estimating the mean µ igc in terms of the radiomics (imaging features) is given as follows: where X i is a p-dimensional of radiomics, β gc is a pdimensional vector of unknown regression coefficients (translate the relationship between X and Y through µ igc ). The estimation of β gc depends on class c and gene g of the ith sample. Hence, if there are two classes, we will need to estimate β g1 and β g2 (one from each class). Plugging Equations (8.a) into Equation (7), yields, Using the estimatedβ gc , andφ g from the training data, we classify a test observation Y t as follows, and,

Radiogenomics-NB Model Parameter Estimation
Estimating Dispersion φ g Using Weighted Likelihood Empirical Bayes Various methods for estimating the dispersion parameter are proposed in the literature. The EdgeR method applies a weighted conditional log-likelihood method to estimate the dispersion parameter (19). The weighted conditional log-likelihood (WL) for φ g is defined as a weighted combination of the individual (per-gene) likelihood l g φ g and common l C φ g likelihood: where α is the weight of l C φ g . In EdgeR,φ g is assumed to be normally distributed with means φ g and known variance τ 2 , and has the following hierarchical model: Under this hierarchical normal model, the maximum weighted conditional log-likelihood estimator is given as: Frontiers in Medicine | www.frontiersin.org where, and,

Computation of the Mean of RNAseq µ igc
The size factor s ic of sample i and class c is the total number of RNAseq read counts of that sample divided by the total number of all RNAseq read counts across all training samples (in class c). The size factor estimation is vital to account for the different sequencing depth (library size) that may be used to sequence different samples and is computed as follows: where, y igc is the RNAseq read count of sample i and RNAseq g in class c, and n c is the total number of samples in class c.
The mean µ igc of sample i and RNAseq g in class c is then estimated as µ igc = s ic λ gc , where λ gc is the total number of reads per RNAseq in class c, and is computed as follows: Using the estimated value of µ igc , the values of β gc are computed using equation 8.a as follows: The algorithm in Figure 2 illustrates the steps of estimating the different parameters in the radiogenomics-NB classification model.

Dataset
The dataset in this study consists of 108 pre-operative lower grade glioma (LGG) patients that are described in Menze et al. (28)

Data Preparation
In this study we first filter RNAseq read counts to remove RNAseq with very low value of read counts before performing any statistical analysis. RNAseq with very low read counts hold very little information because an RNAseq of biological importance needs to be expressed at some minimal level. We utilize a quantile filter (31) with a quantile threshold of 0.25. This step returns each RNAseq that has a mean across all samples higher than the defined quantile threshold of 0.25. Then, we reduce the number of RNAseq that are used in the radiogenomics-NB models, by utilizing EdgeR (18)(19)(20) to extract the differentially expressed RNAseq (DERs). DERs reflect the significance of a gene in a certain biological condition. In this study, we select the top 10, 20, 30, 50, 100, and 150 DERs (see Supplementary Table 1). Furthermore, we use eight volumetric radiomics features as illustrated in Table 1. ANOVA analysis for radiomics in Table 1 shows that feature numbers 1, 3, 4, 7, and 8 are significantly associated (ANOVA test, p < 0.05) with IDH mutations as illustrated in Figure 3A. Our analysis also indicates that feature number 2 is marginally associated (ANOVA test, p = 0.07) with 1p/19q codeletion. Furthermore, our analysis indicates that feature numbers 3, 4, 6, 7, and 8 are significantly associated (ANOVA test, p < 0.05) with ATRX mutations as illustrated in Figure 3B. Additionally, our analysis reveals that thresholding feature number 6 around the mean creates an ordinal feature that is significantly associated (ANOVA test, p < 0.05) with IDH mutations, 1p/19q codeletion, and ATRX mutations. Likewise, thresholding feature numbers 1, 3, 5, 7, and 8 around their means converts these features into ordinal features that are significantly associated (ANOVA test, p < 0.05) with IDH and ATRX mutations. Moreover, thresholding feature numbers 5, 6, 7, and 8 around their median converts these features into ordinal features that are significantly associated (ANOVA test, p < 0.05) with IDH and ATRX mutations.
Few other studies suggest that these volumetric imaging features and their ratios are associated with and predictive of several mutations in gliomas (32)(33)(34)(35).
The 108 LGG cases are randomly split into 80% training and 20% testing sets, and a balanced distribution of the target molecular alteration is ensured in the training and testing sets in each molecular classifier. The trained model classifier is developed using the training set. Model performance prediction is estimated and reported using the testing sets in terms of accuracy, balanced accuracy, F1 score, sensitivity, specificity, negative predictive value, and positive predictive value. The training set is utilized to build our radiogenomics-NB classifier as shown in steps 1-4 in Figure 2. The testing set is used to estimate the performance of the classifier as shown in steps a and b in  (38) repeat training and testing analysis for a specific number of times to ensure the robustness of the model performance. Consequently, in this work, we repeat the whole procedure 100 times independently for the 3 molecular alterations and then report the mean and standard deviation of the classifiers' performance using the testing sets.
Model performance parameters are computed based on the confusion matrix in Figure 4 as follows: Positive predictive value = TP TP + FP; Negative predictive value = TN FN + TN; Balanced Accuracy = Sensitivity + Specificity 2 , and (25) where TP is the true positive, TN is the true negative, FP is the false positive, and TN is the true negative.

Radiogenomics-NB Models Using Different Number of Differentially Expressed RNAs
In this section, we investigate the importance of using different numbers of DERs on the performance of the radiogenomics-NB model.
LGG radiogenomics-NB mutation prediction models are developed based on the top 10, 20, 30, 50, 100, and 150 DERs. The performance of the radiogenomics-NB IDH model using the top 10 DERs achieves slightly higher performance. However, such improvement is not statistically significant (ANOVA test, p > 0.05) when compared to the performance of the IDH models with the other number of DERs ( Figure 5A) except for negative predictive value (NPV) performance when using the top 20 DERs. Using the top 20 DERs in the IDH model achieves significantly worse NPV when compared to the NPV achieved using the top 10 DERs (ANOVA test, p < 0.05).
Radiogenomics-NB IDH model with the top 10 DERs (red line in Figure 5A) achieves an overall accuracy (Acc) of 0.92 ± 0.06, sensitivity (Sens) of 0.94 ± 0.07, specificity (Spec)  worse performance when compared to the performance of using the top 10 DERs (ANOVA test, p < 0.05). Using the top 10 DERs, the radiogenomics-NB codeletion model achieves an accuracy of 0.93 ± 0.06, a balanced accuracy of 0.90 ± 0.10, F1 score of 0.86 ± 0.14, a sensitivity of 0.84 ± 0.19, a specificity of 0.96 ± 0.04, an NPV of 0.95 ± 0.06, and a PPV of 0.90 ± 0.12, respectively. Radiogenomics-NB ATRX model also achieves similar performance (ANOVA test, p > 0.05) using the top 10, 20, and 30 DERs, even though the performance when using the top 10 DERs is slightly better as illustrated in Figure 5C. Using the top 10 DERs, the ATRX model achieves an accuracy of 0.85 ± 0.07, a balanced accuracy of 0.85 ± 0.07, an F1 score of 0.82 ± 0.08, a sensitivity of 0.86 ± 0.13, a specificity of 0.85 ± 0.09, an NPV of 0.91 ± 0.08, and a PPV of 0.80 ± 0.10, respectively. Figure 6 illustrates a graphical performance comparison between our radiogenomics-NB model with that of four different classifiers in the literature: NBLDA (23), VoomNSC (12, 13), RF-genomics where we first log-transformed (20) the RNAseq into a normal distribution, and RF-radiogenomics. Note that the number of DERs that we apply to develop these classifiers is 10 DERs. Moreover, when developing these classifiers, the 108 LGG cases are randomly split into 80% training and 20% testing sets, and balanced distribution is ensured when developing the different classifiers. The trained model classifier is developed using the training set, and 10-fold cross-validation is performed to identify the tuning parameters in the different classifiers. Model performance prediction is estimated and reported using the testing sets. Additionally, to ensure the robustness of the different classifiers' performance, we repeat the whole procedure 100 times independently and every training/testing set is utilized to develop and estimate the performance of each classifier.

Comparative Analysis
The NBLDA (23) classifier is developed by fitting NB to the top 10 DERs; then the mean and dispersion parameter are estimated from these DERs. In RF-genomics, the top 10 DERs of the training sets are first log-transformed into normal distribution and then fed into RF to build the RF-genomics classifier. In RFradiogenomics, radiomics (eight volumetric features described previously in section Data Preparation) are utilized with the is computed across 100 testing sets/splits. Y-axis represents the average performance of the different statistics on the X-axis. Different colors represent the radiogenomics-NB model with different numbers of DERs. The error bar represents one standard deviation. Asterisk "*" represents a statically significant difference between the performance achieved when using the top 10 DERs (in red) and using the number of DER where the star is located.
log-transformed DERs and then fed into RF to build the RFradiogenomics classifier. VoomNSC (12,24) is developed by first applying the Voom-based transformation on the 10 DERs and then applying the NSC classifier as illustrated in Zararsiz et al. (12) and Tibshirani et al. (24).
Comparing the performance of our radiogenomics-NB IDH model with that of NBLDA, RF-genomics, and VoomNSC, the radiogenomics-NB IDH significantly outperforms (ANOVA test, p < 0.05) these methods as shown in Figure 6A and Table 2. Additionally, our radiogenomics-NB IDH model significantly outperforms (ANOVA test, p < 0.05) the F1 score, balanced accuracy, and PPV performance of the RF-radiogenomics method whereas it achieves a similar (ANOVA test, p > 0.05) accuracy, sensitivity, and specificity. Our radiogenomics-NB IDH model archives an accuracy of 0.92 ± 0.06, a sensitivity of 0.94 ± 0.07, a specificity of 0.93 ± 0.18, an F1 score of 0.95 ± 0.04, and a balanced accuracy of 0.88 ± 0.09, respectively. The RF-radiogenomics-IDH model achieves an accuracy of 0.88 ± 0.17, a sensitivity of 0.93 ± 0.07, a specificity of 0.78 ± 0.16, an F1 score of 0.92 ± 0.06, and a balanced accuracy of 0.85 ± 0.08, respectively.
Our radiogenomics-NB codeletion model (Figure 6B and Table 3) performance is similar to NBLDA, RF-genomics, VoomNSC, and RF-radiogenomics models, except for the specificity and NPV performance when using RF-genomics and VoomNSC. The specificity and NPV of our model are significantly higher than those achieved by RF-genomics and VoomNSC. Our radiogenomics-NB codeletion model achieves an accuracy of 0.93 ± 0.06, a sensitivity of 0.84 ± 0.20, a specificity of 0.96 ± 0.5, an F1 score of 0.86 ± 0.14, and a balanced accuracy of 0.90 ± 0.10, respectively.

Gender-Specific Effect Analysis of Radiogenomics-NB
In our LGG dataset, IDH mutated patients, unlike IDH WT patients, have significantly longer survival (65.7 vs. 19.9 months, log-rank test p = 0.004). The association between IDH status and overall survival remains significant after stratifying for gender (likelihood ratio test p = 0.015). However, the association between 1p/19q codeletion and ATRX status and overall survival is not significant. Additionally, the chi-square test shows no significant association (p >0.05) between gender and IDH status, 1p/19q codeletion, and ATRX status. Table 5 shows patient IDH status, 1p/19q codeletion, and ATRX status distribution based on gender.
To explore the gender-specific effect in the performance of the radiogenomics-NB, we build two radiogenomics-NB models based on gender; male-specific radiogenomics-NB and female-specific radiogenomics-NB. Our analysis indicates that female-specific models significantly outperform (ANOVA test, p < 0.05) male-specific models as illustrated in Figure 7. In the radiogenomics-NB IDH, female-specific model achieves an accuracy of 0.93 ± 0.08, a sensitivity of 0.93 ± 0.09, a specificity of 0.91 ± 0.10, a PPV of 0.97 ± 0.05, an NPV of 0.83 ± 0.21, and a balanced accuracy of 0.92 ± 0.11, respectively. The male specific IDH model achieves an accuracy of 0.85 ± 0.08, a sensitivity of 0.97 ± 0.06, a specificity of 0.35 ± 0.33, a PPV of 0.86 ± 0.07, an NPV of 0.55 ± 0.48, and a balanced accuracy of 0.66 ± 0.17, respectively.

DISCUSSION
In this study, we propose a novel radiogenomics-NB model to fuse radiomics (imaging features) with RNAseq (genes) for glioma grading and prediction. NB distribution is appropriate for modeling RNAseq discrete read counts data and for preserving the count-based nature of this data. In the proposed radiogenomics-NB model, log-linear regression modeling is fitted to the estimated mean of the NB distribution and is linked with radiomics. We introduce this step to fuse the continuous radiomics data with the RNAseq count-based data without the need to transform the RNAseq data into a normal distribution.
The NB, unlike a Poisson distribution, has two parameters; the mean (e.g., the expected value of the RNAseq read counts data) and dispersion (e.g., a parameter that helps in capturing FIGURE 7 | Gender-based radiogenomics-NB models performance. (A) IDH mutations, (B) 1p/19q codeletion, and (C) ATRX mutations which are computed across 100 testing sets. The error bar represents one standard deviation. The asterisk * illustrates a significant difference between the two measurements. Y -axis represents the average performance of the different statistics on the X-axis. Different colors represent the female-and male-specific radiogenomics-NB models.
the variability of the RNAseq read counts). If the dispersion of NB is zero, the model reduces to Poisson distribution. In Poisson distribution, the mean is equal to the variance, which makes it rather restrictive. However, variation is usually observed in the real data of RNAseq counts data that the Poisson distribution cannot handle properly. On the other hand, NB has an additional parameter called the "dispersion" that allows the NB distribution of RNAseq counts data to modify its variance without affecting the mean. Thus, NB serves as a practical approximation to model RNAseq count data with variability different from its mean.
The mean of the proposed radiogenomics-NB model is estimated as the size factor multiplied by the total number of reads per RNAseq. Moreover, we utilize EdgeR to estimate the dispersion of the proposed radiogenomics-NB assuming RNAseq variability is assessed using the weighted conditional loglikelihood model. In the weighted conditional model, RNAseq counts data is assumed to have a distinct and individual dispersion for each RNAseq in addition to a common dispersion. Such an assumption can be more reliable when estimating the dispersion of real data of RNAseq counts data.
The performance evaluation of the proposed work indicates that linking simple, clinically feasible radiomics (i.e., tumor volumetric features) to RNAseq improves the performance of IDH and ATRX mutations prediction. The radiomics features utilized in the proposed radiogenomics-NB model that are described in Table 1 mainly depend on volumetric features. Our analysis shows that these features are associated with particular glioma mutations. This outcome supports previous studies that show the association between volumetric features and glioma mutations (32)(33)(34)(35). The efficacy of the proposed radiogenomics-NB model is further investigated using the top 10, 20, 30, 50, 100, and 150 DERs, respectively. Our analysis shows that the smaller the number of DERs (fewer than 30 DERs) utilized in radiogenomics-NB, the better is the radiogenomics-NB model performance. Our analyses indicate that using fewer than 30 DERs in our analysis offers the best performance (statically significant) in the radiogenomics-NB codeletion and ATRX prediction model. This suggests that using large numbers of DERs (more than 30) in the proposed radiogenomics-NB would over parametrize the dataset and create model fitting problems and thus degrade the performance.
Comparing our radiogenomics-NB model to NBLDA, RF-genomics, FR-radiogenomics, and VoomNSC, our model significantly outperforms NBLDA, RF-genomics, and VoomNSC for prediction of IDH and ATRX mutations. Our radiogenomics-NB model offers similar performance as NBLDA, RF-genomics, RF-radiogenomics, and VoomNSC models for prediction of 1p/19q codeletion. Specifically, for prediction of IDH mutations, while the proposed radiogenomics-NB model achieves significantly better balanced-accuracy, F1 score, and PPV than RF-radiogenomics, our model achieves similar accuracy, sensitivity, and specificity. Such results indicate the power of fusing radiomics and genomics data to develop radiogenomics models for classification and prediction models. The findings in this work indicate that the radiomics volumetric features may be vital for the prediction of IDH and ATRX mutations along with the genomics.
Different studies have revealed that gender is a significant factor in identifying cancer survival, prognosis, and treatment response (39)(40)(41). Hence, improved glioma molecular mutation prediction may require the development of gender-specific models. In this study, we explore the gender-specific effect on the radiogenomics-NB models. Our analysis reveals that IDH mutated patients remain significant after stratifying for gender, unlike 1p/19q codeletion and ATRX status. Moreover, our analysis indicates that no association is found between gender and the three specific mutations (IDH mutations, 1p/19q codeletion, and ATRX status) using the Chi-square test. This result is in agreement with the findings in Brat et al. (42), Li et al. (43), and Ebrahimi et al. (44). However, our gender-specific modeling shows that female-specific radiogenomics-NB models significantly outperform the male-specific radiogenomics-NB models for prediction of IDH status, 1p/19q codeletion, and ATRX status, respectively.
In conclusion, we present a glioma mutations radiogenomics-NB prediction model that preserves the count nature of RNAseq counts data in the NB model and utilizes radiomics to develop a complete and a better characterization prediction model of patient data. Our analysis shows the superiority of utilizing both genomics and clinically feasible radiomics data when compared to only genomics models. Use of tumor volumetrics can be more easily and reproducibly implemented in clinical practice compared to more complex radiomics metrics, such as higher order texture analysis features. Finally, this study shows the efficacy of volumetric radiomics features in the radiogenomics-NB model for glioma molecular characterization and prediction. This study is a first step toward implementing joint modeling of RNAseq and MRI patient data for glioma grading. However, further investigation is needed with a larger dataset with both RNAseq and full multimodality MRI dataset for each patient in a cohort. In the future, larger prospective studies may be needed to investigate specific radiomics features and their association with the different mutations and RNAseq read counts data for implementation into clinical workflow. Furthermore, it will be interesting to investigate the cause of superior performance of female-specific radiogenomics-NB models when compared to that of the male-specific radiogenomics-NB models for prediction of IDH status, 1p/19q codeletion, and ATRX status. Also, these models may be further investigated in treatment response and survival prediction in the future.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Old Dominion University IRB. The ethics committee waived the requirement of written informed consent for participation.