Exome-Wide Association Study of candidate pharmacogenes on DOAC haemorrhage

1. Introduction

Direct oral anticoagulants (DOACs) continue to first-line guideline-recommended medications for prevention of stroke in patients with non-valvular atrial fibrillation (AF) [1] and are the cornerstone of treatment for patients with venous thromboembolism [2]. Traditional drugs of the group of vitamin K antagonists (warfarin) and low-molecular-weight heparins have been the gold standard of anticoagulant therapy for decades. However, they still present challenges such as difficulties in drug dosing, a narrow therapeutic range, the need for continuous monitoring of INR values, and a high frequency of drug interactions. Currently, physicians are increasingly transitioning to the use of DOACs that inhibit factor IIa (dabigatran etexilate) or factor Xa (apixaban, rivaroxaban, edoxaban). The ease of use of DOACs compared with vitamin K antagonists has ushered in a new era of anticoagulants. DOACs are administered in fixed doses, do not require routine monitoring of INR values, their bioavailability is less influenced by food intake, and they have limited drug interactions. These benefits have led to an increase in DOACs consumption. In England, for example, DOACs accounted for 61.8% of all oral anticoagulant prescriptions in 2019, up from 16.4% of all prescriptions in 2015 [3]. In Germany, while DOACs accounted for less than 0.1% of anticoagulant prescriptions by 2011, by 2016 they were already at 49.9% [4]. In the US, the trend is similar: the proportion of people receiving DOACs under the Medicare Part D program increased from 7.4% in 2011 to 66.8% in 2019, from 200,000 to 3.5 million [5]. In the United States, apixaban was among the top 3 prescription drugs prescribed by the end of 2022 [6].

High interindividual variability in drug response, which can result in hemorrhagic or thromboembolic events, continues to be a significant concern associated with DOAC medications. Consequently, bleeding is recognized as a critical adverse reaction to anticoagulant treatment. Approximately 2-4% of patients receiving anticoagulant therapy with DOACs encounter major bleeding each year, with an additional 10-12% experiencing clinically relevant non-major bleeding [7]. The increasing utilization of DOACs has been inevitably linked to a rise in the occurrence of DOAC-related bleeding events. In the United States, the incidence of DOAC-related bleeding events showed an average annual increase of 27.9% between 2016 and 2020 [8]. In the UK, a study spanning from 2011 to 2016 revealed that for every additional 10% of DOAC prescriptions out of the total anticoagulant prescriptions, there was a 0.9% rise in hemorrhagic complications. Moreover, there were 4929 additional hospitalizations due to bleeding events associated with the initiation of active DOAC use during this period [9]. Such statistics may suggest that there is a continued need for efforts to enhance the appropriate prescription and monitoring of anticoagulants.

The variability in drug response to DOAC and the corresponding risk of adverse reactions are influenced by various clinical and demographic factors in patients, including sex, age, comorbidities, renal and hepatic function, concurrent medication use, among others. Recently, the role of patient genetics in determining the safety profile of DOAC utilization has gained significant attention. Findings from pharmacogenetic investigations indicate potential associations between genetic polymorphisms in enzymes involved in biotransformation and drug metabolism, and the response to DOAC, particularly for variants in CES1, ABCB1, CYP3A4, CYP3A5, and ABCG2 [10,11,12,13]. There is insufficient conclusive evidence linking specific genetic polymorphisms to an increased risk of bleeding. Further research is needed to identify new markers that may be associated with bleeding risk.

The objective of the current investigation was to identify novel pharmacogenetic biomarkers associated with the risk of bleeding in patients with non-valvular AF using DOACs, specifically apixaban and rivaroxaban.

2. Materials and Methods

2.1. Ethics Statement

This research was approved by the local ethical committee of the Russian Medical Academy of Continuing Professional Education during meetings held on November 25, 2020 (protocol No. 16), October 20, 2021 (protocol No. 12), and October 25, 2021 (protocol No. 15). The study was carried out in compliance with the Declaration of Helsinki. Written informed consent was obtained from all participants involved in the study.

