extfieldsofstudy: Chemistry Medicine

text:
Protein markers of dysfunctional HDL in scavenger receptor class B type I deficient mice

Background Scavenger receptor class B type I (SR-BI) plays a key role in high density lipoproteins (HDL) metabolism. SR-BI deficiency in mice results in enhanced susceptibility to atherosclerosis with abnormal large, cholesterol enriched, and functional impaired HDL. This study was to characterize the protein markers of dysfunctional HDL in SR-BI deficient (SR-BI−/−) mice and to test if the defective of HDL might be affected by probucol treatment. Methods Shotgun proteomics and 2-D gel electrophoresis were performed to examine the profile of HDL protein and distribution of HDL particles isolated from SR-BI−/− mice. HDL’s cell-function, paraoxonase 1 (PON1) and myeloperoxidase activity were assessed. The mice were treated with 1.2 mg/g/day probucol for 6 weeks and the impact on HDL protein markers was analyzed. The differential proteins were quantified by Western blotting. Results The relative amount of protein in SR-BI−/− HDL was decreased by about 25% compared to that in HDL from wild type (WT) mice. Compared to WT HDL, relative protein abundance of representative apoAI and PON1 in SR-BI−/− HDL were significantly reduced, whereas acute-phase protein serum amyloid A (SAA) and apoAIV, proteinase inhibitor proteins α-1-antitrypsin (A1AT) were increased. The distribution of plasma apoAI-containing HDL particles in SR-BI−/− mice was also dramatically altered, although plasma apoAI level was no difference. The protein alterations were accompanied with dysfunction of SR-BI−/− HDL, evidenced by impaired cholesterol homeostasis in macrophages, and reduced anti-oxidative and anti-inflammatory effects. Probucol treatment of SR-BI−/− mice could restored the relative contents of critical proteins including apoAI, PON1, SAA, apoAIV and A1AT on HDL, and improve HDL dysfunction despite decreased HDL-C level. Conclusion SR-BI deficiency leading to dysfunctional HDL is closely related to alteration of HDL protein, suggesting that identification of apoAI, PON1, SAA, apoAIV, and A1AT may serve as the valuable protein markers for diagnosis and therapeutics of dysfunctional HDL-related metabolic diseases. Electronic supplementary material The online version of this article (10.1186/s12967-018-1502-y) contains supplementary material, which is available to authorized users.

Background
Plasma high-density lipoprotein (HDL) is a well-established cardiovascular protective factor with multiple functions, including mediating reverse cholesterol transport, and exerting anti-oxidative, anti-inflammatory, anti-thrombotic and vasodilatory activities [1]. Epidemiological and clinical studies have shown that HDL cholesterol (HDL-C) level is inversely correlated with coronary heart disease risk [2]. However, nearly half of cardiovascular clinical events occurred in subjects with normal or even high levels of HDL-C [3]; and clinical trials with drugs that raised HDL-C levels by cholesteryl ester transfer protein (CETP) inhibitors, had failed to reduce the risk of cardiovascular events [4]. These controversial findings suggest that HDL-associated functional factors, instead the levels of HDL-C, are important for the antiatherogenic properties. Indeed, it has been known that HDL particles have substantial compositional heterogeneity in its lipids, apolipoproteins, lipid transfer proteins, enzymes and numerous other proteins under pathophysiological conditions because of their complex metabolic function. The evidence from transgenic mice indicates that protein composition determines the antiatherogenic properties of HDL. Alterations in protein composition involved in HDL metabolism could promote atherosclerosis, even if plasma HDL-C was normal or elevated [5]. Recent studies come up with a new concept of "dysfunctional HDL", which suggests that antiatherogenic HDL can be converted to proatherogenic HDL under dyslipidemic and inflammatory states in metabolic diseases, and thus may accelerate atherosclerosis [6]. Therefore, the assessment of dysfunctional HDL by measuring the protein composition, and the identification of biomarkers of dysfunctional HDL in different experimental models, will provide important information for the evaluation of cardiovascular risk in patients and for the development of new antiatherogenic therapies [7].
Scavenger receptor class B type I (SR-BI), a welldefined HDL receptor, is a prominent regulator of HDL metabolism, which controls the plasma HDL-C level. It mediates the selective uptake of cholesteryl esters from HDL by the liver and steroidogenic organs, and facilitates the efflux of cholesterol from peripheral tissues to HDL [8]. It has been reported that several mutations in human SR-BI gene, which are associated with low SR-BI protein levels, caused elevated plasma HDL-C levels, a reduction in cholesterol efflux from macrophages, altered platelet function, and decreased adrenal steroidogenesis [9][10][11]. The SR-BI deficient (SR-BI −/− ) mice showed a nearly twofold increase in plasma total cholesterol (TC) and HDL-C levels, accompanied with impaired anti-oxidative function, and lipid deposition in the aorta, contributing to the enhanced susceptibility to atherosclerosis [12][13][14]. However, the protein composition of defective HDL in SR-BI −/− mice has not been well studied.
Therefore, in this study we used shotgun proteomics, a powerful tandem mass spectrometry (MS/MS) approach following the protein separation by liquid chromatography (LC) [15], to compare the protein composition of HDL particles isolated from SR-BI −/− mice and wild type (WT) mice. Meanwhile, HDL-mediating cholesterol flux in macrophages, anti-oxidative and anti-inflammatory function of HDL were detected to verify the functional status of HDL. Considering that probucol as a cholesterol-lowering and antioxidant drug, could alter the abnormal HDL-C level in SR-BI −/− mice [16], probucol treatment was used to investigate whether the alteration of protein markers and HDL function could be affected in SR-BI −/− mice. The goal of this study is to quantify the protein characteristics of SR-BI −/− HDL, providing the valuable protein markers for diagnosis and therapeutics of dysfunctional HDL-related metabolic diseases.

Mice
C57BL/6 WT mice and apolipoprotein E deficient (apoE −/− ) mice were obtained from Vital River Laboratory Animal Technology Company, China. SR-BI heterozygous (SR-BI +/− , 1:1 mixed C57BL/6 × S129 genetic background) mice were kindly provided by Professor MacRae F. Linton (Vanderbilt University, USA). SR-BI −/− mice were generated from heterozygous mating pairs or by mating female SR-BI +/− with SR-BI −/− males. SR-BI genotypes were determined by PCR analysis of DNA extracted from tail biopsies [17]. All of the mice were housed in microisolator cages on a rodent chow diet. Animal care and experimental procedures were performed under the regulations of the Institutional Animal Care and the Ethics Committee for Animal Experiments of Wuhan University.

Isolation of plasma HDL and determination of HDL composition
Blood samples from WT (SR-BI +/+ ) or SR-BI −/− mice (both male and female) were collected with heparin collection tube by retro-orbital venous plexus puncture after overnight fasting. The plasma was immediately separated by centrifugation at 4000×g for 10 min at 4 °C. HDLs were isolated from pooled plasma samples by ultracentrifugation described in our previous work [18]. Briefly, apoB containing lipoproteins in plasma were precipitated by the addition of 20% polyethylene glycol (PEG) in 200 mM glycine (pH 10.0, plasma/PEG ratio, 1.0: 0.04), then the HDL-containing supernatant was dialyzed against 0.9% NaCl, 0.3 mM EDTA (pH 8.0), and the HDL was further isolated at density < 1.21 g/mL by ultracentrifugation, desalted and concentrated in phosphate-buffered saline (PBS) to a final concentration of 1 μg/µL. Total cholesterol (TC), free cholesterol (FC) and triglyceride (TG) levels in HDL were determined enzymatically using Cholesterol and Triglycerides GPO Reagent kits (Raichem, USA). Protein concentration of HDL was determined by Lowry DC kit (BioRad, USA). HDL phospholipid was measured using Wako Phospholipids B enzymatic kit (Wako Laboratory Chemicals, USA). HDL composition was analyzed by calculating the percentage of TC, FC, TG, phospholipid, and total proteins.

