An Explainable Host Genetic Severity Predictor Model for COVID-19 Patients

1. Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has continued to pose a great threat to humanity ever since its first outbreak in late 2019. The SARS-CoV-2 viral strand causes the new coronavirus of 2019 popularly known as COVID-19 which has claimed millions of lives. The disease is widely characterized by a spectrum of clinical severity, suggesting a complex and highly dynamic host response in patients¹. Host (human) genetic variation associated with severity susceptibility or infection might provide clues to effective points to develop therapy or even preventive measures to intervene to develop medicine and vaccine against SARS-CoV-2 infection². Most especially, the scientific community is of kin interest that such findings provided by genetic human variations could give important clues where existing drugs may be repurposed for effective therapy against SARS-CoV-2 infection and life-threatening COVID-19 disease³. Also, we might be able to spot groups of individuals in the human population that might be at unusually high risk and need to be protected or might have innate protection against the SARS-CoV-2 infection^4–7. The SARS-CoV-2 genetic severity and susceptibility can also manifest themselves in rare genetic mutations which can cause healthy individuals to have a life-threatening response to COVID-19 disease⁶. Comorbidities such as diabetes, hepatitis, HIV, kidney-related problems, age, and gender have been observed in several clinical studies to have strong ties with patients’ severity and susceptibility to the disease^7,8. Some hosts are more susceptible to developing a severe disease probably due to modulated influence of genetics, environment, and risk factors.

There is a knowledge gap as to why the response to COVID-19 infection varies so much from patient to patient. The study of human genetics to diseases from several studies has pointed out some links to the severity of the disease among some groups of patients ⁹. For example, in some cases healthy and young patients with no prior existing medical conditions when exposed to the disease developed severe symptoms and some even succumbed to death from the disease. Emerging evidence suggests that asymptomatic patients mount a weaker immune response to the COVID-19 virus^10–12. There are some complex genetic interactions with the disease on the host side that can help to explain the variability in COVID-19 severity susceptibility and outcomes among patients¹³. Vital information as to why the disease differs greatly between people might lie in their DNA (e.g., variations in immune-related genes). Gene expression identifies patterns within human immune cells and may also play a key role in determining how the host immune system interacts with the virus. Examining genomes of patients who have a severe response to COVID-19 becomes necessary to further understand these complex interactions that are crucial to shed more light on understanding the biology of the disease, selecting drugs for repurposing – and knowing patients who are most at risk or providing some sorts of protection against the infection^14,15.

Most scholarly works so far utilized the Geno-Wide Association studies (GWS) approach to identify pointer chromosome loci, variants, and genes that linked COVID-19 susceptibility severity in patients^16,17. These studies have been vital in establishing crucial genetic variants linked to the disease, for example, the study of Q. Zhang et al.,¹⁸ used a candidate gene approach to identify patients with severe COVID-19 who have mutations in genes involved in the regulation of type I and III IFN immunity. Their studies found enrichment of these genes in patients and concluded that genetics may determine the clinical course of the infection. However, the GWS approach usually adopts more stringency measures for filtering genetic variants associated with a disease. For example, the threshold of small p-values is usually lowered to identify highly significant genetic variants with a disease¹⁹. Additionally, there is currently a lack of an organic model explaining how genetic factors concur with COVID-19 susceptibility in patients and possibly predict the genetic severity of the disease²⁰. At the time of writing this report, there is still a lack of Machine learning (ML) techniques available to experts working on Whole-exome sequencing (WES) datasets related to predicting genetic severity and providing explanations for the interaction of identified genetic variants linked to the severity of SARS-CoV-2 among patients. Explainable ML techniques that link identified genetic variant predictors can provide biological insights and interpretation that can assist experts in the field of medical sciences in understanding the complex genetic interactions that may lead to or hinder therapeutic approaches toward a trajectory to personalized medicine^21–23. Model interpretability approaches have been adopted in several studies to discover new knowledge and testable hypotheses ^24–26.