2.2. Patients

The research involved 198 patients who had a verified diagnosis of non-valvular AF along with chronic kidney disease at different stages. Among them, 106 patients were administered apixaban, while 92 patients were given rivaroxaban.

The study sample comprised 97 patients with hemorrhagic complications and 99 patients without such complications. The participants were recruited from the clinical facilities of State Budgetary Institution of Health Care ‘Hospital of Veterans No. 2 of DZM’ and State Budgetary Institution of Health Care ‘S.S. Yudin Municipal Clinical Hospital of DZM’. A detailed description of the demographic, clinical, and laboratory parameters of the patients can be found in Table 1.

2.3. Measurement of Plasma Concentrations of Apixaban and Rivaroxaban

To determine the concentration of the drug, venous blood specimens were gathered using vacuum tubes containing EDTA-K3 anticoagulant. Subsequently, the blood samples were centrifuged at 3000 rpm for 15 minutes to isolate the plasma. The separated plasma was then divided into Eppendorf tubes and stored at –70 °C until further analysis.

The quantification of apixaban and rivaroxaban levels in blood plasma samples was conducted through high-performance liquid chromatography (HPLC) using an Agilent 1200 chromatograph equipped with a four-channel pump, mobile phase degasser, and chromatographic column thermostat. Detection was performed with an Agilent TripleQuad LC/MS 6410 mass spectrometer (triple quadrupole type).

Sample preparation involved the precipitation of blood plasma proteins. Initially, plasma samples were thawed at room temperature. Subsequently, 100 μl of plasma was dispensed into Eppendorf plastic tubes, followed by the addition of 250 μl of a 9:1 mixture of methanol and 0.1% hydrochloric acid (HCl). The mixture was vortexed, allowed to stand for 10 minutes, and then vortexed again. The resulting samples were centrifuged at 10,000 rpm for 10 minutes. The supernatant was then transferred to chromatographic vials and loaded onto the autosampler of the chromatograph.

An Agilent Polaris 3 C18-A column (length 50 mm; inner diameter 3.0 mm; particle size 3.0 μm) was utilized for the separation process. The column temperature was maintained at 40 °C. The mobile phase consisted of two constituents: solution ‘A’ (1 ml of concentrated formic acid diluted with deionized water to a total volume of 1 L) and solution ‘B’ (1 ml of concentrated formic acid diluted with acetonitrile to a total volume of 1 L). The chromatographic separation was conducted using a gradient elution method. The blending of the mobile phase components during the analysis followed the program outlined in Table 2.

The volume of the injected sample was 5 μl, and the analysis was conducted within a 10-minute timeframe. The spectra of apixaban and rivaroxaban were obtained using multiple reaction monitoring with electrospray ionization in positive ionization mode. The atomizer gas pressure was set at 35 psi, with a drying gas flow rate of 11 L/min at a temperature of 350°C. The fragmentation voltage was maintained at 135 V, and the collision cell voltage was set to 25 V.

2.4. CYP3A4 Phenotyping

To determine the concentration of the drug, venous blood specimens were gathered using vacuum

A CYP3A4 activity was determined by estimating the ratio of 6-β-hydroxycortisol to cortisol concentrations in the patient urine collected in the morning. Cortisol is a specific CYP3A4 substrate. By calculating the metabolic ratio of the concentrations of cortisol and its metabolite 6-β-hydroxycortisol, the activity of CYP3A4 is determined: high values of the ratio indicate a high activity of the isoenzyme, while low values indicate a low activity. The methodology for determining a CYP3A4 activity is also generally accepted [14].

Cortisol and its metabolite were determined by HPLC with a mass spectrometric detection. Agilent 1200 liquid chromatograph (Agilent Technologies Inc., USA, 2008) and AgilentTripleQuad LC/MS 6410 mass spectrometer were used. The results were processed using Agilient MassHunter Workstation Software LC/MS Data Acquisition for 6400 Series Triple Quadrupole (version B.08.02).