Shotgun LC-MS/MS
LC-MS/MS (liquid chromatography-tandem mass spectrometry) was performed at State Key Laboratory of Proteomics in the Beijing Proteome Research Center. Equal protein amounts (50 µg) of HDL preparations were separated by 10% SDS-PAGE (120 V, 20 min), and the totality of each lane was then sliced into 10 bands followed by in-gel digestion with trypsin. The resulting peptides were analyzed by reverse phase liquid chromatography coupled with tandem mass spectrometry (RP LC-MS/MS) on an LTQ-Orbitrap Velos mass spectrometer (Thermo Electron, San Jose, CA). All MS/MS spectra were identified by using Sorcerer-SEQUEST Algorithm version 4.0.3 (Sage-N Research, San Jose, CA) against the Refseq mouse protein sequence database. The searching parameters were set up as follows: partial trypsin (KR) cleavage with two missed cleavage was considered; the variable modification was oxidation of methionine; carbamidomethylation was set as fixed modification; the peptide mass tolerance was 20 ppm, and the fragment ion tolerance was 1 Da. Functional annotation of proteins detected in HDL was obtained from UniProt (http:// www.unipr ot.org/). Spectral counts were recorded to quantify relative protein abundance. This approach is based on the correlation between protein abundance and spectral counts. The higher abundance a particular protein is, the more MS/MS spectra are collected for peptides of that protein [19].

Western blotting analysis
To identify the levels of specific proteins in HDL or plasma of SR-BI −/− mice, we further used Western blotting to quantify the different proteins detected in the proteomic analysis. Appropriate amount of HDL proteins (2 µg HDL proteins for checking apoAI, and 15 µg HDL proteins for checking other proteins) or mouse plasma (0.5 µL plasma for checking apoAI and paraoxonase 1 (PON1), 1 µL plasma for checking other proteins) was loaded and separated by 10% SDS-PAGE, and transferred onto nitrocellulose membranes, then target proteins were detected by primary antibodies including goat anti-apoAI (Ab7614), mouse anti-PON1 (Ab24261), α-1-antitrypsin (A1AT) rabbit (Ab133642) antibodies from Abcam, goat anti-apoAIV (Santa, SC-19036), and rat anti-SAA antibodies (Santa, SC-59680), and HRP-conjugated secondary antibody. Signal was detected using an enhanced chemiluminescence (ECL) kit.

2D gel electrophoresis
ApoAI-containing HDL distributions of SR-BI +/+ and SR-BI −/− mice were analyzed by 2D nondenaturing agarose-polyacrylamide gel electrophoresis (2DE) modified from Asztalos et al. [18,20]. Briefly, in the first dimension, 2 µL of fresh plasma lipoproteins were separated on 0.75% agarose gel. Individual agarose strips were transferred to the top of 3-16% nondenaturing polyacrylamide gradient gels and separated in the second dimension. Lipoproteins were transferred to 0.2 µm nitrocellulose membranes, and 2DE patterns of apoAI-containing subpopulations were detected by rabbit anti-mouse apoAI primary and HRP-conjugated secondary antibody. Signal was detected using an ECL kit.

Measurement cholesterol content in macrophage
Mouse peritoneal macrophages from WT C57BL6 mice were harvested by peritoneal lavage 3 days after intraperitoneal injection of 2.5 mL of 3% thioglycollate. Macrophages were suspended in 5% FBS/DMEM and seeded onto 6-well plates at a density of 1 × 10 6 cells/well. After 3 h, non-adherent cells were removed by washing with DMEM, the macrophages were incubated for 48 h in DMEM alone or containing either SR-BI +/+ HDL or SR-BI −/− HDL (12.5, 25, 50 µg protein/mL), then washed twice with DPBS and air dried. Cell lipids were extracted by overnight incubation at room temperature in isopropanol. The FC and TC contents in the lipid extracts were then measured by Cholesterol Reagent (Raichem). Cell proteins were solubilized by adding 1 N NaOH to the wells, and the protein concentration was measured by Lowry DC protein kit (BioRad, USA).

PON1 activity and MPO activity assays
The measurement of PON1 activity in plasma or HDL particle using paraoxon (Sigma, USA) as substrate was described in our previous work [21]. One unit of PON1 activity in serum or HDL was defined as 1 nmol of p-nitrophenol formed per minute under the definition conditions and expressed as U in per ml serum or per mg HDL. Myeloperoxidase (MPO) activity was measured by a MPO determination kit (Jiancheng Biotech Ltd, China), according to the manufacturer's instructions.

Determination of HDL oxidative/anti-oxidative properties by cell-free assay
To determine the anti-oxidative properties of HDL, the cell-free assay was used with modifications [22,23], using apoE −/− mice plasma as the fluorescence-inducing agent. HDL was isolated by ultracentrifugation or from plasma after apoB-lipoprotein precipitation using dextran sulfate [23]. The change in fluorescence intensity as a result of the oxidation of dichlorodihydrofluorescein (DCFH) by apoE −/− mice plasma in the absence or presence of HDL was tested. Briefly, 0.3-0.5 µL plasma from apoE −/− mice and 15 µL of isolated HDL (final concentration of 15 µg/ mL) or 15 µL of HDL-containing dextran sulfate supernatant were incubated with PBS to a total volume of 200 µL per well in black polystyrene microplates. After 90 min of incubation on a rotator at 37 °C, 2 µL of 2.0 mg/mL fresh DCFH-DA in methanol solution was added to each well, and incubated for additional 30 min with rotation. Fluorescence signal was determined (excitation/emission wavelength = 485 nm/530 nm, cutoff of 515 nm) using a Fluorescent Plate Reader (GENIOS VA2000 plate reader, TECAN). Assay was performed in triplicates and values in the absence of HDL were normalized to 1. After the addition of the test HDL, values < 1.0 indicated anti-oxidative HDL, values > 1.0 indicated oxidative HDL.
Meanwhile, the cultured macrophages were collected using Trizol reagent, total RNA of macrophages was isolated and reverse transcribed into cDNA with reverse transcriptase (Invitrogen, USA). The target region in mouse MCP-1 and TNF-α gene was amplified by quantitative real-time PCR (Q-PCR) using universal PCR master mix (Invitrogen, USA). The primers for MCP-1 were 5′-AGG TCC CTG TCA TGC TTC TG-3′ and 5′-TCT GGA CCC ATT CCT TCT TG -3′. The primers for TNF-α were 5′-CGT CAG CCG ATT T GCT ATC T -3′ and 5′-CGG ACT CCG CAA AGT CTA AG -3′. The primers for 18s rRNA were 5′-CGC GGT TCT ATT TTG TTG GT-3′ and 5′-AGT CGG CAT CGT TTA TGG TC-3′. The relative expression levels of these genes were detected by normalization of the cycle threshold (Ct) to that of the control 18S RNA. 2 −∆∆Ct was used to represent the relative expression of target genes.

Effects of probucol treatment on SR-BI −/− mice
Eight male SR-BI −/− mice at 19-20 weeks of age (body weight: 22.7 ± 1.1 g) were randomly divided into control group and probucol treatment group. All mice were fed a normal chow diet. Probucol (1.2 mg/g/day) was intragastrically given to SR-BI −/− mice for 6 weeks and 0.5% CMC-Na as the control. The levels of plasma TC and TG were determined enzymatically. Lipoprotein cholesterol profiles were made by fast protein liquid chromatography (FPLC) as described our previous study [18]. The expression levels of specific plasma HDL-associated proteins in SR-BI −/− mice after probucol treatment were detected by Western blotting. Plasma PON1 activity, MPO activity and oxidative/anti-oxidative properties of HDL were measured using above methods.

Statistical analysis
Statistical analysis was performed using Prism 5. Continuous variables were expressed as mean ± SD (standard deviation). The two-tailed Student's t test for parametric variables was used to assess difference between two groups. P < 0.05 were considered statistically significant.

Lipid and protein contents in SR-BI −/− HDL
HDLs isolated from WT (SR-BI +/+ ) and SR-BI −/− mice by ultracentrifugation were analyzed for the composition of lipids and proteins. In lipid classes, SR-BI −/− HDL displayed a significantly increased percentage of FC and TC. FC content increased to 20.90% by 3.7-fold as compared with only 4.42% FC in SR-BI +/+ HDL, and TC content increased from 18.83% (SR-BI +/+ HDL) to 29.11% (SR-BI −/− HDL); while phospholipid content was not significantly altered. The relative amount of protein in SR-BI −/− HDL was decreased by about 25% compared with SR-BI +/+ HDL (Additional file 1: Table S1), suggesting that SR-BI deficiency in mice markedly affected HDL protein and lipids ratio.

Proteomic alterations of SR-BI −/− HDL
To identify HDL-associated proteins, we performed proteomic analysis of purified HDL using LC-MS/MS. After removing predicted contaminants such as keratin, 50 HDL-associated proteins in WT mice and 78 proteins in SR-BI −/− mice were identified as true-positives (proteins were required to be present at > 0.05% of the total spectral counts). To assess the alterations of the contents of the proteins involved in HDL function, identified HDL-associated proteins were grouped into functional categories. The result showed that some proteins involved in immune response like beta-2-microglobulin, some in cell adhesion such as insulin-like growth factor-binding protein, integrin αIIb, integrin β1, and some others in proteinase inhibition such as antithrombin-III, were only been identified in SR-BI −/− HDL, but not in SR-BI +/+ HDL. As the majority of HDL-associated proteins were common to both control SR-BI +/+ and SR-BI −/− HDL, relative protein abundance was quantified by spectral counts, and resulting outputs were expressed as a percentage of the total spectral counts obtained per group. The analysis revealed that among the 78 proteins identified in SR-BI −/− HDL, 26 proteins were appeared or increased and 10 proteins were decreased, when compared to those in SR-BI +/+ HDL. The relative levels of apoAI, apoAII, apoCI, apoCII, apoM and PON1 were significantly reduced, whereas apoE, apoH, LCAT, acutephase proteins apoAIV, SAA, complement C3, several proteinase inhibitor proteins, including A1AT, inter alpha-trypsin inhibitor and α-2-macroglobulin were increased in SR-BI −/− HDL, as compared to those in SR-BI +/+ HDL (Additional file 2: Table S2).
Gene ontology classification was applied to estimate the contribution of various functional protein classes present in HDL. Compared to SR-BI +/+ HDL, HDL from SR-BI −/− mice had a significant decrease in proteins involved in lipid metabolism (37.34% versus 57.98%) and antioxidant (1.94% versus 7.09%), while had a significant increase in inflammatory and immune response (22.14% versus 9.19%) and proteinase inhibition (15.67% versus 8.49%) (Fig. 1a).
To further validate the proteome results, Western blotting analysis of representative HDL proteins were performed and demonstrated that apoAI and PON1 levels were decreased, whereas apoAIV, A1AT and SAA levels were increased in SR-BI −/− HDL (Fig. 1b).

Differential profile of plasma HDL proteins
We further measured the levels of those proteins in mouse plasma, the data showed that the changes of PON1, SAA, A1AT and apoAIV levels in plasma were consistent with those in HDL, but levels of apoAI in plasma were only slightly reduced in SR-BI −/− mice compared to SR-BI +/+ mice (P > 0.05) (Fig. 2a).
Using 2D gel electrophoresis, we investigated the effects of SR-BI deficiency on apoAI-containing HDL subpopulations in the plasma (Fig. 2b). The profile showed that the majority of apoAI in plasma of SR-BI +/+ mice was in α-mobility and relatively small amount in pre-β migrating HDL. Compared to SR-BI +/+ mice, SR-BI −/− mice had an accumulation of large HDL migrating between pre-β and α, suggesting that the compositions and structure of apoAI-containing HDL in SR-BI −/− mice was dramatically altered.

Functional alterations of SR-BI −/− HDL HDL-mediating cholesterol flux in macrophage
Considering that modification of apoAI and changes of HDL composition may impair the capacity of HDL in maintaining macrophage cholesterol homeostasis, we measured the cellular cholesterol content of WT macrophages after incubated with isolated HDL at various protein concentrations (12.5, 25 or 50 mg/mL) for 48 h. Compared to the cells treated with SR-BI +/+ HDL, the cells treated with SR-BI −/− HDL had significantly more TC and FC remained in the cells in a concentration dependent manner (Fig. 3a).

Anti-oxidative activity of HDL
PON1, as a specific HDL-binding protein, contributes to the anti-oxidative activity of HDL. In agreement with proteomic data, PON1 activity was significantly reduced in both HDL and plasma of SR-BI −/− mice (Fig. 3b). The PON1 activity in SR-BI −/− HDL was decreased by 50% as compared with that in SR-BI +/+ HDL (55.7 ± 6.4 U/mg versus 111.9 ± 21.1 U/mg, P < 0.01), and plasma PON1 activity of SR-BI −/− mice was significantly lower than that of WT mice (108.2 ± 17.5 U/mL versus 189.3 ± 33.8 U/ mL, P < 0.01). Furthermore, we found that the plasma activity of MPO, an oxidative enzyme related to inflammation, was significantly higher in SR-BI −/− mice.
In addition, the result of the oxidative/anti-oxidative properties of HDL detected by the cell-free assay showed that isolated HDL by ultracentrifugation or plasma HDL after apoB-lipoprotein precipitation from SR-BI −/− mice had oxidative status with values > 1.0, compared with the anti-oxidative HDL from SR-BI +/+ controls, whose values were < 1.0 (Fig. 3c).

Anti-inflammatory capacity of HDL
Macrophages respond to ox-LDL stimulation by expressing and secreting inflammatory mediators. We tested the effects of SR-BI +/+ HDL and SR-BI −/− HDL on ox-LDL-induced cytokine production in macrophages. As shown in Fig. 4a, b, ox-LDL resulted in a strong induction of both MCP-1 and TNF-α expression and secretion in macrophages, and SR-BI +/+ HDL significantly diminished the induction at both mRNA and protein levels, while SR-BI −/− HDL showed no effects. This data suggests that SR-BI −/− HDL lost its anti-inflammatory property.

Remodeling structure and function of SR-BI −/− HDL by probucol treatment
After SR-BI −/− mice were treated with probucol, plasma TC and FC levels were significantly decreased, mainly in the HDL-C fractions, compared with the control group (Fig. 5a). However, the HDL-binding proteins apoAI and PON1 levels in plasma were significantly increased, and acute-phase proteins apoAIV, SAA and A1AT levels were reduced in SR-BI −/− mice treated with probucol (Fig. 5b); and the anti-oxidative function of HDL in SR-BI −/− mice could be improved, evidenced by increases in plasma PON1 activity and decreases in HDL oxidative property and plasma MPO activity (Fig. 5c).

Discussion
In current study, we demonstrate that the defective HDL functions in SR-BI −/− mice are associated with the altered protein composition by shotgun proteomic profiling and biochemical analyses, despite the higher HDL-C level in SR-BI −/− mice. We observed that a significantly lower abundance of lipid metabolism and antioxidant related proteins, such as apolipoproteins (apoAI, apoAII and apoCI) and HDL-associated PON1 in SR-BI −/− mice, whereas levels of proteins involved in inflammatory and immune response, proteinase inhibition and blood coagulation, including SAA, apoAIV, complement C3 and A1AT, were increased. These protein composition  16:155 alterations result in impairment of HDL function, including impaired function in maintaining cholesterol homeostasis in macrophages, and the anti-oxidative and anti-inflammatory capacity. Furthermore, the lipidlowering and antioxidant drug probucol can improve the HDL functions in SR-BI −/− mice by reversing the contents of some critical proteins, including increasing plasma apoAI and PON1 levels, and decreasing plasma SAA, apoAIV and A1AT levels.
There is accumulating evidence that alterations in HDL protein quality and composition result in HDL function impairment [24,25]. In this study, we employed shotgun proteomics methods to identify the quantitative and qualitative alterations in HDL-associated proteins that may render SR-BI −/− HDL dysfunctional. The proteomic profiling of HDL particles shown that some of essential protective apolipoproteins and enzymes were decreased, and the oxidative and inflammatory related proteins identified in SR-BI −/− HDL appeared or were increased, compared to SR-BI +/+ HDL. We observed that relative levels of HDL-associated apolipoproteins involved in lipid metabolism, such as apoAI, apoAII, and apoCI, were significantly decreased in SR-BI −/− HDL. As the key apolipoprotein on HDL, apoAI is a well-known anti-atherogenic molecule by mediating reverse cholesterol transport (RCT) and exerting anti-inflammatory effects [26][27][28]. ApoAII, as the second most abundant protein on HDL, also acts as a primary acceptor and efficiently removes cholesterol from macrophages in vivo [29]. ApoCI, as a component of HDL and VLDL, was considered to play a beneficial role in lipid metabolism, but recently study reported that apoCI could increase the LPS-induced inflammation in macrophages in vitro and in apoE −/− mice [30]. In the proteomic profiling, we found that the relative abundance of apoCI in SR-BI −/− HDL was decreased. Therefore, the role of  16:155 apoCI in SR-BI −/− mice remains to be further studied. An interesting finding is that the percentage of apoAI is decreased in HDL by spectral counting, while the plasma apoAI content is not reduced in SR-BI −/− mice. Meanwhile, we also identified by 2D gel electrophoresis that the SR-BI −/− mice had an accumulation of large HDL migrating between pre-β and α, with a higher level of apoAI. These suggested that apoAI in abnormal SR-BI −/− HDL might have underwent dramatic modification, and became loosely associated with the larger HDL, thus, part of apoAI was likely to lost from isolation of HDL by precipitation and ultracentrifugation. Considering that SR-BI −/− mice had comparable plasma apoAI concentration as WT mice (108 versus 100 mg/ dL), as reported by a previous study [12], it is conceivable that the relative level of apoAI in SR-BI −/− HDL Besides involved in lipid metabolism, apoAI also exerts the anti-oxidative and anti-inflammatory functions of HDL. A recent study has demonstrated that impaired antioxidant HDL function is independently associated with the development of premature acute myocardial infarction [31]. In addition to apoAI, several important enzymes such as PON1, PON3 and PAFAH in HDL play anti-oxidative role. As an essential HDL-associated enzyme, PON1 contributes substantially to the antiinflammatory properties of HDL through preventing LDL oxidation by hydrolyzing oxidized lipids. Our previous study has reported that lower plasma PON1 activity is associated with increased atherosclerotic lesions in apoE −/− mice and in coronary artery disease (CAD) patients [21]. This present study showed that the abundance and activities of PON1 in HDL and plasma of SR-BI −/− mice were markedly reduced, indicating that anti-oxidative activity of SR-BI −/− HDL was hampered. Therefore, oxidative stress was increased in SR-BI −/− mice, which is evidenced by increased plasma MPO activity and HDL-related oxidative property, consistent with a previous study [14].
We also found some other proteins which may contribute to the pro-inflammatory properties of SR-BI −/− HDL, such as acute-phase protein SAA, apoAIV, several cell adhesion molecules, complement C3 and immune response proteins. SAA is demonstrated as a clinically useful marker for inflammation and is strongly associated with increased risk of cardiovascular events [32]. A Fig. 4 The effect of HDL isolated from mice on macrophage inflammatory response induced by oxidized LDL. The mRNA expression of MCP-1 and TNF-α in macrophages were measured by qPCR (a) and the cytokine concentrations in the medium were measured by ELISA (b), n = 3, **P < 0.01 versus ox-LDL group previous report showed that SAA could result in altered metabolic properties of HDL by displacing apoAI, rendering HDL proatherogenic [33]. Human apoAIV has been reported with multiple functions related to lipid metabolism through activation of LCAT and modulation of lipoprotein lipase, and anti-oxidative activity [34]; however, mouse apoAIV, as a positive acute phase protein, was increased during inflammation [35]. Our study identified that SR-BI −/− mice had higher HDLassociated apoAIV levels, supporting that SR-BI −/− HDL in the inflammatory state. The previous study showed that SR-BI −/− mice have a significant lower LCAT activity [36], Our results suggested that decreased activity of LCAT maybe related with apoAI modification and higher  16:155 apoAIV levels in SR-BI −/− mice, although the relative LCAT level of SR-BI −/− HDL was increased in proteomic profile. In addition, proteins involved in cell adhesion like intercellular adhesion molecule 1 (ICAM1), integrin β1, and vitronectin were found in SR-BI −/− HDL, but not in WT HDL. Several proteins with serine proteinase inhibitor domains were detected with increased levels in SR-BI −/− HDL, such as A1AT, antithrombin-III, inter-αtrypsin inhibitor, angiotensinogen, α-2-macroglobulin and murinoglobulin-1, compared to SR-BI +/+ HDL. A recent study of HDL proteomics indicated that protease inhibition ability modulates the anti-inflammatory properties of HDL [37]. Our proteomic data supported that accumulation of the damaged proteins in SR-BI −/− HDL could initiate inflammation, coagulation, and complement activation responses. Furthermore, we observed that HDL isolated from SR-BI −/− mice exhibited less capability of suppressing inflammatory cytokine (MCP-1, TNF-α) expression induced by ox-LDL in macrophages. Taken together, the present study demonstrated that SR-BI deficiency results in a significant alteration in the HDL proteome, which contributes to impair the critical HDL functions, may explain increased risk of cardiovascular disease despite higher HDL-C level in SR-BI −/− mice.
Certainly, mouse does not express CETP, which transfers cholesteryl ester from HDL to VLDL/LDL for further uptake by the liver, it means that HDL metabolism and its related proteins of mice are not exactly the same as those of humans. For example, the changes of apoAIV and A1AT in defective HDL in SR-BI −/− mice were different from those in patients with HDL-associated metabolic diseases [35,38]. However, the previous studies showed that dysfunctional HDL with increased plasma HDL-C levels in SR-BI −/− mice mirrored that in those patients who harbor SR-BI mutations [9][10][11]. Accordingly, our observation made in SR-BI −/− mice studies also supported that the identification of significant proteins in dysfunctional HDL using SR-BI −/− mice as a model, might provide a powerful clue for exploring the biomarkers in clinical diagnosis and therapeutics of dysfunctional HDL-related metabolic diseases. SR-BI −/− mice could be considered as a good preclinical animal model for studying the structure-function relationship of HDL.
Probucol is lipid-lowering and antioxidant drug that has been used in clinic for the treatment of CAD. It could promote RCT, normalize the FC level of HDL [39,40], and reverse female infertility in SR-BI −/− mice [16]. Therefore, we explored the effects of probucol on several critical HDL protein components and HDL dysfunction in SR-BI −/− mice. The current study indicated that accompanied with decreased plasma TC, FC and HDL-C levels, the alterations of some representative proteins of SR-BI −/− HDL could be reversible, and anti-oxidative function was also improved by probucol treatment. These are in agreement with the other studies, that small, protein-enriched, cholesterol-depleted HDL particles are more effective in mediating cholesterol efflux, and with more potent antioxidant, anti-inflammatory, and antiapoptotic capacity [41,42]. The derived protein biomarkers in dysfunctional HDL may eventually help to exploit diagnostics targets and therapeutic strategies for HDLrelated metabolic diseases.

Conclusion
Our proteomic and functional assessments of HDL provide evidence that defective HDL in SR-BI −/− mice have significant decreased levels of apoAI and PON1 proteins and increased levels of pro-inflammatory proteins, such as SAA, apoAIV, and A1AT. The alterations in protein markers of SR-BI −/− HDL result in the dysfunction of HDL including impaired capability in maintaining macrophage cholesterol homeostasis, and the reduced anti-oxidative and anti-inflammatory capabilities. These HDL dysfunctions in SR-BI −/− mice could be improved by probucol treatment. This observations suggest that quantification of apoAI, PON1, SAA, apoAIV, and A1AT may serve as the valuable protein markers for diagnosis and therapeutics of dysfunctional HDL-related metabolic diseases.


Epigenetic basis of hepatocellular carcinoma: A network-based integrative meta-analysis

AIM To identify the key epigenetically modulated genes and pathways in HCC by performing an integrative meta-analysis of all major, well-annotated and publicly available methylation datasets using tools of network analysis. METHODS PubMed and Gene Expression Omnibus were searched for genome-wide DNA methylation datasets. Patient clinical and demographic characteristics were obtained. DNA methylation data were integrated using the Ingenuity Pathway Analysis, a software package for visualizing and analyzing biological networks. Pathway enrichment analysis was performed using IPA, which also provides literature-driven and computationally-predicted annotations for significant association of genes to curated molecular pathways. RESULTS From an initial 928 potential abstracts, we identified and analyzed 11 eligible high-throughput methylation datasets representing 354 patients. A significant proportion of studies did not provide concomitant clinical data. In the promoter region, HIST1H2AJ and SPDYA were the most commonly methylated, whereas HRNBP3 gene was the most commonly hypomethylated. ESR1 and ERK were central genes in the principal networks. The pathways most associated with the frequently methylated genes were G-protein coupled receptor and cAMP-mediated signalling. CONCLUSION Using an integrative network-based analysis approach of genome-wide DNA methylation data of both the promoter and body of genes, we identified G-protein coupled receptor signalling as the most highly associated with HCC. This encompasses a diverse range of cancer pathways, such as the PI3K/Akt/mTOR and Ras/Raf/MAPK pathways, and is therefore supportive of previous literature on gene expression in HCC. However, there are novel targetable genes such as HIST1H2AJ that are epigenetically modified, suggesting their potential as biomarkers and for therapeutic targeting of the HCC epigenome.

INTRODUCTION
Hepatocellular carcinoma (HCC) arises in the context of chronic liver disease, where there is ongoing injury over decades. HCC incidence in North America has been increasing in recent years, in the setting of a higher prevalence of cirrhosis secondary to hepatitis C and fatty liver disease [1] . It is the fifth most common cancer worldwide, and five-year survival is the second worst worldwide among all cancers at 8.9%. HCC is often diagnosed at later stages, and there is an inability to tolerate chemotherapy in patients with cirrhosis [1] . Curative treatment with resection, radiofrequency ablation or transplantation is only possible in early stage disease [2] . When diagnosed at a later stage, the first-line chemotherapeutic agent is sorafenib, which extends survival only by 3 mo [2] . Several other trials of chemotherapeutic regimens have been developed based on studies of genomic and transcriptomic data, with no further improvement in overall survival [3] .
There is therefore a dire need to better understand HCC pathogenesis, elucidate screening biomarkers in patients at risk, and develop more optimal therapeutic agents.
Epigenetic changes are of significant interest to this malignancy, given that it results from mutations accumulating over time with exposure to various insults such as viral hepatitis, alcohol or fatty liver. Epigenetic modifications are heritable states of gene expression without altering DNA sequences [4] . They encompass processes such as DNA methylation, histone modifications, non-coding RNAs and nucleosome positioning. These changes are passed along faithfully to daughter cells during cell division [4] .
Among these, DNA methylation has been the most studied, regulating gene expression through a stable silencing mechanism [5] . Covalent modification by DNA methyltransferases of cytosine residues with methyl groups in CpG dinucleotides occurs preferentially at the 5' end in promoter regions. Transcriptional gene silencing results from this through two mechanisms: steric hindrance of transcription factors being able to access their cognate binding sites on gene promoters [5] , and direct binding of methyl CpG binding domain containing proteins to the methylated DNA causing transcriptional repression [6] . The recent advent of genome-wide methylation analysis has enabled an appreciation of methylation status in genes of interest to cancer: Hypermethylation of tumor suppressor genes, hypomethylation of oncogenes, and methylation of repetitive elements [7] . The extensive reprogramming of the epigenome in cancer has led to a growing interest in epigenetic therapy. Specifically in HCC, aberrant DNA methylation of tumor suppressor gene promoters has been documented [8] . These epigenetic changes have been closely correlated with disease stage and clinical outcome [9] . There has been significant variability in the reported frequency of hypermethylated loci in HCC [10][11][12] . CDKN2A is methylated in 30%-70% of HCCs [10,12,13] , RASSF1A in up to 85% [11,12] , GSTP1 in 50-90% [14] and MGMT in 40% [15] . DNA methylation loci have also been reported as significantly enriched in the signaling networks of cellular development, gene expression, cell death, and cancer [16] .
Our study represents the first comprehensive network-based attempt to integrate all relevant, publicly available, high-throughput genome-wide DNA methylation data to better understand the epigenetic landscape in HCC. Network and pathway analysis tools can enable identification of the most commonly methylated genes across studies and associated pathways, and propose novel treatment options using network-based analysis [17,18] . The goal of this study was to identify key epigenetically modulated genes and pathways in HCC by integrating all major, well-annotated and publicly available methylation datasets datasets using tools of network analysis.
The study workflow is illustrated in Figure 1A. Results were retrieved from both databases: GEO and PubMed. The exclusion criteria listed in Figure 1A were applied to identify papers reporting quantitative results of methylation profile performed on HCC patients and the relative adjacent tissue as control.
Available patient data, including etiology of liver disease (HCV, HBV, alcohol, fatty liver disease) on the basis of which the HCC tumors developed, presence of cirrhosis, the Model for End-stage Liver Disease score (MELD score, an assessment of the severity of liver dysfunction), tumor histology, stage of cancer, alphafetoprotein (AFP) level, overall and recurrence-free survival following treatment were also documented.

Demographic and clinical characteristics
Demographic and clinical patient information pertaining to each dataset are presented in Tables 2 and 3.
Only 6 out of 11 papers (55%) included details regarding presence/absence of cirrhosis, and 10 out 11 papers (91%) provided details regarding the etiology of liver disease. MELD or Child-Pugh score were provided in 27% of papers. Less than half of the selected publications, 5 out of 11 (45%), included details regarding the stage of cancer, although all studies were performed using hepatectomy patient samples. The same trend is identified for information regarding the histologic grade of the tumor (well-, moderately-, or poorly-differentiated tumors). Only 5 out of 11 papers (45%) had alpha-fetoprotein levels available. Overall survival and HCC recurrence statistics as follow-up data were available in 3/11 studies (27%).

Genomic region selection
For the final network-based integrative analysis, we selected 11 datasets ( Table 1). Out of these, raw data from eight studies were available on the GEO website (https://www.ncbi.nlm.nih.gov/geo/). Except for the GSE60753 [25] dataset, for which we have performed our own analysis with R [31] (due to comparison between sample groups required being different from the main paper), we selected the CpG sites or genes reported to be hyper-or hypo-methylated in the corresponding Only genes are reported, without information on CpG sites; 2 Differential methylation data analysis performed in-house.  Table 2 Availability of clinicopathological information on the 11 datasets used for the integrative analysis of genome-wide DNA methylation publications. In 6/11 datasets, CpGs and the mapped genes were selected, whereas in the remaining 5/11 datasets, only the genes without CpG sites were found. Therefore, we separated our analysis into two parts: (1) Taking into account approximately 13500 CpG sites provided or obtained with R analysis ( Figure 1B); and (2) considering only the 4122 differentially methylated genes without information on the corresponding CpG  sites ( Figure 1C). In both cases, the genomic region is considered differentially methylated between HCC tissue and the adjacent non-tumoral sample, if the FDR [32] corrected p-value < 0.01. Furthermore, we filtered out everything that did not satisfy the criteria: ∆ β ≥ 0.20 or ∆β ≤ -0.20, where ∆β = βHCCβadjacent was the difference in methylation between above specified groups. When the CpG sites were considered, the Illumina HumanMethylation450K and 27K platforms were used for mapping to the genes. When multiple sites or genes were found having the same sense of differential methylation, the mean value of ∆β was calculated. The CpGs in the 5'UTR, 1 st Exon, TSS200, TSS1500 or in CpG islands were considered in the promoter and all other CpGs were considered to be in the body of the gene.  [33] , which takes into account the number of focus molecules (input genes) in the network and the total number of molecules in the IPA database that could be included in the corresponding networks.

RESULTS
Based on datasets with known CpG sites, the most frequently hypermethylated genes in the promoter region included HIST1H2AJ, which is a histone protein and SpDYA, a cell cycle regulator known to trigger transition from G1 to S phase. The HRNBp3 gene, an RNA-binding protein, was the most commonly hypomethylated in the promoter region. Further details are provided in Supplementary Tables 1 and 2.

Using the differentially methylated CpG sites in the promoter
Canonical pathways: The most significantly associated pathways with our list of differentially methylated genes are given in Table 4. G-protein coupled receptor signaling, Transcriptional Regulatory Network in Embryonic Cells, cAMP-mediated signaling were the top hits.  (1), HCV (7), alcohol (9), other (10)

Table 3 Clinicopathological information of the individual datasets and methodology used for methylation analysis
differentially methylated in the promoter region, three of them captured our attention: Organismal Development, Organismal Injury and Abnormalities, Cellular Development (Figure 2A), Lipid Metabolism, Small Molecule Biochemistry, Cell Death and Survival ( Figure 2B) and Cell-to-Cell Signaling and Interaction, Drug Metabolism, Small Molecule Biochemistry ( Figure  2C). We noted that networks 1 and 3 were mainly formed by the hyper-methylated genes, whereas network 2 is constituted by both hyper-and hypomethylated genes approximately equally.

Using the differentially methylated CpG sites in the body of the gene
Canonical pathways: G-protein coupled receptor signaling and cAMP-mediated signaling were among the top 20 hits (Table 5), which was similar to the pathways generated by the genes with methylated CpG sites in the promoter region.

Diseases and functions:
Our study shows that DNA methylation in HCC patients leads to the same diseases and functions, regardless of the CpG site position (promoter or body), with cancer and, in particular, abdominal/digestive system cancer among the top listed by IPA.
Networks: From the top detected networks, the Cell-To-Cell Signaling and Interaction, Cellular Assembly and Organization, Cellular Function and Maintenance ( Figure  3A) and Drug Metabolism, Glutathione Depletion in Liver, Small Molecule Biochemistry ( Figure 3B) are closely related to HCC.

Validation using only the reported differentially methylated genes
Having identified the differentially methylated genes corresponding to the CpG sites in the promoter or the body of the genes through this integrative meta-analysis, we wanted to verify how many of these overlapped with the reported differentially methylated genes in those studies that did not include information on CpG sites. The Venn diagram in Figure 4 shows 165 genes reported genes in common with CpG sites in the promoter and 82 genes in common with CpG sites in the body (Supplementary Table 3).

DISCUSSION
The literature on the epigenome in HCC has grown since the advent of tools permitting genome-wide methylation analysis. Epigenetic changes in HCC arise in the context of various etiologies of chronic liver disease, and have been revealed to contribute to tumorigenesis and cancer progression. Therapeutic targeting of HCC has not been as successful as in other malignancies, and requires exploration of a different approach [34,35] . Given that this cancer is driven by various known environmental factors, targeting epigenetic changes in HCC represents a potentially promising therapeutic avenue [36] .
Network-based tools offer a different and unique perspective into the key genes and pathways implicated in disease pathogenesis and progression [37] .
Network-based medicine is critical to a broader understanding of HCC, whose pathogenesis has been difficult to elucidate given the multiplicity of underlying liver disease etiologies [30] . Epigenetic changes impact genetic networks, and a networkbased integrative meta-analysis is ideally suited to integrating and exploring effects of networks on disease pathogenesis [38] . Using IPA, we performed this integrative analysis in order to identify the most commonly aberrantly methylated genes and associated pathways. The most commonly hyper-and hypomethylated genes were identified. These included HIST1H2AJ, which is a histone-coding cell cycle gene previously also identified as hypermethylated in patient lung adenocarcinoma samples [37] and head and neck squamous cell carcinomas [39] . In a study investigating the genetic-and-epigenetic cell cycle network in HeLa cancer cells, methylation of HIST1H2AJ (among other genes) was found to result in cell proliferation and anti-apoptosis through NFκB, TGF-β, and PI3K/Akt/ mTOR pathways [40] . HIST1H2AJ has not previously been highlighted as a gene of interest in HCC, which is a novel finding of our integrative analysis of methylation datasets. This illustrates the power of integrating all available high-throughput data to better understand important genes in cancer. SpDYA, a cell cycle regulator known to trigger transition from G1 to S phase, was differentially hypermethylated. The HRNBp3 gene, an RNA-binding protein, was the most commonly hypomethylated, as had been reported in the integrative analysis of epigenetic data by Song et al [16] in 2012. One would thereby anticipate increased gene expression of HRNBp3 in HCC. These genes with differential methylation have not previously been highlighted in the HCC literature, and serve as potential new biomarkers and therapeutic targets [41] . Using IPA, we then determined the most commonly affected networks in HCC.
G-protein coupled receptor signaling, Transcri-ptional Regulatory Network in Embryonic Cells, cAMPmediated signaling were the top hits, which was in perfect agreement with the work of Song et al [16] , wherein they used IPA to analyze methylation profiling for a set of 27 HCC tumors compared with 20 normal patients. G-protein coupled receptor signalling is common to various principal pathways known to be implicated in HCC, including the PI3K/Akt/mTOR and Ras/Raf/MAPK pathways based on genomic and gene expression analyses. Therefore, our results reinforce the biological rationale of targeting these pathways. We also elucidated the crosstalk between proteins within the networks of interest to HCC. This analysis revealed ESR1 and ERK to be proteins central to they key networks.
A unique aspect of our study was the analysis of methylation at CpG sites in both the promoter and body of genes. Whereas methylation in the gene promoter is known to cause transcriptional repression, methylation in the body has the opposite effect, promoting gene expression. A novel finding was that the genes with the greatest differential methylation in the promoter were the same as those with the greatest differential methylation in the body, further confirming the importance of these genes. We were also able to validate the identity of several genes with data on CpG sites within datasets without such data available.
Limitations of our study include the lack of methylation data on individual HCC samples. Given the relatively recent advent of genome-wide methylation analysis methods, with the earliest dataset in HCC being released in 2010, this analysis was representative of only 354 samples in comparison to a similar number of non-cancerous liver samples. Nonetheless, our study is the largest integrative network-based analysis of DNA methylation in HCC. Clinicopathological characteristics such as grade, stage and survival were available only for half of the datasets, thereby limiting the ability to correlate these data points with the most aberrantly methylated genes. Finally, these data were most representative of hepatitis B and C, as described above.
In conclusion, our integrative analysis of genomewide DNA methylation represents the largest such study in HCC. By integrating all genome-wide DNA methylation data with network-based tools, we have systematically elucidated the landscape of epigenetic DNA modifications in HCC and identified novel potential biomarkers and targetable genes within known pathways of interest to HCC. Therapeutic targeting of the epigenome in HCC is a potential avenue to address this malignancy that arises in the context of various etiologies of chronic liver disease.

Research background
The advent of high-throughput technologies in epigenetics has led to improved characterization of methylation status and its impact on development of Hepatocellular carcinoma (HCC).

Research motivation
HCC is a malignancy that arises in the context of ongoing liver injury from various causes, such as hepatitis B, hepatitis C, alcoholic and non-alcoholic liver disease. Therefore, epigenetic changes are very likely to contribute to the pathogenesis of this malignancy.

Research objectives
We aimed to identify the key epigenetically modulated genes and pathways in HCC by performing an integrative meta-analysis of all major, well-annotated and publicly available methylation datasets using tools of network analysis.

Research methods
PubMed and Gene Expression Omnibus were searched for genome-wide DNA methylation datasets. Patient clinical and demographic characteristics were obtained. DNA methylation data were integrated using the Ingenuity Pathway Analysis, a software package for visualizing and analyzing biological networks. Pathway enrichment analysis was performed using IPA, which also provides literature-driven and computationally-predicted annotations for significant association of genes to curated molecular pathways.

Research results
From an initial 928 potential abstracts, we identified and analyzed 11 eligible high-throughput methylation datasets representing 354 patients. A significant proportion of studies did not provide concomitant clinical data. In the promoter region, HIST1H2AJ and SPDYA were the most commonly methylated, whereas HRNBP3 gene was the most commonly hypomethylated. ESR1 and ERK were central genes in the principal networks. The pathways most associated with the frequently methylated genes were G-protein coupled receptor and cAMPmediated signalling.

Research conclusions
Using an integrative network-based analysis approach of genome-wide DNA methylation data of both the promoter and body of genes, we identified G-protein coupled receptor signalling as the most highly associated with HCC. This encompasses a diverse range of cancer pathways, such as the PI3K/Akt/mTOR and Ras/Raf/MAPK pathways, and is therefore supportive of previous literature on gene expression in HCC. However, there are novel targetable genes such as HIST1H2AJ that are epigenetically modified, suggesting their potential as biomarkers and for therapeutic targeting of the HCC epigenome.

Research perspectives
Our integrative analysis of genome-wide DNA methylation represents the largest such study in HCC. By integrating all genome-wide DNA methylation data with network-based tools, we have systematically elucidated the landscape of epigenetic DNA modifications in HCC and identified novel potential biomarkers and targetable genes within known pathways of interest to HCC. Therapeutic targeting of the epigenome in HCC is a potential avenue to address this malignancy that arises in the context of various etiologies of chronic liver disease.


extfieldsofstudy:
Medicine
Computer Science

text:
Critical assessment of approaches for molecular docking to elucidate associations of HLA alleles with adverse drug reactions

HighlightsAll software assessed could dock Abacavir back into the risk allele structure but not always predict the exact binding mode.Most docking software assessed can distinguish between risk and control alleles.Docking performance can be degraded by using a homology model.Receptor flexibility can negatively affect the docking performance for complex HLA examples.Using AutoDockFR cannot compensate for the added difficulty of docking to the unbound target. &NA; Adverse drug reactions have been linked with genetic polymorphisms in HLA genes in numerous different studies. HLA proteins have an essential role in the presentation of self and non‐self peptides, as part of the adaptive immune response. Amongst the associated drugs‐allele combinations, anti‐HIV drug Abacavir has been shown to be associated with the HLA‐B*57:01 allele, and anti‐epilepsy drug Carbamazepine with B*15:02, in both cases likely following the altered peptide repertoire model of interaction. Under this model, the drug binds directly to the antigen presentation region, causing different self peptides to be presented, which trigger an unwanted immune response. There is growing interest in searching for evidence supporting this model for other ADRs using bioinformatics techniques. In this study, in silico docking was used to assess the utility and reliability of well‐known docking programs when addressing these challenging HLA‐drug situations. The overall aim was to address the uncertainty of docking programs giving different results by completing a detailed comparative study of docking software, grounded in the MHC‐ligand experimental structural data – for Abacavir and to a lesser extent Carbamazepine ‐ in order to assess their performance. Four docking programs: SwissDock, ROSIE, AutoDock Vina and AutoDockFR, were used to investigate if each software could accurately dock the Abacavir back into the crystal structure for the protein arising from the known risk allele, and if they were able to distinguish between the HLA‐associated and non‐HLA‐associated (control) alleles. The impact of using homology models on the docking performance and how using different parameters, such as including receptor flexibility, affected the docking performance were also investigated to simulate the approach where a crystal structure for a given HLA allele may be unavailable. The programs that were best able to predict the binding position of Abacavir were then used to recreate the docking seen for Carbamazepine with B*15:02 and controls alleles. It was found that the programs investigated were sometimes able to correctly predict the binding mode of Abacavir with B*57:01 but not always. Each of the software packages that were assessed could predict the binding of Abacavir and Carbamazepine within the correct sub‐pocket and, with the exception of ROSIE, was able to correctly distinguish between risk and control alleles. We found that docking to homology models could produce poorer quality predictions, especially when sequence differences impact the architecture of predicted binding pockets. Caution must therefore be used as inaccurate structures may lead to erroneous docking predictions. Incorporating receptor flexibility was found to negatively affect the docking performance for the examples investigated. Taken together, our findings help characterise the potential but also the limitations of computational prediction of drug‐HLA interactions. These docking techniques should therefore always be used with care and alongside other methods of investigation, in order to be able to draw strong conclusions from the given results.

and HLA-A allele (A*01:01) which could be assumed to not be associated as they are seen at a high 1 9 4 frequency across European origin (Caucasian) populations (average frequencies obtained from AFND 1 9 5 using gold standard populations [9]: B*57:01 = 0.03, B*07:02=0.10 and A*01:01=0.14). 1 9 6 The allele structures were obtained from the PDB database, where available. Models were made for 1 9 7 those alleles where the structure is not publicly available (Table 1). Target and template sequences 1 9 8 were aligned with ClustalX [42]. For each modelling exercise, ten models were generated using 1 9 9 Modeller 9.9 automodel class and the model with the lowest objective function score was chosen as 2 0 0 the model for docking. This objective function is a score generated from the spatial restraints and the 2 0 1 CHARMM energy terms, reflecting stereochemistry within the structure [43].In these simple cases 2 0 2 there was no need to explore alternative target-template alignments since sequences could be aligned 2 0 3 unambiguously with no insertions or deletions.  The Abacavir risk and control alleles were used to evaluate the homology modelling as the known 2 0 9 structures are available for each allele investigated and so can be compared with the model structure 2 1 0 and docking results. Two models were created for B*57:03, one with one template allele (B5703_m) 2 1 1 and another with two template alleles (B5703_m2). These models could then be compared to the 2 1 2 known structure of B*57:03 (B5703_s), as could the docking predictions, to understand the influence 2 1 3 of these steps when employed in a typical docking protocol. The structure of B*57:01 was also 2 1 4 modelled (B5701_m), from two similar sequences identified to make similar comparisons and 2 1 5 evaluate the reliability of using homology modelling.

1 6
For the Carbamazepine risk associated allele B*15:02, there are four differing residues with control 2 1 7 allele B*15:01. Three of these lie in the peptide binding groove, with only one of these being vital to 2 1 8 the D-pocket architecture, where the Carbamazepine is predicted to bind (pos 156). Only a single 2 1 9 template, the structure of B*15:01, was therefore used to model B*15:02 (B1502_m).

0
The quality of each model was investigated using Ramachandran plots and QMEAN scores. For the 2 2 1 B5701_m and B5703_ m and m2, RMSDs were used to give a measure of how well the models 2 2 2 represent the known structures. The percentage sequence identity for each of the templates used for 2 2 3 each model are shown in S1 In order to identify the flexible side chains for the target structure, the relax function of Rosetta was 2 2 6 used to explore the conformational properties of each residue. By looking at 10 different relaxed 2 2 7 structures for each target, flexible residues can be identified. This allows us to consider the flexibility 2 2 8 of the target structure when using AutoDockFR. When using the ROSIE server, a similar sampling is 2 2 9 incorporated into the docking procedures [34]. Using SwissDock [31], the default parameters search the whole target structure but by setting the 2 3 3 search space parameters, it is possible to restrict binding to perform a local docking assay using the 2 3 4 known binding pocket. Here, the area of interest is the peptide binding groove, including residues 1-2 3 5 180 of the alpha chain. The search space was therefore restricted to this area of interest. The file was 2 3 6 processed to ensure it was in the correct format to be uploaded to SwissDock. This included removing 2 3 7 the ligand and peptide from the structure. The PDB was then passed through the Prepdock server [50] 2 3 8 to prepare the structure for docking. This prepared file was then submitted to SwissDock with the 2 3 9 relevant known drug structures. SwissDock used the ZINC database to obtain the known structures of Abacavir and peptide from 3VRI and this time also removing the water molecules from the target 2 4 4 structure, as Rosetta is unable to natively recognise these residues. The Abacavir drug structure was 2 4 5 extracted from the relevant PDB (3VRI) and converted to SDF format. This was then submitted to the 2 4 6 server to be docked, with the search space specified to centre on the peptide binding groove. This 2 4 7 process was repeated for each of the risk and control alleles to be investigated. Using SwissDock, it can be seen that the B5701_s example gives a docking solution close to the 3 1 4 native position (Fig 2 a-c) with all poses showing binding in the F-pocket. When using AutoDockFR 3 1 5 assuming the receptor to be rigid, similar results were seen to those using SwissDock (Fig 2 d-f). All 3 1 6 poses for the B5701_s risk allele were shown bound in the same pocket as the known binding 3 1 7 position. Little difference was seen between the docking shown for each search space, with only one 3 1 8 pose from the B*57:01 run showing binding outside of the peptide binding groove when the search 3 1 9 space was extended around the whole protein (data not shown). This indicates that the peptide binding 3 2 0 groove is the most favourable region for binding. It can be concluded overall that the Abacavir docks 3 2 1 in the expected binding pocket for these packages, but they cannot predict the exact native pose. As the ROSIE server allows movement of side chains within the target, the modelled structure of the 3 3 2 B*57:01 risk allele target was compared to the B5703_s structure submitted. This resulted in the 3 3 3 Table 2 with B5701_s giving similar average RMSDs to the control alleles.

4
Comparing the lowest RMSDs, the control structures give poses with lower RMSDs than B5701_s. 3 3 5  When docking the Abacavir structure using AutoDock Vina, the exhaustiveness was investigated 3 4 0 using the B*57:01 example (B5701_s) to optimise the protocol before docking the other non-risk 3 4 1 alleles. Using the default exhaustiveness of 8, the RMSD values were shown to be quite variable 3 4 2 ( Table 2). The exhaustiveness was then increased starting at exh=18 and doubling the exhaustiveness 3 4 3 to see the effect. It was found that the exhaustiveness of 112 gave the lowest RMSDs overall and 3 4 4 these did not improve with further increasing of exhaustiveness 3 4 5

Most docking software assessed can distinguish between risk and control alleles 3 4 6
Two methods were used to investigate if the docking software can distinguish between the risk alleles. When the flexible residues were incorporated for the 3VRI docking with Abacavir example, poses 4 1 5 occupied the whole peptide binding groove and did not favour the F-pocket as expected (Fig 6). When 4 1 6 the scores for the poses found inside the F-pocket were compared to those found outside the F-pocket, 4 1 7 it was also seen that these scores did not favour the F-pocket (data not shown). Thus, although the 4 1 8 complex algorithms developed for AutoDockFR have been shown to improve the success of docking 4 1 9 the difficulties of using docking to investigate these challenging HLA-ADR cases as in general 4 4 7 docking for new ADRs will be performed, for example, on structures that have a peptide already 4 4 8 bound but not the drug that is to be docked. Docking Carbamazepine with B1502_m using SwissDock, the poses were predicted to sit in the D-4 5 6 pocket previously identified as of interest. Looking at the docking results (Fig 8a), it can be seen that 4 5 7 Carbamazepine is predicted to bind in the D-pocket of B*15:02, close to Leu156, identified by the 4 5 8 study as important, with only one pose predicted out of this pocket. Using AutoDockFR, the same 4 5 9 general pattern was seen with the D-pocket generally being favoured for B*15:02 (Fig 8b). Using The purpose of this exercise was to compare multiple docking programs to assess their performance 4 7 2 with these challenging HLA-ADR cases. We used the Abacavir example as the benchmark for the 4 7 3


extfieldsofstudy:
Computer Science
Medicine

text:
A new chemoinformatics approach with improved strategies for effective predictions of potential drugs

Background Fast and accurate identification of potential drug candidates against therapeutic targets (i.e., drug–target interactions, DTIs) is a fundamental step in the early drug discovery process. However, experimental determination of DTIs is time-consuming and costly, especially for testing the associations between the entire chemical and genomic spaces. Therefore, computationally efficient algorithms with accurate predictions are required to achieve such a challenging task. In this work, we design a new chemoinformatics approach derived from neighbor-based collaborative filtering (NBCF) to infer potential drug candidates for targets of interest. One of the fundamental steps of NBCF in the application of DTI predictions is to accurately measure the similarity between drugs solely based on the DTI profiles of known knowledge. However, commonly used similarity calculation methods such as COSINE may be noise-prone due to the extremely sparse property of the DTI bipartite network, which decreases the model performance of NBCF. We herein propose three strategies to remedy such a dilemma, which include: (1) adopting a positive pointwise mutual information (PPMI)-based similarity metric, which is noise-immune to some extent; (2) performing low-rank approximation of the original prediction scores; (3) incorporating auxiliary (complementary) information to produce the final predictions. Results We test the proposed methods in three benchmark datasets and the results indicate that our strategies are helpful to improve the NBCF performance for DTI predictions. Comparing to the prior algorithm, our methods exhibit better results assessed by a recall-based evaluation metric. Conclusions A new chemoinformatics approach with improved strategies was successfully developed to predict potential DTIs. Among them, the model based on the sparsity resistant PPMI similarity metric exhibits the best performance, which may be helpful to researchers for identifying potential drugs against therapeutic targets of interest, and can also be applied to related research such as identifying candidate disease genes.

Background
A key component in the drug discovery process is to accurately identify the drug-target interactions (DTIs). Traditionally, experimental determination of DTIs is both costly and time consuming. In addition, to fully explore the growing chemical and genomic (for drug targets) spaces being discovered, it becomes impractical to experimentally validate all possible combinations of drug-target pairs. Thus, effective computational algorithms used for predicting potential DTIs are increasingly in demand. Typically, docking simulation is often used to probe the interactions between a series of small molecules and a target under study [1] at a molecular level. However, docking methods require accurate three-dimensional structures of target proteins, making such studies challenging for membrane proteins due to the challenge of protein crystallization. Quantitative structure-activity relationship (QSAR) is another method to depict possible DTIs. However, QSAR typically requires molecular structures with similar scaffolds [2][3][4] for stronger performance. Nowadays, the technology advancement of next-generation sequencing (NGS) and small molecule high-throughput screening (HTS) is accelerating the identification of potential therapeutic targets and drug compounds, which presents great challenges as well as opportunities for chemogenomic research to explore both chemical and genomic spaces simultaneously. In line with this, Yamanishi et al. [5] proposed a bipartite graph learning method correlating the chemical/genomic spaces with the interaction space (i.e., pharmacological space) for predicting potential DTIs, which was followed by several algorithms with improved performance. For example, Bleakley et al. [6] proposed a novel supervised inference method to predict unknown DTIs by using several bipartite local models (BLM). Specifically, BLM transformed the edge prediction problem into the binary classification problem of points with labels. van Larrhoven et al. [7] used a regularized least squares algorithm combined with the Gaussian interaction profile (GIP) kernel calculated only from the topological information of the drug-target network for inferring DTIs. Mei et al. [8] introduced a neighbor-based interaction-profile inferring method and integrated it into BLM, enabling the model for predicting new drugs/targets. Hao et al. [9] employed a nonlinear kernel diffusion (KF) technique to infer DTIs. Liu et al. [10] proposed a neighborhood regularized logistic matrix factorization (NRLMF) algorithm to partly overcome the imbalanced problem in the DTI prediction process. Later on, Hao et al. [11] designed a dual-network integrated logistic matrix factorization (DNILMF) technique by incorporating an idea for modeling social ensemble into the DTI prediction model. Recently, Olayan et al. [12] proposed a method (called DDR), which is based on heterogeneous graph including known DTI network and multiple similarities from both targets and drugs, to predict unknown DTIs by using Random Forest as a classifier. By adding a heuristic selection of similarity matrices and nonlinear KF technology, DDR outperformed other state-of-the-art priors [12]. Additionally, many other DTI prediction algorithms developed previously can be found in the reviews [13][14][15][16][17].
Among popular DTI prediction algorithms, the most reliable and accurate ones are those based on similarities. However, the used similarity information is derived either from protein sequences or from drug structures. Despite of the importance of the DTI graph, little studies considered using its similarity information as the main source when building the model with the exception of previous work [7]. In fact, the DTI bipartite network itself contains extremely important information, which will be beneficial to the model performance. The success of recommender system in e-commence has provided a proof of concept [18][19][20][21], which explores the bipartite network solely. Inspired by this technology, we in this work make an effort to apply and extend it for DTI predictions. Herein, we adopt a technology, called neighbor-based collaborative filtering (NBCF), which is one of the most successful technologies in the community of recommendations. For applying NBCF to DTI predictions, a fundamental step is to accurately compute the pairwise drug similarities based on the drug interaction profiles (DIPs) with targets in the bipartite interaction network, rather than based on drug structures as used in the previous studies. In fact, a similar idea has been reported whereas the protein similarities were measured based on their associated ligands but not based on amino acid sequences [22,23]. With the DIPsbased drug similarities, an intuitive model is built using NBCF by making the following assumption: if drug A and drug B are highly similar (again, as indicated by similarity from DIPs), and if drug A interacts with the current target, then drug B has a high probability of interacting with the same target, though there may be exceptional cases [22,23]. However, it is well-known that the experimentally validated interaction information is extremely limited compared to the whole drug-target interaction space, which will introduce noise when computing similarity from such a sparse network (sparseness, defined as the number of links divided by the total number of possible target-drug pairs) using the conventional similarity calculation methods. To tackle this challenge, we in this work propose three strategies to remedy the issue, i.e., by designing a new similarity metric to mitigate noise, performing low-rank approximation (LRA) of the original prediction scores, and incorporating the auxiliary information into the model.
It is critical to select an appropriate evaluation method in order to assess the strength of a developed DTI prediction algorithm as well as to identify rooms for further improvement. Instead of adopting the commonly used evaluation metrics [i.e., area under curve (AUC) and area under precision-recall (AUPR) curve], we introduce a recall-based metric, namely mean percentile ranking (MPR), which is under-studied in DTI predictions [17] but routinely used in the recommender system studies [18,24]. The reason for selecting MPR as the evaluation criteria is that one only knows about the one-class experimentally validated information (i.e., a drug interacts with a target, which is considered as the positive information) but does not know about the negative information (i.e., a drug does not interact with a target) due to the lack of comprehensive experimental data on a drug-target pair. Thus, a recall-based metric is suitable to such a scenario. Finally, we validate our method in three large publicly available datasets and compare the proposed algorithm with the prior art based on MPR. We conclude that the proposed NBCF algorithm with the improved strategies is both effective and computationally efficient for DTI predictions, which outperforms the previously developed algorithm for identifying potential drugs against therapeutic targets under a study.

Datasets
Three large benchmark datasets were used to evaluate the current proposed NBCF algorithm for DTI predictions. The first dataset (denoted by DATASET-H) was derived from our previous work [11], which consists of 733 targets and 829 drugs with 3688 known DTI pairs. This dataset was obtained based on the DrugBank database [25] followed by several pre-processing operations including removing duplicated molecules, mapping to unique identifiers and a few other steps as described previously [11]. The second dataset (denoted by DATA-SET-K) was retrieved from the study of Kuang et al. [26], which includes 809 targets and 786 drugs forming 3681 known DTI interactions. In this dataset, the drugs were approved by FDA, assigned with at least one ATC code, and the drug data were deposited in the KEGG database [27]. The third dataset (denoted by DATASET-Y) includes 664 targets and 445 drugs with 2926 experimentally validated interactions, which was studied by Yamannishi et al. [5]. Specifically, DATASET-Y was retrieved from multiple databases including KEGG BRITE [27], BRENDA [28], SuperTarget [29] and DrugBank [25]. All benchmark datasets used in this work consist of three matrices: (1) drug-target interaction (adjacency) matrix, denoted by Y ∈ R M×N with M targets and N drugs; (2) target sequence similarity matrix, denoted by S T ∈ R M×M , calculated from target sequences; and (3) drug structural similarity matrix, denoted by S D ∈ R N ×N , computed from drug chemical structures. Matrix Y is often filled by binary numbers, where Y ij = 1 if target i is targeted by drug j validated by the previous experiment, and otherwise Y ij = 0 (indicating that drug-target interaction information for the specific pair is unknown). Table 1 shows the benchmark datasets as well as corresponding properties used in this work.

Workflow of the proposed algorithm
The task of DTI predictions considered in the work is to identify drugs that have larger possibilities of interacting with the targets of interest. Specifically, given a series of targets and drugs, as well as a very small number of known (experimentally determined) interactions, a bipartite network was constructed as shown in Fig. 1a. The bipartite network was converted into an adjacency matrix (also called drug-target interaction matrix), which is very sparse due to the extremely low number of experimentally validated interactions compared to the whole drug-target pair space (shown in Fig. 1b). While "1" is used to indicate a known (positive) interaction, "0" is used to indicate that it is unknown whether the corresponding drug and target interact with each other, because an experiment has not been performed. Based on the sparse interaction matrix, we proposed to use NBCF to infer the potential interactions for those drug-target pairs labelled as 0 s. The development of NBCF was based on a hypothesis that if a query target T 1 has been reported to interact with drugs of D 1 , D 2 and D 3 that are very similar to D N , then T 1 has a large probability for interacting with D N . While it is true that chemicals with similar structures do not always exert the same biological properties depending on the similarity degree (e.g., activity cliff) [22,23,30], chemical similarity is still a significant principle used when searching for compound candidates for the desired biological activity in drug design and development [31]. Evidently, the key step of NBCF is to accurately assess the pairwise similarity between drugs. Being different from the previous algorithms such as BLM [6] whereas the prior similarity information from drug structures and protein sequences, such as S T or S D , was used as the input (kernel) matrix of support vector machine (SVM) and a conventional binary classification was performed with the fixed regularization parameter C of 1, the NBCF technique proposed in the work mainly depends on the similarity information calculated from DIPs in the drug-target interaction matrix, which is denoted by S DIP as shown in Fig. 1c. It is well-known that there are multiple where S cos jj ′ denotes the COSINE similarity between drug j and drug j' with the range from − 1 to 1, and M is the number of targets. The TANIMOTO similarity (coefficient) is defined as follows: where S tan jj ′ denotes the TANIMOTO similarity between drug j and drug j′ with the range from 0 to 1. In addition to the two commonly used similarity calculation methods based on the binary data, we also proposed to use positive pointwise mutual information (PPMI) to measure the similarity between a drug pair. The PPMI approach, which is under-studied in DTI research, has been reported to be a similarity metric, which can mitigate the data sparsity issue to some extent [32]. In the sparse DTI network, the PPMI similarity is defined as follows: where the probabilities P Y .j , Y .j ′ and P Y .j are estimated empirically as follows: where co Y .j , Y .j ′ is the number of times that drugs j and j′ co-occur calculated by summing both co-occurred ones and zeroes in the matrix Y, and N is the number of drugs.
It should be noted that S ppm jj ′ is non-negative by replacing Fig. 1 Workflow of the proposed NBCF algorithm with strategies designed for improving DTI predictions negative values to zeroes, and hereby the base 2 logarithm was used in Eq. (3). After yielding three similarity matrices ( S cos jj ′ , S tan jj ′ and S ppm jj ′ ) calculated from the interaction matrix Y solely, the proposed NBCF algorithm was used to calculate the prediction scores (Fig. 1d, e), which is defined as follows: where Ŷ ij denotes the predicted interaction scores between the target i of interest and the query drug j. S jk denotes the similarity values (i.e., those from either S cos jj ′ , S tan jj ′ or S ppm jj ′ ) between the query drug j and drugs with known interaction information for the current target i. It should be pointed out that, while being simple and intuitive, the proposed algorithm is effective and computationally efficient for DTI predictions due to the model-free property similarly as reported by the previous studies [22,23]. In fact, a similar algorithm has been successfully applied in the field of recommender systems [19,21,33]. However, it should be emphasized that the DTI interaction matrix is extremely sparse, therefore the calculated similarity matrix may include noise, which will decrease the model performance [19]. Thus, we proposed three strategies in the work to overcome such a dilemma. Strategy 1: we designed a similarity calculation algorithm, which is immune to the data sparsity issue to certain degree, with the final generated similarity S ppm jj ′ . As shown in Fig. 1d, the final prediction scores were obtained directly by using Eq. (6) based on S ppm jj ′ . Moreover, if the commonly used similarity calculation algorithm is used with the generated matrix (i.e., S cos jj ′ or S tan jj ′ ), despite that the prediction scores are calculated by using Eq. (6), the scores would be considered as temporary ones as denoted by Ŷ t ij as shown in Fig. 1e, due to that Ŷ t ij may be sub-optimal because of the noisy similarity information. Thus, we proposed two additional remedy strategies (i.e., Strategy 2 and Strategy 3) to improve the performance on the basis of temporary prediction scores, Ŷ t ij . Strategy 2: as reported in the community of recommender systems [19], LRA of original prediction scores can help to partially mitigate noise. Thus, we incorporated this technique into the DTI prediction domain. Specially, the singular value decomposition (SVD) algorithm as one of the most popular LRA techniques was adopted to factorize the temporary score matrix, Ŷ t ij (Fig. 1f ). The final prediction scores were formed according to the following equation: where U ∈ R M×R is the left singular vector matrix with rank R (empirically set to 100), S ∈ R R×R is the diagonal matrix, and V ∈ R N ×R is the right singular vector matrix with V T denoting the transpose of V. Strategy 3: while S DIP remains as one of the key components of NBCF derived from the DTI network, the auxiliary similarity information (e.g., S T and S D ) may be attributed as complementary sources that are beneficial to the model performance. In fact, several previous studies have demonstrated its effectiveness [11,34]. Therefore, we also explored to include auxiliary information in the NBCF method for the final DTI predictions as defined below (Fig. 1g): where α , β , and γ are the smoothing coefficients (empirically set to 0.025, 0.95 and 0.025, respectively).

Evaluation method
In this work, we used tenfold cross-validation to evaluate the proposed algorithm. Specifically, we removed randomly a subset of 10% of the links (known interaction pairs) in the drug-target interaction matrix Y as the test set and trained models on the remaining links (i.e., 90% of the known interaction pairs). In addition, we ensured each drug has at least one interaction with a target (and vice versa that each target has at least one interaction with a drug as well) similarly to the previous work [35]. We adopted a recall-based evaluation metric, MPR [18,24], to evaluate the algorithm performance. In detail, for each target i in the test set, we generated a ranked list of potential drugs, sorted by a decreasing order according to the final prediction scores for the potential interaction between target i and each of the drugs in the dataset. Let rank ji denote the percentile ranking (PR) of target i for drug j. This way, at rank ji = 0%, drug j is predicted as the drug with the highest probability of interacting with target i, while at rank ji = 100%, drug j is predicted as the drug with the lowest probability of interacting with target i. Herein, the definition of MPR is described as follows: where N test T denotes the number of targets in the test set, and R i is computed as follows: where N test D denotes the number of drugs in the test set for the current target i. It should be pointed out that the lower MPR is, the more desirable performance the model exhibits, as a lower MPR value indicates the drug-target pair is predicted as interacting with each other with a higher possibility. Evidently, randomly generated lists have an expected MPR of 50% [24]. Using this metric, one can obtain a recommended list of candidate drugs, with top predictions recommended to be given higher priority for experimental validation.

Properties of benchmark datasets
We validate our algorithm using three benchmark datasets (Table 1). (1) DATASET-H: in this dataset which was obtained from our previous work [11], there are 733 unique targets and 829 unique drugs extracted from the DrugBank database following several preprocessing steps. On average, DATASET-H has about 4 known targets for each drug and 5 drugs for each target. Among them, looked from the drug side, the minimum and maximum number of interacted targets are 1 and 48, respectively. From the target end, the minimum and maximum number of interacted drugs are 1 and 75, respectively. The sparsity value (calculated from known interactions divided by the totally possible interaction pairs between drugs and targets; the lower the value, the sparser the dataset is) is 0.006, indicating the dataset is very sparse. (2) DATASET-K: the dataset was retrieved from the publication of Kuang and co-workers [26]. This dataset is similar to DATASET-H, but has more targets than drugs. This dataset also has the sparsity value of 0.006, making the DTI predictions extreme difficult.
(3) DATASET-Y: this dataset is a subset of the previous work with the largest number of possible interaction pairs published by Yamanishi and co-workers [5]. Similar to DATASET-K, this dataset also has more targets than drugs. Compared to the first two datasets, the sparsity value of DATASET-Y is relative higher (0.010) indicating it is relatively less (but still very) sparse and has more known interactions within the dataset. In summary, all these three benchmark datasets have a very low sparsity value leaving a larger room for challenging the algorithms for DTI predictions.

Results of the proposed algorithm
In this section, we evaluate the proposed NBCF algorithm for predicting DTIs using the three extremely sparse benchmark datasets. As shown in Table 2, in Strategy 1 (i.e., results are totally based on S DIP as shown in Fig. 1c-e), results based on PPMI give MPR values as of 0.054, 0.049 and 0.020 for DATASET-H, DATASET-K and DATASET-Y, respectively. COSINE-based MPR values are 0.081, 0.068 and 0.037 for the same datasets, respectively, while TANIMOTO-based MPR values are 0.092, 0.070 and 0.035. From these results, we conclude that the proposed NBCF algorithm has generated promising results which largely outperform the random recommendation accuracy (i.e., 0.5) in terms of MPR [18,24], and evidently, PPMI-based NCBF significantly outperforms both the COSINE-based and TANIMOTObased counterparts (P < 0.01, t test). The observation is not surprising because the similarity information used in the PPMI-based NBCF technique is intentionally designed for overcoming noise from the sparse DTI network, while the NBCF methods based on the COSINE and TANIMOTO similarity metrics exhibit the suboptimal results due to the noise-prone properties in such commonly used calculation methods. It should be emphasized that while results from the PPMI-based NBCF algorithm are used as the final prediction scores (Fig. 1d), the COSINE and TANIMOT-based ones are considered as the temporary results, denoted by Ŷ t ij as shown in Fig. 1e, which can be further improved by our proposed strategies as described in the following. Since the previous study has reported that the LRA operation of original prediction scores can reduce noise to some extent [19], we thus adopt one of the most popular LRA techniques (i.e., SVD) to factorize Ŷ t ij and yield the final prediction scores according to Eq. (7), which belongs to Strategy 2. As shown in Table 2, the strategy largely enhances the performance for COSINE-based and TANIMOTO-based NBCF. For example, in DATA-SET-H, COSINE-based NBCF improved MPR from the where the performance actually decreases when LRA was applied. This is because similarity based on PPMI used in NBCF has already successfully reduced noise from the sparse DTI network, thus an extra LRA operation by using SVD might be "over killing" and may even affect the results adversely. Moreover, it is interesting to note that in DATASET-K, COSINE-based NBCF in Strategy 2 generated comparable results with PPMI-based one in Strategy 1. Though one of the key ideas of NBCF is to accurately construct the similarity matrix (i.e., S DIP as shown in Fig. 1c) only from the DTI profiles (Fig. 1b), the auxiliary similarity information, such as S T and S D used in this work, may be beneficial to the model performance by incorporating complementary information appropriately, as expected for models from COSINE-based and TANIMOTO-based NBCF especially. Thus, in Strategy 3, we explore such auxiliary information by adding them into the original NBCF model with a similar approach used by the previous studies [11,34] (Eq. 8). As shown in Table 2, it is evident that both COSINE-based and TANIMOTO-based NBCF models exhibit enhanced performance in all benchmark datasets, with the exception that in DATASET-H, COSINE-based NBCF gave slightly lower performance than that in Strategy 1. However, PPMI-based NBCF in Strategy 3 does not show appreciation for such auxiliary similarity information at all. On the contrary, decreased performance is observed with the incorporation of S T and S D , which indicates that the NBCF model on the basis of PPMI can generate the most optimal performance, while extra operations may have an adverse effect on the model. In summary, we conclude that the proposed strategies are undoubtedly playing a central role in improving the DTI prediction performance based on the NBCF model. Among them, PPMIbased NBCF gives the best results in all three benchmark datasets due to the well-designed similarity measurement method, which can effectively tackle the sparsity issue in the DTI network. Moreover, both the LRA operation and incorporation of auxiliary information are helpful to enhance the performance for models that are based on the commonly used similarity metrics such as COSINE and TANIMOTO. Figure 2 shows the corresponding boxplots of all these results.

Comparison to counterpart and further consideration
We compared the proposed NBCF algorithm to DT-Hybrid proposed by Alaimo and co-workers [35]. We select DT-Hybrid for comparison because (1) both NBCF and DT-Hybrid are derived from network based recommendation technology [19,20,33]; (2) both algorithms adopt a recall-based metrics; and (3) they are both effective and computationally efficient for DTI predictions. For DT-Hybrid, we adopt the default parameters according to the reported values (i.e., lambda set to 0.5 and alpha set to 0.4) [35]. As shown in Table 2, in all three datasets, PPMI-based NBCF in Strategy 1, and both COSINE-based and TANIMOTO-based NBCF models in Strategy 2 generated much better results than those from DT-Hybrid. Similarly, models from COSINE and TANIMOTO in Strategy 3 consistently outperform those from DT-Hybrid. Therefore, all the results indicate that our proposed algorithm with the improved strategies demonstrated stronger prediction ability for inferring potential DTIs. Though the NBCF algorithm combined with similarity from either PPMI, COSINE or TANIMOTO proposed in this work has been successful for DTI predictions, we were interested to explore the effect of other similarity methods on the model performance. Since the GIP kernel was reported to be a useful similarity metric for predicting potential DTIs in the previous work [7,9,11], we performed an experiment based on GIP. However, no satisfied results were obtained in all three benchmark datasets in terms of MPR. When we tested another similarity metric called DICE coefficient, results showed similar trend with those based on COSINE and TANIMOTO. Furthermore, we experimented the proposed algorithm with the IC dataset from the previous study [5], and noticed that the PPMI-based model still exhibits the best performance. To further validate the model effectiveness, we also performed fivefold cross-validations, which generated similar results as those from the tenfold cross-validation. It should be noted that the current work mainly focuses on inferring the potential drugs for interesting targets. In fact, the inverse operation (i.e., inferring the potential targets for interesting drugs) is also possible, which will be further explored in the future. Moreover, we plan to improve the current algorithm to make it scalable to larger datasets, and suitable to the new targets (or new drugs) scenarios [8,10,11,36,37].

Conclusions
In this work, we propose a straightforward yet effective and computationally efficient algorithm, NBCF, for inferring potential DTIs. For overcoming data sparsity inherently existing in the known DTI network, we designed three strategies to tackle the difficult issue. In Strategy 1, we propose to use a sparsity resistant similarity metric, PPMI, to measure the correlation between drugs from the DTI network solely, which as a result exhibits the best performance in the current work. In Strategies 2 and 3, we apply the low-rank approximation technique and incorporate additional auxiliary similarity into noise-prone models (i.e., COSINE-based NBCF and TANIMOTO-based NBCF) respectively, which have been shown to enhance the prediction accuracy to identify drug candidates for therapeutic targets.