This study, therefore, sheds new light on genetic variants discovered from a previous study linked to the severity of COVID-19 in patients of European descent. This study provides a host genetic severity prediction and post-hoc model explanations using 16 identified candidate genetic variants and clinical covariates (age and gender) from a prior study using a 2000 cohort patients’ WES dataset²⁰. The post-hoc explanation of the HGSP model predictions using the ExplainerDashboard. The ExplainerDashboard is an open-source python library used by social scientists to explain global or local prediction results of their model. We adopt this approach to interpret plausible complex genetic interactions that may alongside patients’ clinical covariates, interplay with the host severity of the disease outcome. The HGSP model was developed by combining several trained decision tree-based models (Random Forest and XGBoost classifiers) across a 5-fold cross-validation (CV) splitting computational strategy we employed on the original problem dataset to improve the model stability abilities ²⁰. We further carried out domain interpretations of the genetic variants linked to the severity of the disease via an enrichment database (Transcription, Pathways, Ontologies, Diseases/Drugs, and Cell Types)²⁷. We also linked the genetic variant via association mapping studies using the Phenome-wide Association (PheWAS) technique by leveraging an open-source bioinformatics PheWAS tool called OpenTarget to associate the 16 identified genetic variants with disease traits that could lead to plausible clinical trajectories of the COVID-19 disease severity in patients ²⁸.

To this end, this study aimed to use the HGSP model developed from a previous 2000 cohort study of COVID-9 severity among European patients using their WES and covariates information ²⁰ for out-of-sample prediction (3000 cohort WES and covariates information) and importantly, to provide post-hoc model interpretations and explanations that linked disease traits and COVID-19 severity in patients using a genetic dataset.

The rest of this paper is organized as follows: In Section 2 we explore related works, Section 3 discusses the methods, Section 4 details the methods employed, Section 5 presents the results, Section 6 discusses the results, and Section 7 concludes the study.

2. Related Work

Several studies have established the roles of genetic risk factors and susceptibility to COVID-19 disease, for example, the study of Choudhary & Sreenivasulu et al.,⁸ observed that the host genetic factors, along with other risk factors may help determine susceptibility to respiratory tract infections. It is hypothesized that the Angiotensin Converting Enzyme 2 (ACE2) gene, encoding ACE2, is a genetic risk factor for SARS-CoV-2 infection and is required by the virus to enter cells ²⁹. Together with ACE2, Transmembrane Protease Serine 2 (TMPRSS2) and Dipeptidyl Peptidase-4 (DPP4) also play an important role in disease severity. According to Debnath, & Banerjee et al.,⁴ considering the implications of host genes in the entry and replication of SARS-CoV-2 and in mounting the host immune response, it showed that multiple genes might be crucially involved in shaping COVID-19 patients’ severity dynamics. Their study proposed three potential genetic gateways: variations within the ACE2 gene influencing the spatial transmission dynamics of the disease, the Human Leukocyte Antigen locus, a master regulator of immunity against infection which in turn seems to influence the susceptibility and severity of the disease in patients, and the genes regulating Toll-like receptor and complement pathways. This subsequently triggers the cytokine storm which induced exaggerated inflammatory pathways that seem to underlie the severity of COVID-19 in patients.