To perform a chromatographic determination, the sample preparation technique and conditions of chromatographic analysis presented in the work by Smirnov V.V. et al., were used [14].

2.5. Exome SEQUENCING

Venous blood samples were obtained from patients using vacuum tubes containing a minimum of 4 ml of EDTA-K3 anticoagulant. Genomic DNA was extracted from venous blood samples utilizing the QIAamp DNA Blood Mini Kit (Qiagen). Subsequently, libraries were generated from 500 ng of genomic DNA using the MGI Easy Universal DNA Library Prep Set (MGI Tech) following the manufacturer’s instructions. The DNA fragmentation process was carried out through ultrasonication on a Covaris S-220 instrument, resulting in an average fragment size of 250 bp. Enrichment of DNA libraries was achieved through pre-pooling as described in [15], employing SureSelect Human All Exon v8 probes (Agilent Technologies) designed to cover the entire human exome. Quantification of DNA and libraries was conducted using a Qubit Flex instrument with the dsDNA HS Assay Kit in accordance with the manufacturer’s guidelines (Thermo Fisher). Subsequent quality assessment of the libraries was performed on a Bioanalyzer 2100 utilizing a High Sensitivity DNA kit (Agilent Technologies) following the manufacturer’s protocol.

The libraries were subsequently circularized and subjected to paired-end sequencing on the DNBSEQ G-400 platform utilizing the DNBSEQ-G400RS High-throughput Sequencing Set PE100 following the protocol provided by the manufacturer (MGI Tech), achieving an average coverage of 100×. FastQ format files were produced through the application of the manufacturer’s basecallLite software (MGI Tech).

Bioinformatic quality control and quality correction of the resulting sequencing data were conducted utilizing FastQC v0.11.9 [16] and bbduk v38.96 software [17]. Subsequent bioinformatic processing of the sequencing data for each sample involved aligning to the GRCh38 human reference genome using bwa-mem2 v2.2.1 [18] and SAMtools v1.9 software [19], acquiring quality metrics for exome enrichment with Picard v2.22.4 [20], collating variants with bcftools v1.9 [21] and deepvariant v1.5.0 software [22], annotating variants using AnnoVar [23], Intervar v2.2.2 24] software, and employing several in-house scripts to enhance the quality of the final variant annotation files. Following the bioinformatic pipeline, MultiQC v1.16 [25] software was executed as a final step for quality control.

2.6. Exome-Wide Association Study (EWAS)

The study cohort comprised 97 patients with bleeding and 99 without bleeding, totaling 461,938 unique polymorphic markers with a genotyping quality of 99.47%. Genetic variants were filtered using PLINK v1.90b6.24 [26]. Multi-allelic variants were treated as false heterozygotes (--vcf-half-call r). Markers genotyped in less than 98% of samples were excluded, resulting in the removal of 16,542 markers. Samples with a genotyping quality below 98% were also excluded. Variants on the X and Y chromosomes were removed, and those with a minor allele frequency below 0.2 were filtered out, leaving 66,720 variants. A Hardy-Weinberg equilibrium test was conducted on 60,142 markers (p<0.00001), and the final genotyping quality was 99.99%. Association analysis was performed using logistic regression with additive, dominant, and recessive models, with covariates including the first principal components or drug. The resulting p-values were adjusted for multiple testing using Bonferroni corrections (P_adj<0.05).

We used PRSice-2 [27] with parameters --stat OR and --binary-target F to generate polygenic risk scores (PRS) for predicting the occurrence of bleeding events, with no file for --pheno. Subsequently, we conducted additional quality control measures and eliminated duplicates, resulting in a final set of 11,938 variants.

To exclude hereditary forms of bleeding we employed the HP:0001892 and HP:0000790 panels [28]. The diplotyping of pharmacogenes was conducted utilizing the Panno resource, with a focus on the European population [29].

2.7. Statistical Analysis

