Genetic Profiles in Pharmacogenes Indicate Personalized Drug Therapy for COVID-19

Data collection

Drugs and their clinical trial information for COVID-19 treatment were col lected from Guideline of Diagnosis and Treatment of Pneumonitis Caused by C OVID-19 (version 7.0) (http://www.nhc.gov.cn/yzygj/s7653p/202003/46c9294a7dfe4cef80dc7f5912eb1989.shtml), Clinical Trails (https://clinicaltrials.gov) and Chine se Clinical Trail Registry (http://www.chictr.org.cn). The pharmacogenes related to all drugs were obtained from DrugBank (https://www.drugbank.ca/) and Phar mGKB (https://www.pharmgkb.org/) databases. Actionable PGx biomarkers were determined by Clinical Pharmacogenetics Implementation Consortium (CPIC: http://cpicpgx.org/guidelines/) and DPWG (http://clivar-dpwg.iri.columbia.edu/) guid elines. Genetic variation data on each of the genes were retrieved from Gnom AD database (http://gnomad.broadinstitute.org/, version: 2.1.1). Finally, single-cel l RNA sequencing (sc-RNA seq) data from two organs (liver (GSE115469, n=5) and lung (GSE130148, n=4)) of 9 healthy adults was download from GEO da tabase (https://www.ncbi.nlm.nih.gov/gds).

Sc-RNA seq data analysis

The raw count matrix (unique molecular identifier counts per gene per cell) was processed by Seurat v3 ^[12]. Firstly, low-quality data was filtered out as following: (1) cells expressed less than 10 genes; (2) genes expressed in less than two percent of cells. Secondly, we calculated log₂(count-per-million + 1) expression, and followed by normalization and standardization. Thirdly, the cells were clustered. In detail, the highly variable genes were identified and cells were then clustered through embedding them into a graph structure in Principal Component Analysis space. The parameter resolution was set as 0.5 to identify only major cell types. The clustered cells were then projected onto a two-dimensional space and visualized. Finally, the differentially expressed genes on each cluster were identified. They were used to annotate and merge the cell clusters according to curated known cell markers (http://biocc.hrbmu.edu.cn/CellMarker/index.jsp).

Genetic variations annotation and drug-gene network construction

All genetic variations were annotated by allele frequency, location and function in different populations using ANNOVAR (version: 2019Oct24) (Table S1). In the current study, populations were divided into eight categories: African, Latino, East Asian, South Asian, Finnish, Non-Finnish European, Ashkenazi Jewish and Other. Based on the location and functions, all mutations were divided into 12 categories: frameshift, non-frameshift, synonymous single nucleotide polymorphism (SNP), nonsynonymous SNP, stop-gain, untranslated regions, intronic and others (including splicing, upstream of a gene, downstream of a gene, stop-loss and ncRNA_exonic mutations). Functional nonsynonymous mutations were predicted by PROVEAN (http://provean.jcvi.org/). The drug-gene network was constructed by Cytoscape software (version:3.7.1) ^[10].

Candidate drugs and their clinical trials for COVID-19 treatment

In the current study, a total of 67 drugs were collected in our study, including agents being used clinically, undergoing clinical trials, only confirmed in vitro experiments or simply with the therapeutic potential. They could be classified into nine categories according to their functions: viral toxicity (38.8%) and their sensitizers (3.0%), inhibiting virus invasion (6.0%), improving lung function (4.5%), maintaining tissue and organ function balance (25.4%), anti-hypoxic (4.5%), anti-septic shock (7.5%), anti-secondary infection (6.0%) and others (4.5%) (Figure 1A). Among them, antiviral drugs were mostly important and account for 47.8% in all drugs (including virus invasion inhibitors, viral toxicity and their sensitizers). As of March 18, 2020, a total of 125 clinical trials have been conducted to study the safety and efficacy of antiviral drugs. And most of these clinical trials are all ongoing and have not been completed. There are six clinical trials involving remdesivir, the earliest closed one involving 308 participants will be completed on April 27, 2020 (NCT04252664). As many as 16 clinical trials related to chloroquine, and two clinical trials for mild and common patients and severe patients will be completed on April 30, 2020 (ChiCTR2000029898, ChiCTR2000029899), each clinical trial includes a total of 100 participants. As for hydroxychloroquine, two of the nine clinical trials have ended on February 29, 2020 (ChiCTR2000029559, ChiCTR2000029740), and the results have not yet been available. In addition, the clinical trials of lopinavir (18) and ritonavir (23) will both be completed as early as March 25, 2020. Information of other categories of drugs can be found in Figure S1.

• Download figure
• Open in new tab

Figure 1: Candidate drugs for COVID-19 therapy.

(a) All potential drugs which may be utilized for COVID-19 therapy were divided into three categories and nine sub-categories as indicated. (b) All clinical trials conducted for antivirus drugs were listed in the figure. The number in each column indicated the number of trials registered now (until Mar 3, 2020). The earliest initial date and the earliest endpoint which was expected by trial conductors of trials for a same drug were also indicted here.

Drug-gene network for COVID-19 treatment

Genetic variation is one of the most powerful biomarkers to guide personalized therapy. Thus, it is important to systematically identify pharmacogenes which could affect response and toxicity of these 67 drugs. A drug-gene network which demonstrating the connection of each drug and pharmacogene was constructed (Figure 2A). After excluding drugs without pharmacogene data, a total of 44 drugs and 313 genes were connected.

• Download figure
• Open in new tab

Figure 2: Drug-gene network in COVID-19 treatment.

(a) The drug-gene network links 44 drugs and 313 candidate pharmacogenes as reported. All drugs were marked as diamond in the outermost, while all genes were marked as circles in the inner space. The colors of diamonds, circles and lines indicated the categories of related drugs in the drug-gene interactions. The kind of color of a circle reflects the number of drugs categories it links to. The size of the gene circle is determined by the number of linked drugs, and the transparency of the line it linked is proportional to the size of circles. (b) The DDIs predicted by drug-gene network. Solid lines indicate validated DDIs. The size of the diamond indicated the number of linked drugs. Drugs may interact with multi drugs should be used more carefully than other drugs.

Based on the nodes in this network, we could easily find genes shared by multiple drugs. They could be the mediator of drug-drug interactions (DDIs) (Figure 2B). In addition, their genetic variations could potentially affect response and toxicity of a number of different drugs. Thus, these genes should be paid attention during COVID-19 patient treatment. Among them, the top 2 genes both ligated 16 different drugs simultaneously. cytochrome P450 family (CYP) 3A4 is one of the major members of the cytochrome P450 enzymes superfamily, it involved in the metabolism of about half DCTs. The other one is ATP binding cassette subfamily B member 1 (ABCB1), which is an important transporter pumping out many drugs across the cellular membrane. Other top 10 genes included: solute carrier family solute carrier organic anion transporter family member 2B1, ALB, CYP3A5, CYP2C9, solute carrier family (SLC) 22A6, SLC22A1, CYP2C19 and CYP1A2. Most of them were related to absorption, distribution, metabolism and excretion. It is interesting to note that almost all of these gene are highly genetically variable and differ between populations

Another information obtained from this network was that there were actionable PGx biomarkers for three drugs based on CPIC or DPWG guidelines. They were G6PD and chloroquine hemolysis toxicity, vitamin D receptor (VDR) and ribavirin efficacy, inosine triphosphatase (ITPA) and ribavirin/α-interferon anemia risk, angiotensin I converting enzyme (ACE) and captopril response. We recommend that these genes should be preemptively tested before drug treatment.

Mutation profiles of pharmacogenes

We next investigated the allelic frequency, location and function of genetic mutations for all pharmacogenes in different populations. A total number of 295,897 variations in these 313 genes were found in 125,748 subjects from GnomAD databases. Although the mutations distributed in all gene regions, the most common mutation types were intronic, non-synonymous and synonymous mutations (Figure 3A). It should be noteworthy that 32.8% of the mutations were non-synonymous, which is the major types of functional PGx variants. They could be the explanation of individual and ethnic difference for COVID-19 drug treatment. To learn the distribution of mutations in pharmacogenes in more detail. The fractions of different mutation types were provided for each gene. They were summarized based on the gene functions in Figure S2. The frequencies of these variations were indicated in Figure 3B, they showed similar distributions. However, there were remarkably differences for the total number of mutations.

• Download figure
• Open in new tab

Figure 3: Variations profiles of pharmacogenes of drugs for COVID-19 therapy.

(a)The distribution of all variations in pharmacogenes, which were classed according to their locations and functions. (b) The distribution of all variations in pharmacogenes, which were classed according to their frequency. (c) The proportion of variations with different frequency in eight different populations. (d) The number and distributions of variations in four different categories of important pharmacogenes.

Then, we explored the allelic frequencies of all mutations. As indicated in Figure 3B, 98.52% of the mutations were rare genetic variants with minor allele frequency (MAF) lower than 1%. This result emphasized again the importance of pharmacogenomic scale test before treatment. We further analyzed mutations with different MAF stratified by populations. Except African population, there was no significant difference of MAF distribution between the remained seven races (Figure 3C). The African individuals were higher genetically variable for these pharmacogenes, since they showed more percentages of common mutations (MAF>1%). This result indicated that they had higher requirement for PGx test.

Taken together, these results indicated that COVID-19 pharmacogenes were highly genetically variable, especially for African population. Thus, personalized drug treatment and PGx test is needed for different patients. In addition, these mutations distributed in all gene regions and most of them were rare variations (Figure 3D). Therefore, sequencing technologies is recommended for the preemptive test.

Actionable and key PGx biomarkers

To provide clinical guide for COVID-19 treatment, we analyzed the actionable and key PGx biomarkers in detail. There are four paired actionable pharmacogenes and drugs: VDR and ribavirin, ITPA and ribavirin/α-interferon, G6PD and chloroquine, ACE and captopril.

Ribavirin is a classic antivirus drug which has been involved in the guideline for 2019-nCoV therapy as well. VDR polymorphism is associated with its efficacy (Figure 4A). There is only a common non-synonymous polymorphism M1T (rs2228570) that could hamper the VDR activity and the M1T carriers showed lower efficacy. The frequency of this mutation is over 0.5 in all populations (Figure 4A), which means more than half patients will be resistant to this drug.

• Download figure
• Open in new tab

Figure 4: Common nonsynonymous variations in important pharmacogenes.

Common nonsynonymous variations were indicated as amino alteration in the protein. The role of pharmacogenes in the treatment was indicated in the pathway. The frequencies of common nonsynonymous variations in different populations were also presented in the dot plot. Drug and pharmacogenes including: (a) VDR and ribavirin, (b)ITPA and ribavirin/α-interferon, (c) G6PD and chloroquine, (d) CYP3A4 and lopinavir, (e) ABCB1 and lopinavir, (f) ACE and captopril were presented here.

α-interferon is usually utilized in combination with antivirus drugs such as ribavirin in COVID-19 treatment. ITPA polymorphism was reported to be associated with the anemia risk in ribavirin/α-interferon treatment (Figure 4B). The functional mutation P32T (rs1127354) reduced the activity of the enzyme encoded by ITPA. The P32T carriers showed lower anemia risk. Another common non-synonymous polymorphism E93G (rs34982958) is also found in ITPA (Figure 4B), while its effect has not been revealed. The P32T was more identified in Asians (both of East Asians and South Asians, Figure 4B). Therefore, Asians may have a higher anemia risk in ribavirin/α-interferon treatment.

chloroquine has been utilized for 2019-nCoV therapy as guideline referred. G6PD polymorphisms has been reported to be associated with hemolysis risk (Figure 4C). There are two important variants V98M (rs1050828) and N156D (rs1050828), which could decrease the function of G6PD. The V98M and N156D carriers showed higher hemolysis risk for chloroquine treatment. In addition to these two mutations, we further identified three other common non-synonymous polymorphisms in all populations. Among them, the mutation S218F (rs5030868) was also reported to decrease the activity of G6PD, which indicated that this mutation could increase the hemolysis risk as well. What’s more, both V98M and N156D were more common in the African population (Figure 4C), suggesting that chloroquine is not recommended in African patients.

Lopinavir and ritonavir are used together to reduce their metabolism in liver through inhibiting the activity of CYP3A4, a very important enzyme which could metabolize about half of DCTs. The CYP3A4 polymorphism L292P (rs28371759, CYP*18B) was linked with increased CYP3A4 activity, thus facilitate the metabolism of drugs like lopinavir. L292P carriers may need to increase the dosage of these drugs compared with wild type patients (Figure 4D). Another common non-synonymous polymorphism R162Q in CYP3A4 was still not be reported with CYP3A4 phenotype now (Figure 4D). L292P was more in East Asians (Figure 4D), indicated that East Asians may metabolize lopinavir and ritonavir more rapidly.

As we mentioned above, ABCB1 can pumping out many drugs across the cellular membrane thus alleviate the efficacy. In our drug list, drugs like lopinavir can be transported out of cells by ABCB1, and its efficacy can be reduced (Figure 4E). We found there are 10 common non-synonymous polymorphisms in ABCB1 (Figure 4E). Among them, S893T, S893A (rs2032582), N21D (rs9282564) and S400N (rs2229109) were reported as deleterious mutations, which can increase the drug concentration through decreasing the efflux of ABCB1. S893A carriers were found in all populations with high frequency, which is up to 90% in Africans. S893T carriers were almost East Asians, while N21D carriers were Europeans (Figure 4E). These patients may be more response to these drugs transported by ABCB1.

In addition, ACE inhibitors (ACEI) can be used to improve the lung function in COVID-19 treatment. It was widely reported that ACE polymorphisms could affect ACEI therapeutic response (Figure 4F). The most important variant was I/D (rs1799752) polymorphism, which was a low frequency indel mutation causing low enzyme expression. The D allele carriers showed better response for ACEI treatment. In addition to this mutation, we further identified 10 common non-synonymous polymorphisms in all populations (Figure 4F). Although there is no related report about these variations, G2V (rs558593002), A154T (rs13306087), Y244C (rs3730025), R324W (rs35141294) and R561W (rs4314) were predicted to be loss-of-function mutations by PROVEAN. Thus, we speculated that patients harboring these single nucleotide polymorphisms (SNPs) were poor responders of this drug. Based on their MAFs in different populations, rs1799752 were more identified in the East Asian population (Figure 4F). Thus, we concluded that East Asian patients could be more resistance to this drug.

In summary, African patients could be more hazardous to chloroquine, Asian patients have a higher risk of anemia in ribavirin/α-interferon treatment and are more resistant to captopril. Meanwhile, drugs like lopinavir and ritonavir which metabolized by CYP3A4 or transported by ABCB1 should be carefully utilized in East Asians.

COVID-19 treatment based on sc-RNA seq data

Recently, sc-RNA seq technology was rapidly developed. It is capable to specifically profile cell populations at the single-cell resolution. Thus, it transformed many fields of genomic research. Here, we tried to explore the utility of sc-RNA seq on personalized COVID-19 treatment. We collected healthy adults sc-RNA seq data from 4 lung and 5 liver tissues.

Lung is the major 2019-nCoV attacking organ. Thus, we analyzed the expression of virus or drug target genes in different lung cells (Figure 5A). Based on the current findings, ACE2 and transmembrane serine protease 2 (TMPRSS2) were the two major targets of both virus and some drugs. As indicated in Figure 5B, TMPRSS2 was almost expressed in all kinds of cells. In our enrolled drugs, camostat mesylate act as a trypsin like protease inhibitor and could be potentially used to anti-virus by inhibiting TMPRSS2 activity. Due to the ubiquitous expression of its target, camostat mesylate could also cause potential toxicity when exerting therapeutic effect. On the other hand, ACE2 was mainly expressed in type I pneumocyte, type II pneumocyte, FOXN4⁺, club and mast cells. chloroquine phosphate can change the ACE2 structure, or inhibit the binding of coronavirus S protein with ACE2. Based on these data, chloroquine was more effectively targeting drugs with less side effects (Figure 5C). Sc-RNA seq data in lung tissue suggested that ACE2 targeting drugs (chloroquine) were superior for COVID-19 treatment than TMPRSS2 targeting drugs (camostat mesylate).

• Download figure
• Open in new tab

Figure 5: Discrepancy of pharmacogenes’ expression in target genes.

(a) Lung cells were clustered by UMAP method. Different cell populations were marked in different colors. (b) Violin plot of ACE2 and TMPRSS2 expression in different cell types of lungs. (c) Mechanism of chloroquine and Camostat mesylate acting on different lung cells by identifying targets.

Most drugs can be metabolized in the liver. We next analyzed the expression of drug metabolic enzymes in different liver cells (Figure 6A). CYP1A2, CYP2E1, CYP2C8, CYP2C9, CYP2B6, CYP2D6, CYP3A4 and CYP3A5 were the major CYP450 enzymes responsible for 67 drugs metabolism. Their expression was investigated (Figure 6B, 6D). In addition to hepatocytes, it was unexpected to found that some enzymes were also found to be highly expressed in immune cells. Based on the expression distribution, we categorized all enzymes into two groups: (1) narrow distribution enzymes (NDE) which were only expressed in hepatocytes and scarcely expressed in other cells, including CYP3A4, CYP1A2, CYP2B6, CYP4F2 and CYP4F12; (2) wide distribution enzymes (WDE) which is expressed in both hepatocytes and immune cells in the liver, including CYP2C8, CYP2C9, CYP2D6 and CYP3A5 (Table 1). For our enrolled drugs, anti-viral agents such as α-interferon, ambroxol, baloxavir and darunavir were mainly metabolize by NDEs, while thalidomide and valsartan were mainly metabolized by WDEs, and some drugs such as chloroquine and leflunomide were metabolize by both NDEs and WDEs. Currently, the consequence of drug metabolism in immune cells still remains unknown. However, it was rational to speculate that the toxically drug metabolites in the immune cells would affect liver immunologic function. Thus, in case of producing toxically metabolites, the drugs that can be metabolized by WDE would be more tending to cause hepatotoxicity. They should be cautiously used when alternative drugs were available. Leflunomide can be metabolized by WDEs CYP2C9 which expressed both in hepatocyte and CD4 ⁺ T cells. Its active metabolite, A771726, inhibits cell survival and proliferation by inhibiting pyrimidine synthesis. The aggregation of A771726 in CD4 ⁺ T cells inhibits cell survival and proliferation, and damages the liver immunologic function (Figure 6C). Therefore, leflunomide may be not recommended while other immunosuppressive drugs are available.

• Download figure
• Open in new tab

Figure 6: Discrepancy of pharmacogenes’ expression in CYP450 enzymes.

(a) Liver cells were clustered by UMAP method. (b) Violin plot of CYP2C9 expression in different cell types of livers. (c) The mechanism of Leflunomide inhibiting cell proliferation after metabolism by CYP2C9. (d) Violin plot of CYP2D6, CYP3A4, CYP2B6, CYP2C8, CYP1A2, CYP2E1, CYP3A5, CYP4F2 and CYP4F12 expression in different cell types of livers.

Candidate Strategies for COVID-19 Therapy

Since there are many candidate drugs for COVID-19 therapy, it is important to make a better choice for each patient. For example, ribavirin and α-interferon are not suggested in patients with M1T variation in VDR or without P32T variation in ITPA, then other anti-virus drugs can be used in these patients. Similarly, chloroquine is not suggested to be utilized in patients with V98M and N156D variants. To avoid these abuses of drugs in COVID-19 therapy, a genetic detection panel is suggested. In this panel, SNPs with clinic evidences in pharmacogenes should be considered. The mutation frequencies are not the same in all populations, and we provided the frequency of candidate pharmacogenes SNPs in all populations (Figure 7). Mutations like rs117648444 is only found in FIN, so that the panels for FIN should include this variation. Meanwhile, suggestions should be given based on these SNPs. Examples can be found in Table 2.

• Download figure
• Open in new tab

Figure 7: The frequency of important SNPs in eight populations.

Drugs were listed in the top of the panel, while the SNPs were listed in the bottom of the panel. The color of the column is as same as the color of the drug, which indicated the category of the drug.

As a large part of drugs are metabolized by CYP450 in liver, the activity of CYP450 family may impact many drugs metabolization rate. Due to the activity difference exists among populations, the genotypes of important alleles in CYP450 family should be detected after diagnosis, then treatment strategies for these drugs can be suggested in these patients. Meanwhile, DDI should also be considered as many drugs may be utilized together during COVID-19 treatment as mentioned above. For instance, drugs like lopinavir and ritonavir can inhibit the activity of CYP3A4, thus the drugs dosage should be adjusted when combined with drugs metabolized by CYP3A4. Drugs that were metabolized by CYP450 family in COVID-19 therapy were listed in Table 3. In addition, when chloroquine is utilized, drugs like losartan which can increase the concentration of chloroquine may be replaced by valsartan. Validated DDIs between drugs for COVID-19 treatment are given in table 3.

Finally, based on current researches, 2019-nCoV can only infect cells with ACE2. As indicated above, ACE2 is not expressed in all cells, while TMPRSS2 can be found in large part of cells. So that drugs like camostat mesylate is not recommended compared to ACE2 specific inhibitors due to the side effect caused by TMPRSS2 in uninfected cells.