According to a study carried out by Lassau, et al..,³⁰ observed that the SARS-CoV-2 pandemic has continued to overwhelm hospitals’ intensive care units, and the need to identify predictors of disease severity is a priority. They conducted a study on 58 clinical and biological variables, and collected a chest CT scan dataset, from 1003 coronavirus-infected patients from two French hospitals. They trained a deep learning model based on CT scans to predict SARS-Cov-2 severity in patients. Constructing a multimodal AI-severity score that includes 5 clinical and biological variables (age, sex, oxygenation, urea, platelet) coupled with the deep learning model developed an AI-severity model that improves prognosis performance. The findings of their study demonstrated that AI approaches such as neural network analysis of CT scans can bring unique prognosis information to help clinicians. The work of Muhammad et al.,²⁴ utilizes supervised ML algorithms (Logistic Regression, Decision trees, Support Vector Machine, Naïve Bayes, and Artificial Neural Networks) using an epidemiological labeled dataset from positive and negative COVID-19 patients in Mexico. The results of their studies showed that the decision tree model outperformed all the other models in predicting the disease outcome. Yan et al., ³¹ utilize interpretable ML approaches to model COVID-19 patients’ mortality, and the results of their findings identified crucial predictive biomarkers of disease mortality. For this purpose, the ML tools selected three biomarkers that predict the mortality of individual patients more than 10 days in advance with more than 90% accuracy: lactic dehydrogenase (LDH), lymphocyte, and high-sensitivity C-reactive protein (hs-CRP). Relatively high levels of LDH alone seem to play a crucial role in distinguishing most cases that require immediate medical attention.

Patel, Kher, et al., ¹² uses socio-demographic data, clinical data, and blood panel profile data to develop ML algorithms for predicting the need for intensive care and mechanical ventilation. The findings of their study established that the three categories of information are crucial for intensive care and mechanical ventilation allocation to patients by healthcare systems in planning for surge capacity for COVID-19. The study of Doewes, & Sharma³², combined an ensemble genetic algorithm and ML classifier to filter relevant features and perform classification to diagnose COVID-19 through patients’ blood samples. Their findings showed that ML techniques can augment current medical practices and tools, and improve modern methods deployed in the fight against the disease. The work of Albadr, et al.,³³ employs an Optimised Genetic Algorithm-Extreme Learning Machine (OGA-ELM) with three selection criteria (i.e., random, K-tournament, and roulette wheel) to detect COVID-19 using X-ray images. The results from their findings showed that the OGA-ELM can be used to achieve 100.00% accuracy with fast computation time. This demonstrated that OGA-ELM is an efficient method for COVID-19 detection using chest X-ray images.

The study by Marcos, Belhassen-GarcÍa, et al.,³⁴ demonstrated that efficient and early triage of hospitalized Covid-19 patients is crucial to detect those with a higher risk of severe Covid-19 infection. Their study trained, validated, and externally tested an ML model to early identify patients who will die or require mechanical ventilation during hospitalization from clinical and laboratory features obtained at admission. Their findings strongly supported evidence that ML models could be useful in medical emergencies such as pandemics in determining patients for hospital admission and predicting the risk of severe Covid-19 infection. Different ML techniques have been utilized by researchers in response to mitigating the SARS-Cov-2 virus threat since the onset of the pandemic some of which have been useful in the diagnosis of COVID-19 and prediction of mortality risk and severity, using readily available clinical and laboratory data^7,30. Fallerini, et al.³⁵ studied common and rare variants from WES data of about 4000 SARS-CoV-2-positive patients. This was utilized to define an interpretable machine-learning model for predicting COVID-19 severity. The variants were converted into separate sets of Boolean features, depending on the absence or the presence of variants in each gene. They developed an ensemble model of LASSO logistic regression which was used to identify the most informative Boolean features concerning the genetic bases of severity. Their study formed the basis of our study, however, in our study, we took a different approach to the development, interpretations, and explanations of the HGSP model. First, we used a computational strategy of a 5-fold CV to split the clinical phenotype information of the patients (i.e., each of the folds in the 5-fold CV has 80% training and 20% testing sets) and screened the 80% training sets in each fold for genetic variants (from our prior study ¹⁹) linked with the disease using an odd-ratio statistics approach other than the GWAS. We then remapped the identified significant variants for each of the folds for the 20% unscreened testing sets in each fold. We developed our count matrices (training and testing for each fold) following the identified variants and employed several traditional ML algorithms (Logistic Regression, Support Vector, Random Forest, and XGBoost classifiers) to identify genetic variants via feature importance from the decision tree models (i.e., Random Forest and XGBoost classifiers) that provided the highest support (i.e., non-zero feature importance weighted scores across the 5-fold). 16 genetic variants and the clinical covariates (age and gender) were identified with the highest support, and we retrain the models using these features. An ensemble voting classifier was used to combine the 5-fold decision tree models to develop the HGSP model. We validated the HGSP model on an external dataset and provided further domain interpretations and post-hoc explanations.