Statistical analysis and graphing of the data were conducted utilizing the R programming language (RStudio Pro 2023.09.1) and the GraphPad Prism 8.0.1 software (GraphPad, La Jolla, CA, USA). Any missing values were imputed using the amelia package [30].

3. Results

3.1. Clinical Features of the Individuals

Table S1 presents the clinical characteristics of the patients. Boruta’s algorithm [31] was utilized to assess the significance of the predictors determining C_{ss min}/D (ng/mL/mg). Due to the replacement of two missing values with negative values (–5.173 and –1.724) by amelia, the natural logarithm was calculated for each C_{ss min}/D value and 1 was added. Body mass index (BMI), the medications themselves, the presence of a history of coronary heart disease (CHD) and chronic kidney disease (CKD), and the ratio of 6-β-hydroxycortisol to urinary cortisol concentrations were identified as the most crucial factors for medication. Nevertheless, variations were observed among the drug groups: for apixaban, the concentration of urinary 6-β-hydroxycortisol (ng/mL) and its ratio to cortisol were significant, while for rivaroxaban, in addition to these parameters, age, the presence of CHD, and group membership (bleeding and non-bleeding) were also found to be important (Figure A1).

Interestingly, the concentration of urinary 6-β-hydroxycortisol and urinary cortisol (ng/mL) emerged as predictors of bleeding development, with the addition of the stage of CKD (Figure A2). In this instance, we utilized the unadjusted minimum concentration at steady state (C_{ss min}/D) value.

Utilizing the R package ‘pROC’ [32], an additional analysis was conducted to assess the predictive capability of C_{ss min} and C_{ss min}/D (Figure 1). The results indicated statistical insignificance, with low sensitivity and specificity values, highlighting the necessity of incorporating other variables that impact the risk of bleeding.

The residual equilibrium concentration adjusted for daily dose (C_{ss min}/D) did not exhibit a normal distribution, as depicted in Figure 2. A similar analysis was conducted for C_{ss min} and blood platelet counts, as illustrated in Figure A3.

Factorial analysis of variance was utilized to examine the differences between patients with and without bleeding in terms of C_{ss min}/D. The results indicated that there were no statistically significant variances (F(1,192)=0.0367, p=0.5454). However, when considering drug intake within groups, C_{ss min}/D levels were notably higher in patients using apixaban (p<0.0001 in both groups) (Figure 3). Furthermore, among patients using apixaban, the difference in C_{ss min}/D between subgroups with and without bleeding was found to be significant at p=0.9896, while in patients using rivaroxaban, the significance level was p=0.9480.

The study revealed that patients with bleeding disorders had a higher BMI compared to those without bleeding disorders (F(1,187)=4.6283, p=0.0327), with no significant differences based on drug intake. Furthermore, patients using apixaban exhibited a more advanced stage of chronic renal disease (F(1,192)=9.5407, p=0.0023), although there were no variations observed among groups with and without bleeding. Post-hoc analysis indicated a potential statistical significance within the bleeding group (p=0.0809). Factorial analysis of variance demonstrated a significant impact of the drug factor on HAS-BLED scores (F(1,192)=26.6144, p=6.164 x 10^-7). Specifically, it was found that bleeding patients taking apixaban had higher HAS-BLED scores compared to those taking rivaroxaban (p<0.0001), with similar differences noted in the non-bleeding group (p=0.0433). Notably, age, blood hemoglobin concentration, and platelet count did not exhibit variations across the groups and subgroups in the study.

3.2. The EWAS Did Not Reveal any Statistically Significant Variances

The study participants underwent screening for hereditary diseases associated with abnormal bleeding, such as hematuria. The analysis identified 14 carriers of pathogenic variants in patients with bleeding disorders. These variants were located in autosomes and affected the following genes: F5 (rs6025-T), F11 (rs121965064-C), VWF (rs41276738-T, rs900907976-A), ATP7B (rs76151636-T, rs28942076-T, rs780558532-TG, rs780558532-T), F7 (rs36209567-T), ABCC6 (rs745900279-delAG, rs72653706-A, rs1423674851-C), FANCA (rs397507553-delAGA), GP6 (rs368858591-A, rs199588110-A), FYB1 (rs769310295-insA), TNXB (rs200718357-T), SIK3(NM_001366686.3):c.741+1G>A and VPS33B(NM_018668.5):c.1105+2T>C. Consequently, based on the full-exome sequencing data, none of the participants were deemed to be at risk of hemorrhagic events with a sterile clinical presentation or late onset.

Based on the results of the EWAS conducted on the samples from the two groups, the average sequencing coverage depth of the target region was determined to be 86.88×, with an average coverage completeness of 97.15% at a depth of 10×.

Figure 4 illustrates the population affiliations of the samples, where it was not possible to definitively identify two or more clusters.

The variant rs7299963-C demonstrated the highest level of statistical significance in the recessive model (p=2.282 x 10^-4, C allele frequency 0.3454, T allele frequency 0.5707; OR 0.3968 [95% CI, 0.2638-0.597]). In the additive and dominant models, rs539002185-G yielded significant results (p=9.345 x 10^-6 in both instances, G allele frequency 0.1237, A allele frequency 0.2828; OR 0.358 [95% CI, 0.2112-0.6067]). The inclusion of covariates such as the first six principal components and medication did not modify the statistical significance. Both variants were located in intronic regions and exhibited a protective effect. The genomic factor prior to adjusting for different inheritance patterns was 1.03156.

The variant most strongly associated with bleeding in all three inheritance models, with drug covariate, was rs4298115-T (p=0.0007236 in the additive model, p=0.06448 in the recessive model, and p=0.0001805 in the dominant model, before correction; T allele frequency 0.5876, C allele frequency 0.404: OR 2.102 [95% CI, 1.405-3.145]). But Bonferroni correction increased the p value above 0.05. The protective variant with the lowest P-value before correction for multiple testing was rs343122-C (C allele frequency, 0.1443, A allele frequency, 0.3485; OR 0.3153 [95% CI, 0.1921-0.5177]).

Filtering the .vcf file based on regions containing the specified genes and rerunning the EWAS with the previously described algorithm did not reveal statistically significant differences, even with a reduction in the burden of multiple testing.

In order to establish the PRS for predicting bleeding, further quality control measures were implemented, resulting in the removal of duplicates and leaving a total of 11,938 variants.

We developed PRS focusing only on the presence or absence of bleeding (Figure 5). In one male participant, the specific type of bleeding remained undetermined. The determined PRS had a strong correlation with the investigated trait, explaining 98.2943% of its variance. This robust link was further supported by the PRS (61.7823) regression coefficient, which showed a highly significant statistical value (p<1 x 10^-174). The estimation’s precision was demonstrated by the small standard error (0.5843). An intriguing observation was the absence of overlap in PRS values between the two groups.

We assume that such unusual results might have been due to the small sample sizes and the minor effects of each variant. Additionally, we used a non-standard method for determining PRS, based on the logic of additive variant effects. For this reason, we decided to use PRS only as one of the independent variables, subsequently predicting the risk of bleeding by incorporating clinical parameters.

3.3. Frequencies of Variants in the Target Gene

Genotype frequencies of genes encoding enzymes metabolizing apixaban and rivaroxaban, including CYP3A4, CYP1A2, CYP2C8, CYP2C9, CYP2C19, CYP2J2, CYP3A5, as well as genes for the DOAC transporters ABCB1 and ABCG2 [33], were analyzed using plink (Table S2).

No initial quality control measures were conducted on the .vcf file in this instance. The sole multiallelic variant identified was the rs546527484 variant.

The rs2472304-G variant exhibited the most significant association with the risk of hemorrhagic events (prior to adjustment for multiple testing). In the additive model, the p-value was 0.008736, while in the dominant model, the p-value was 0.007099 before adjustment. The G allele frequency was 0.4794, the A allele frequency was 0.3434, with an odds ratio of 1.76 [95% CI, 1.172-2.644].