3. Materials and Methods

Data Resources and Description

The genetic and clinical data type we used for this study was extracted from an observational study dataset coordinated by the GEN-COVID Multicenter Study group led by the University of Siena with an estimated enrolment of 2000 participants which constitutes the basis of the sample from which the genetic information was obtained (https://clinicaltrials.gov/ct2/show/NCT04549831). The initial 2000 cohort WES and covariates (age and gender) dataset were used to train and develop the HGSP model. The HGSP model was developed based on 16 identified candidate genetic variants and clinical covariates (age and gender)²⁰. The out-of-sample dataset used for external prediction in this study was a follow-up dataset information provided by the same GEN-COVID group to validate the HGSP model and this formed the basis of this study. In this study, we maintained the same data preparation procedures we used for the training dataset study ¹⁹²⁰. To perform the downstream analyses (external model validation and post-hoc interpretations), we considered the two classification grading schemes: adjusted-by-age and unadjusted-by-age. The adjusted-by-age grading classification scheme refers to a refinement made on the grading classification based on an ordinal logistic model which uses age as an input feature for sex-stratified patients^3,11,35, while the unadjusted-by-age grading classification scheme was not refined. To build the feature matrices for the external model validation, we curated the allele frequencies from the WES dataset (out-of-sample) for the 16 identified genetic variants together with their clinical covariates.

We binarized the patients’ COVID-19 outcome severity grading classifications: gradings 5, 4, and 3 were grouped as “severe” and coded as 1; grading 0 was grouped as “asymptomatic” and coded as 0. We excluded patients’ information whose severity grading was classified as 1 & 2. This exclusion was done purposefully to minimize noise signals during the filtering of genetic variants that are linked with the disease severity or protection in patients. We further refined the grading classification based on an ordinal logistic model which uses age as an input feature for sex-stratified patients ³⁶ and we retained only those patients whose grading classification was concordant with the one adjusted by age. The age distribution of the participants ranges from 18 years and older (Adult, Older Adult). The Gender this study considered were all sexes. Patients’ genotypic information contains a reference (Ref) or alternative (Alt) alleles in either severe or control groups which were defined by employing an additive model, whereby homozygous genotype (1/1) has twice the risk (or protection) of the heterozygous type (0/1 or 1/0).

Further details and descriptions of the WES dataset and preprocessing can be found in our prior study of the 2000 cohort dataset²⁰. A total of 618 patients met the filtering and selection criteria from the out-of-sample dataset and were therefore used for the external model validation. The feature matrix for the external validation of the HGSP model contained 18 features (16 genetic variants and covariates). Each of the 16 variants that constituted the feature matrix has allelic frequency counts for each patient’s genotype information assigned (0/0 counted as 0, 0/1 or 1/0 counted as 1, and 1/1 counted as 2). The covariates were further merged into the feature matrix from the patients’ phenotype information. The outcome variable (grouping) is binarized as 0 and 1, where 0 represents asymptomatic and 1 depicts the severity of the disease as shown in figure 1.

• Download figure
• Open in new tab

Figure 1:

Classification of patients’ SARS-CoV-2 severity adjusted by age grading scheme. We considered the most severe cases patients in groups 3, 4, and 5 versus control cases patients in group 0.

• Download figure
• Open in new tab

Figure 2:

Patient phenotype information (adjusted-by-age grading scheme) in the follow-up WES dataset. We employed a four-scenario approach to validate the HGSP model using the out-of-sample dataset. The baseline scenario was mainly focused on validating the model using the adjusted-by-age grading classification scheme feature matrix. While scenarios 1 – 3 focused on validating the model using the unadjusted-by-age grading classification scheme.

Domain Interpretation

4. Results

We present results from a validation of the trained HGSP model we developed using the 16 identified genetic variants and clinical covariates (age and gender) considering the adjusted-by-age grading scheme (see Fig. 1). We further provided post-hoc model interpretations for the HGSP model prediction outcome considering the unadjusted-by-age grading scheme.

Results of Analysis Post-hoc HGSP Model Interpretations and Explanation

Next, we performed the post-hoc model agnostic interpretations and explanations at a local level for the HGSP predictions. We employed the ExplainerDashboard interpretation and explanation approach focusing on the SHAP dependence plots, and feature importance plots, results to further shed new light on understanding complex interactions of genetics with clinical covariates. Here we seek to unravel hidden insights such as patients whose COVID-19 severity predictions are not driven by covariates (age and gender) but as a result of some complex genetic interactions from the 16 identified consistent features. The HGSP does not perform basic EDA approaches such as descriptive statistics summary and bar or histogram plots. The ExplainerDashboard provides options to view different post-hoc interpretation results. The ExplainerDashboard results can be saved in PDF, HTML, JPEG, or PNG file formats. Displayed in Figure 5(a) – 5(e) are the results of the HGSP model from the ExplainerDashboard interactive interface.

• Download figure
• Open in new tab

Figure 3:

Proposed methodological design of COVID-19 HGSP model post-hoc explanations. We first load the saved HGSP model developed by training fully supported variants and covariates (age and gender) identified from a simple stratified 5-fold CV splitting strategy adopted from the 2000 cohort dataset study. We proceeded to validate the model on the out-of-sample dataset. Our method used the HGSP for prediction purposes and employed the ExplainerDashboard python library to provide a post-hoc model explanation of the prediction outcome at an individualistic level using the adjusted by-age information (baseline scenario). Note full details on how the HGSP model was developed and trained are provided in a study by Onoja et al. ¹⁹.

• Download figure
• Open in new tab

Figure 4:

HGSP model performance considering out-of-sample model validation in all case studies.

• Download figure
• Open in new tab

Figure 5(a): ExplainerDashboard displayed SHAP feature importance plot.

The plot shows the features sorted from most important to least important. The features were sorted based on the absolute SHAP values (average absolute impact of the 18 features on the final prediction outcome). The features can also be shuffled and sorted based on their permutation importance.

• Download figure
• Open in new tab

Figure 5(b): ExplainerDashboard displayed Classification Stats.

In figure 5(b) we show the ExplainerDashboard displayed the classification statistics results of the HGSP prediction: confusion matrix plot; precision plot and classification plot. The observations were binned together in a group of roughly equal predicted probabilities and the percentage of positives is calculated for each bin. A perfectly calibrated model would show a straight line from the bottom left corner to the right corner. A strong model would classify most observations correctly and close to 0% or 100% probability the classification plot displayed the fraction of each class above and below the probability cut-off of 0.50: the ROC curve performance of the ensemble model on an external prediction dataset.

• Download figure
• Open in new tab

Figure 5(c): ExplainerDashboard displayed Classification Stats

The plots displayed the Precision-Recall Area Under Curve (AUC) performance on an external follow-up cohort dataset; the lift curve plot is used to depict the percentage of positive classes when one selects only observations with a score above the cut-off Vs selecting the observations randomly. It aimed to help us evaluate how much better our developed ensemble voting classifier is than a random (the lift).

• Download figure
• Open in new tab

Figure 5(d): ExplainerDashboard Individual predictions plot

The dialogue box is pulled down to select a sample_ID directly by choosing it from the dropdown list or hitting the random sample_ID button to randomly select a sample_ID that fits the constraints. This aimed to help us assess in general the false positives and false-negative rates of our prediction. The doughnut prediction plot shows the predicted probability for each grouping label for the selected sample_ID of interest.

• Download figure
• Open in new tab

Figure 5(e): ExplainerDashboard Feature dependence plot

The SHAP dependence plot displays the relationship between feature values and SHAP values. This aimed to allow us to investigate the general relationship between feature value and impact on the prediction. One can ascertain whether the model uses features as expected or uses the plots to learn more about the relationships that the model has learned between the input feature and the predicted outcome.

Domain Interpretations and Explanations

We linked the genetic variants used to validate the HGSP model and ExplainerDashboard with the OpenTarget genetics and Enrichr bioinformatic web-based tools. The PheWAS results from OpenTarget genetics are presented in table 2 while Figure 6 presented a snapshot interface of the Enrichr results for the top 15 genetic variants.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2:

Feature importance description associated with PheWAS analysis

• Download figure
• Open in new tab

Figure 6: Snapshot of Enrichr web-based results linking the genetic variants used for the HGSP Model.

The detailed results of the Enrichr domain interpretations can be found using the link below: https://maayanlab.cloud/Enrichr/enrich?dataset=26f3365c99e0255115dd818c11aba294#

Next, we presented the results of the SHAP value feature importance from the ExplainerDashboard visualizing it via the hierarchical clustering heatmap (See Fig. 7).

• Download figure
• Open in new tab

Figure 7: ExplainerDashboard clustering heatmap visualizations of SHAP features importance output.

The covariates and variants (Fully supported variants) at the local explanations level we visualized the SHAP feature importance output for plausible interactions interplaying with COVID-19 severity predictions. The Shapley values for each feature’s important interpretations and explanations range from negative to positive. Negative values were colored blue while positive were colored red. A positive value means the feature is pushing the prediction output in a forward or positive direction while a negative value means the feature is pushing the output backward. Meaning features with positive pull force will favor grouping 1 (severe) while features with backward pull force (negative) will favor grouping 0 (asymptomatic).

5. Discussion

There is still a lot we are yet to unravel when it comes to the severity manifested by different SARS-CoV-2 patients. For example, why are certain patients even though not advanced in age and with no comorbidity susceptible severity to the disease while others are not? The significance of our study is three folds, first, we identified 16 candidate genetic variants that are likely determinants of COVID-19 severity outcomes in patients from a prior study we carried out using the 2000 cohort WES and clinical datasets. We further developed by combining several traditional interpretable ML algorithm models (decision tree-based models – Random Forest and XGBoost classifiers across a simple stratified trained 5-fold split CVs) to form the HGSP model and carried out an external validation using a follow-up dataset. Secondly, we utilized the explainer dashboard open-source python library to carry out a post-hoc explanation of our model predictions particularly considering the out-of-sample dataset. We employed the domain knowledge of the Phenome-wide Association technique by leveraging the OpenTarget web tool to associate the 16 identified genetic variants with disease traits that could lead to plausible clinical trajectories of the COVID-19 disease in patients. Lastly, we used the hierarchical clustering visualization heatmap to unravel hidden insights from the SHAP feature importance values from the explainerdashboard.

From the heatmap (see Fig. 7) the covariates of age and gender were separated along two lines of magnitude impacts (positive and negative directions). This means they can be a push forward (severe) or pull back (asymptomatic) to the model severity prediction. This is in line with existing findings from the literature that the male gender is more at risk of severe COVID-19 disease than the female gender. More so, research has shown that as one advance in age, they are more likely at risk of severe COVID-19 than younger patients.

Visualizing the SHAP value feature importance via a hierarchically clustered heatmap helps us to further understand the directionality contributions of the features at an individualistic level for each patient’s severity predictions. For example, the variants PLEC & PCSK5 (see Fig. 7) are strongly associated with abnormalities of breathing, symptoms, and signs involving the circulatory and respiratory systems. Most of the patients’ genetic severity prediction is strongly linked to age and gender, however, there is a small portion of patients whose severity predictions are topmost contributed by genetic variants; for example, patients with sample IDs S46_hg38, COV3908-1549_hg38, COV5958-2221_hg38, COV6351-2318_hg38, COV6807-2447_hg38, COV7658-2744_hg38, COV7878-2817_hg38 and COV8603-3159_hg38, whose severity contribution is coming from the variant PLEC, and not gender or age covariates. The variant PLEC has been linked with abnormalities in breathing phenotype traits and may be a suspected candidate disease variant whose interactions may interplay with the disease outcome in the identified patients. This implies that not all the patients’ severity is purely driven by the covariates’ age and gender, some of which come from the genetic interactions of the host with the disease. This is in line with our previous study which utilizes the Principal Component Analysis (PCA) on the training dataset and identified a small cluster group of patients (about 29 in number with 98% homogeneous COVID-19 severity). When these patients’ gene lists were further investigated for biological implications via the pathway enrichment, biological processes such as the JAK-STAT signaling pathway, Cytokine-cytokine receptor interaction, and Interleukin-6 family signaling were detected. This revealed that genetics may have inter-play with the severity of the SARS-CoV-2 virus of patients in this cluster grouping.

Also, the HGSP model is further deployed as a web application to assist experts working with WES datasets of COVID-19 patients and seek to evaluate their model based on the 16 identified genetic variants and clinical covariates. For example, an individual patient (see Fig. 5) can zoom in to understand how genetics is interacting with clinical covariates in predicting the severity outcome of the disease. The patient even though young and possibly healthy could be at susceptible risk to the disease if exposed to the virus. However, it is worth noting that the HGSP model tends to favor severity outcome predictions compared to the asymptomatic cases, partially due to inherent imbalanced class distributions in the problem dataset. This reflects the reality that during the pandemic’s first and second waves, there was a high surge of patients infected by the disease which led to a high influx of admitted severe patients and the collapse of healthcare facilities in many countries. The HGSP model, therefore, could serve as a good decision tool to lessen healthcare cost burdens by decision-makers to identify patients with true COVID-19 severity. Also, the HGSP model can be used to explain the rationale behind the patients’ severity prediction outcomes as to why they were selected or not. Developing a user-friendly explanation ML system can augment the efforts of healthcare decision-makers and clinicians to further build trustworthy and explainable models that will promote efforts toward personalized medicine.

6. Conclusion

Providing explanations and insights to the 16 identified genetic candidate variants we used to develop the HGSP model and validated it on an external follow-up dataset is crucial and formed the major contribution of this study. The findings from this study shed more insights into our understanding of the complex genetic interactions that the 16 identified genetic variants may be interplaying in the host genetic severity of COVID-19 disease. Most of the patients’ severity has been contributed by age, gender covariates, and other modulated factors such as comorbidities, however, there is a small portion of the patients whose severity is purely driven by genetics. This group of individual patients spotted in the population might be at unusually high risk and need to be protected against the SARS-CoV-2 infection. In the future, we seek to investigate the WES patients’ dataset using statistical and bioinformatic tools such as linear regression SKAT-O analysis, polygenic variant scoring, and deep neural knowledge prior networks to provide domain experts with biologically inspired interpretability abilities for robust data-driven solutions.

Data and code availability

We confirmed that all materials used for this study are readily available from the authors upon request. Further information, sample feature count matrix, codes, and figures are available in the Supplementary Information. The code implementation is available at https://github.com/raimondilab/COVID-19-severity-host-genetic-predictor-model-explanation Enrichr results: https://maayanlab.cloud/Enrichr/enrich?dataset=26f3365c99e0255115dd818c11aba294#