No significant differences were found among the prevalence of CYP2B6, CYP2C19, CYP2C8, CYP2C9, CYP2D6, CYP3A4, CYP3A5, CYP4F2, DPYD, NUDT15, SLCO1B1, TPMT and UGT1A1 haplotypes between the two sample groups (Table A1).

The study examined the diplotypes of the CYP3A4 gene, responsible for metabolizing rivaroxaban and apixaban, as presented in Table 3. Conjugacy tables were constructed, and statistical analysis was conducted using Fisher’s exact test, except for CYP3A4*1 and CYP3A4*36 where a chi-square test with Yates correction was applied. The results indicated no significant difference in the distribution of diplotypes. Notably, a trend towards statistical significance was observed in patients with the CYP3A4*1, CYP3A4*36 diplotype (chi-squared 2.0471, df=1, p=0.1525). Conversely, a less pronounced difference was noted for the CYP3A4*36 haplotype compared to other CYP3A4 haplotypes across the four groups (chi-squared 0.6854, df=1, p=0.4077). Furthermore, a slightly higher incidence of bleeding was observed in patients receiving apixaban, while a opposite pattern was noted for those on rivaroxaban.

3.4. Linear Regression Models

We observed a statistically significant correlation between PRS and urinary 6-β-hydroxycortisol and cortisol concentrations in patients using both medications. However, this correlation was particularly notable for rivaroxaban (Figure 6). Additionally, the random forest algorithm indicated that the PRS for apixaban was the least influential factor in determining ln C_{ss min}/D + 1 compared to other predictors.

A linear regression model (the lm() function in R) was developed to analyze patients using apixaban, incorporating predictors such as PRS, urinary 6-β-hydroxycortisol concentration, the ratio of 6-β-hydroxycortisol to cortisol, and the presence of chronic kidney disease (refer to Figure S4). The model demonstrated homoscedasticity (BP=3.472, df=4, p=0.4822), meaning that the predictor values had no effect on the residuals’ variance, which stayed constant. However, as the low multiple R-square (0.08972) and F-statistic value (2.439, p=0.05187) show, the model’s explanatory power was restricted. Notably, the ratio of urine 6-β-hydroxycortisol to cortisol concentrations was the only variable to exhibit a significant coefficient (p=0.0216), with an increase of one unit in the ratio resulting in a 0.1038 decrease in ln C_{ss min}/D + 1. These results suggest a poor predictive capacity of PRS, implying that genetics was not a major factor in apixaban-using patients.

The development of a linear regression model for patients receiving rivaroxaban revealed significant heteroscedasticity of the residuals, as indicated by the results of the Breusch-Pagan test (BP=14.021, df=5, p=0.01548). Predictors such as age, PRS, CHD, urine 6-β-hydroxycortisol concentration, and its ratio to cortisol were all included in the model. The generalized least squares (GLS) method was used, applying different standard deviations for each age value to address the heteroscedasticity issue. Significant improvements were seen after applying GLS: all predictor coefficients attained significance (p<0.05), and there were no longer any indicators of heteroscedasticity in the residuals (p=0.8567). The model estimates’ accuracy increased as the residual standard error dropped below 0.0138. Notably, ln C_{ss min}/D + 1 increased by 2.835 when PRS increased by one unit (t=20.4307, p<0.0001). On the other hand, age had a nonsignificant positive effect; an increase of one year led to a rise of 0.00066 in ln C_{ss min}/D + 1 (t=3.2296, p=0.0018). The presence of coronary ischemia was associated with a 0.867 decrease in ln C_{ss min}/D + 1 (t=–458.8183, p<0.0001). Furthermore, a unit increase in urinary 6-β-hydroxycortisol concentration resulted in a 0.00384 decrease in ln C_{ss min}/D + 1 (t=-57.8738, p<0.0001), while its ratio to cortisol was linked to a 0.157 increase in ln C_{ss min}/D + 1 per unit (t=49.0953, p<0.0001). By the way, these findings could be explained.

4. Discussion

The existing scientific literature contains a considerable number of studies that focus on investigating the associations between the presence of specific single nucleotide variants (SNV) genes related to transporter proteins and enzymes involved in the primary metabolic pathways of apixaban (ABCB1, CYP3A4) and rivaroxaban (ABCB1, CYP3A4, CYP3A5). The authors emphasize the impact of SNV markers of CYP3A4 (rs35599367), CYP3A5 (rs776746), and ABCB1 (rs1128503, rs2032582, rs1045642, rs4148738) genes on the pharmacokinetic parameters of POAK [12,13,34,35,36,37,38,39]. The relationship between genetic variations influencing DOAC pharmacokinetics and the risk of bleeding remains unclear. Research has indicated a higher incidence of bleeding in patients with AF who have elevated peak drug concentrations during DOAC administration, reaching 156% and 266% of the baseline for apixaban and rivaroxaban, respectively [40]. Pharmacogenetic investigations have suggested that specific genotypes and allele variants of ABCB1 (rs4148732) and ABCG2 (rs2231142) may lead to increased peak plasma concentrations of apixaban and rivaroxaban, although no direct associations with elevated bleeding risk have been established [13,40,41].

The search results from the PharmGKB database [42] for apixaban and rivaroxaban reveal the presence of pharmacogenetic markers with a level of evidence 3 (low). This level of evidence does not support the consideration of these markers for the integration of pharmacogenetic testing into clinical practice. The limitation is attributed to the lack of reproducibility of results across different studies [41]. Currently, there are no recommendations with a high level of evidence concerning the utilization of SNV of metabolizing enzymes (such as CYP3A4, CYP3A5, etc.) and transporters genes (such as ABCB1, ABCG2, etc.) for optimizing the therapy of patients with indications for prevention of atrial fibrillation-related stroke. Findings from association studies suggest that the carriage of pharmacogenetic markers does not significantly impact the exposure to DOACs, which would otherwise manifest in an increased risk of hemorrhagic complications.

Several reasons for the absence of these associations can be contemplated. Initially, while the examined variations impact the functional activity of apixaban and rivaroxaban biotransformation catalysts, additional enzymes (such as CYP1A2, CYP2C8, CYP2C9, CYP2C19, and CYP2J2) might also contribute to their metabolism [39], thereby offsetting the potential impact of SNVs in the primary pathway enzymes and maintaining the overall drug metabolism.

Secondly, within our context, the occurrence of bleeding as a clinical outcome cannot solely be attributed to genetically determined variations in the functional activity of DOAC metabolizing enzymes. The risk of hemorrhagic complications in individuals undergoing anticoagulant treatment is impacted by various clinical and demographic factors, with genetics representing just one aspect. Conversely, blood coagulation involves a intricate network of interactions among cascades of factors and cofactors, making it challenging to pinpoint the singular most critical component within this system.

Finally, studies using the candidate gene approach are guided by specific hypotheses and offer distinct advantages in terms of the robustness of statistical inference for identifying associations between genes and drugs. In the context of apixaban and rivaroxaban metabolism, the candidate gene approach proposes that variations in the ABCB1 P-glycoprotein gene or the CYP enzymes responsible for DOAC metabolism (CYP3A4/5, CYP1A2, CYP2C8, CYP2C9, CYP2C19, and CYP2J2) may be linked to drug concentration levels in the bloodstream and, consequently, to variations in clinical outcomes. Nonetheless, such studies may overlook numerous other genes and millions of SNVs that could potentially impact the variability in the pharmacokinetics and pharmacodynamics of these medications. Another limitation of candidate gene approach studies is their inability to fully consider the combined effects of multiple genetic variants simultaneously, known as polygenic effects. Research indicates that the variability in drug response among individuals may be attributed to polygenic inheritance [43,44].

To date, there is a limited number of GWAS and EWAS conducted to investigate genetic risk factors linked to drug response or pharmacokinetics of DOACs in patients [13,45].

In our investigation, an EWAS analysis was conducted to identify novel markers of bleeding risk. The results of our research revealed that none of the genetic variations in our sample exhibited a statistically significant correlation with the risk of hemorrhagic episodes at the exome level in both the apixaban and rivaroxaban treatment groups. Notably, the most noteworthy genetic variations in both groups were located within the CIT and ADPRHL1 genes, which have not been previously linked to coagulation disorders. Upon adding the medication as a covariate, robust correlations were observed in the ATP10D and UTP15 genes, which also have not been mentioned in the existing literature in the context of bleeding risk. In a separate analysis of SNVs of candidate enzyme metabolizer genes of the studied DOACs, only the rs2472304-G variant of the cytochrome CYP1A2 gene exhibited a significant association with bleeding risk under the additive and dominant inheritance models (prior to adjustment for multiple comparisons). Furthermore, a tendency towards a notable correlation with hemorrhagic incidents was identified in patients with CYP3A4*1 and CYP3A4*36 diplotypes who were prescribed apixaban.

The data obtained suggest that clinical and demographic parameters of patients, such as BMI, age, concomitant CKD and CHD, were predictors of C_{ss min}/D values. Patients taking apixaban, regardless of bleeding status, exhibited higher average C_{ss min}/D values (p<0.0001). However, the association between C_{ss min}/D values and bleeding risk was observed only in the rivaroxaban group.

An additional discovery indicated that the urinary levels of cortisol, its metabolite, and their ratio played a significant role in predicting C_{ss min}/D values for both apixaban and rivaroxaban. Cortisol is a specific substrate of CYP3A4, the metabolizing enzyme for these drugs. The determination of CYP3A4 activity can be achieved by calculating the metabolic ratio between cortisol and its metabolite 6-beta-hydroxycortisol: higher ratios reflect increased isoenzyme activity, while lower ratios suggest reduced activity. However, the reliability of this approach for CYP3A4 phenotyping is a subject of debate among researchers, with some questioning its practical utility in clinical decision-making [46,47,48].

PRS demonstrated an inverse correlation with the concentrations of 6-β-hydroxycortisol and cortisol in urine, showing a more significant impact in patients undergoing rivaroxaban treatment. However, the ratio of 6-β-hydroxycortisol to cortisol had varying effects on the ln C_{ss min}/D + 1 for the two medications: an inverse relationship was observed for apixaban, while a direct relationship was noted for rivaroxaban. Additionally, in prognostic models, it was found that in patients treated with rivaroxaban, the ln C_{ss min}/D + 1 index increased with PRS, age, the ratio of 6-β-hydroxycortisol to cortisol concentration in urine, and decreased with higher levels of 6-β-hydroxycortisol in urine, as well as in the presence of CHD as concomitant disease. However, detailed explanations for these findings cannot be provided.

5. Conclusions

Thus, no new polymorphic variants associated with an increased risk of bleeding following the use of apixaban and rivaroxaban were identified in the recent study using EWAS analysis. The research findings indicate that the residual equilibrium concentration of anticoagulants, including dose-adjusted, cannot be considered as a predictor of the risk of hemorrhagic complications in patients with non-valvular atrial fibrillation.

6. Limitation of Our Study

This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.

The study is subject to several limitations. The samples used in the study consisted of patients with AF who were receiving DOAC and were treated in various institutions within a single city, Moscow. These genetic samples were obtained as part of distinct observational pharmacogenetic studies focusing on apixaban and rivaroxaban. The observational design of the studies conducted across different clinics may have led to variations in routine clinical practices, potentially introducing biases in the collected data.

Patient data were extracted from medical records, which limiting the scope of the study to specific datasets (real-world data). The data collection and analysis focused mainly on standard clinical parameters typically obtained from patients diagnosed with atrial fibrillation. Nevertheless, it is important to acknowledge that there were discrepancies in clinical protocols among various healthcare institutions.

The clinical sample size was constrained by limited resources for patient recruitment and genetic analysis. Future studies with larger samples and improved control over inclusion criteria homogeneity may provide further insights into the relationship between DOAC and genetic factors.