Individualized Drug Repositioning For Rheumatoid Arthritis Using Weighted Kolmogorov-Smirnov Algorithm.

Rheumatoid arthritis (RA) is a chronic autoimmune disease characterized by joint synovial tissue inflammation associated with the disability of affected joints. 1 Patients with RA have an increased mortality, and the expected survival of RA patients is likely to decrease 3–10 years. 2 Early diagnosis and effective therapy are critical to prevent joint deterioration and unfavorable disease outcome. Currently, the treatment of RA primarily rests on the use of disease-modifying antirheumatic drugs, and has improved outcomes in RA patients significantly. 3 – 5 Despite significant therapeutic advances in improving the lives of RA patients, RA remains a hard clinical problem because of the accumulated and persistent disease. 1 The administration of RA patients needs new drugs for preventative or curative therapies. Presently, drug repositioning has drawn broad attention from the pharmaceutical companies and research institutes. Relative to the traditional drug development process, drug repositioning replenishes the drying out drug pipelines by reusing marketed drugs and clinical candidates for new uses, such as treating another disease. 6 These repositioned drugs with known bioavailability, safety profiles and well-characterized pharmacology can enter clinical trials for alternative indications more rapidly and less risk. 7 Currently, multiple computational approaches have been established for drug repositioning. 8 – 10 Connectivity map (CMap) has been widely used in drug repositioning by measuring the similarity in gene expression profiles between compounds in mammalian cell lines. 11 While current methods mainly focus on discovering drug candidates targeting huge populations, 12 , 13 personalized therapeutic decisions are scarce. Human genetics provides insight into disease pathogenesis and guides drug discovery for complex traits. 14 , 15 A large body of evidence points out the influence of inherited genetic factors on both susceptibility and resistance to the disease. 14 – 16 Currently, high-throughput genome-wide association studies have resulted in a paradigm shift in the way that researchers treat complex diseases. Several lines of evidence has revealed numerous genes influencing the likelihood of developing RA by genome-wide analysis. 17 , 18 While most of the biological functionality of the cell arises from complex interactions among genes, and interpreting the consequences on a pathway level has more powerful in understanding how gene activity perturbations account for disease. 19 , 20 Moreover, personalized pathway analysis has been proposed to perform personalized or customized interpretation of disease data, 21 making it possible to develop personalized therapeutic decisions. Here, we attempted to predict candidate drugs for RA treatment from CMap database by mining the similarities of pathway aberrance induced by disease and various drugs, on a personalized or customized basis. Our study will provide novel insights into personalized drug discovery for RA.

In the present study, 15 healthy subjects were denoted as nRefs (reference) of 18 subjects diagnosed RA. The genes were subjected to quantile normalization to evaluate the gene-level statistics. Meanwhile, a total of 281 pathways were screened from the KEGG pathway database after data preprocessing. By iPAS algorithm, we obtained the pathway aberrance scores of individual RA subjects. Using the mean value of iPAS as the pathway aberrance level of RA, 170 pathways were up-regulated and 111 pathways were down-regulated in RA. The conditions of the altered pathways in all RA subjects were elucidated by individual pathway analysis. The pathway statistics for each individual was tested by wilcoxon-test and the p-value was adjusted by FDR. Under p-value <0.05, a total of 4 pathways were regarded as differential pathways ( Table 1 ). In our study, cardiac muscle contraction pathway showed the highest altered frequency, which altered in 14 of 18 RA subjects, followed by amoebiasis, amino sugar and nucleotide sugar metabolism, and protein processing in endoplasmic reticulum pathway. Table 1 The Differential Pathways In Rheumatoid Arthritis Based On The Individualized Pathway Aberrance Score Category Pathway p Altered Frequency Hsa00520 Amino sugar and nucleotide sugar metabolism 0.019 11 Hsa04260 Cardiac muscle contraction 0.021 14 Hsa04141 Protein processing in endoplasmic reticulum 0.024 10 Hsa05146 Amoebiasis 0.047 12 The Differential Pathways In Rheumatoid Arthritis Based On The Individualized Pathway Aberrance Score Similar to the identification of pathway aberrance in RA, we analyzed the function aberrance induced by each drug in CMap using iPAS algorithm. For each drug, we calculated the specific gene expression profiles to detect the drug-induced iPAS status of each pathway, using untreated human cells as nRefs (reference). After data preprocessing, we obtained 888 drugs from CMap database and 281 pathways from KEGG database for further analysis. After drug iPAS analysis, we obtained a 281 pathways × 888 drugs matrix. Based on the iPAS matrix of RA and drugs, a KS-weighted ES analysis was implemented to compare similarities between RA-induced pathway aberrances and drug-induced pathway aberrances. Then, a non-parametric, rank-based procedure was employed to yield a rank list to select the drug candidates that might mimic or treat RA. The pathway profiles for 888 drugs were merged into a specific PRLs by ranking the aberrant pathways. Then, pathway-oriented drug matrix was converted into an RA subject-oriented drug matrix by measuring the ES values. After that, each drug was given a specific ES value for each RA subject. A drug with high ES values indicated that the drug showed related genomic response to RA. The candidate drugs for RA subjects were assessed by applying the rank-based procedure, and the drug-subject matrix was sorted row-wise to select RA-related drugs from the most inhibiting one to the most mimicking one. Here, we identified 9 candidate drugs for RA treatment, as shown in Table 2 . Also, nine good mimics were identified based on the similarity of pathway profiles to disease ( Table 3 ). Some candidates with current indications to treat other diseases might be repurposed to treat RA. Table 2 The Candidate Therapeutic Drugs Of Rheumatoid Arthritis Candidate Drug Mean ES Rank Score CAS Number Stage Current Drug Indication Suloctidil 0.245 0.783 54063-56-8 Withdrawn Vasodilator Prenylamine 0.245 0.783 390-64-7 Withdrawn Vasodilator Mebendazole 0.253 0.780 31431-39-7 Approved Anthelmintic Danazol 0.247 0.778 17230-88-5 Approved Endometriosis Piperlongumine 0.249 0.775 20069-09-4 Investigated Anticancer Harmine 0.254 0.773 442-51-3 Investigated Osteogenic differentiation Celecoxib 0.247 0.770 169590-42-5 Approved Arthritis Sanguinarine 0.246 0.769 2447-54-3 Investigated Anticancer Parthenolide 0.250 0.767 20554-84-1 Investigated Anti-inflammatory Abbreviation: ES, enrichment score. Table 3 The Good Mimics Of Rheumatoid Arthritis Candidate Drug Mean ES Rank Score CAS Number Stage Current Drug Indication Camptothecin 0.608 0.088 2114454 Experimental Anticancer Doxorubicin 0.583 0.104 23214-92-8 Approved Anticancer Daunorubicin 0.582 0.132 20830-81-3 Approved Anticancer Ampicillin 0.531 0.221 69-53-4 Approved Anti-infection Cloxacillin 0.518 0.241 61-72-3 Approved Anti-infection Streptomycin 0.500 0.271 57-92-1 Approved Anti-infection Cortisone 0.501 0.271 53-06-5 Approved Anti-inflammatory Estradiol 0.499 0.271 50-28-2 Approved Estrogen deficiency Mitoxantrone 0.487 0.278 65271-80-9 Approved Anticancer Abbreviation: ES, enrichment score. The Candidate Therapeutic Drugs Of Rheumatoid Arthritis Abbreviation: ES, enrichment score. The Good Mimics Of Rheumatoid Arthritis Abbreviation: ES, enrichment score.

Here we retrieved the transcriptome data of RA from ArrayExpress database ( http://www.ebi.ac.uk/arrayexpress/ ), under the accession number of E-GEOD-15573. 17 In the study of Teixeira et al, 17 a complete genome-wide transcript profiling of peripheral blood mononuclear cells from 18 RA patients and 15 controls was conducted using the Illumina Human-6v2 Expression BeadChips. The raw data and the annotations were obtained from the manufacturer’s documents, and the probes were re-annotated to genes symbols. The CMap is a collection of genome-wide transcriptome data from cultured human cells treated with bioactive small molecules and simple pattern-matching algorithms that discover connections between drugs, genes expression changes and disease phenotypes. 11 The CMap contains more than 7000 gene expression profiles for 1309 compounds. All raw data and the annotations were downloaded from the CMap and the expression values for all samples were calculated by affy package R 22 with MAS 5.0 normalization. The probes were re-annotated to genes symbols using Brainarray CDF packages. 23 Finally, the samples corresponding to the same drugs were merged, and the gene-drug matrix was obtained for subsequent analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/ ) is a knowledge database for systemic analysis of gene functional information. 24 All 300 human pathways (covering 6919 genes) were downloaded from KEGG database. In the present study, we retrieved transcriptome data of RA, CMap data and KEGG pathway data. Prior to analysis, a preprocessing procedure was performed. Firstly, the genes obtained from three data were intersected to gain the common genes. Then, KEGG pathways containing 100 genes were removed, because pathways with too many genes might be too complex to understand, and pathways with too few genes may not have sufficient biological content. Finally, a total of 6919 genes, 888 drugs, and 281 pathways were selected for subsequent analysis. In this section, the pathway levels in each sample were calculated using iPAS algorithm by making use of the accumulated normal data. In this study, 15 normal control samples were combined and regarded as references (nRef). For individual RA cases, an uniformly normalization was performed after combining the single RA data with all nRef samples. The gene expression value of individual RA sample was standardized by mean and standard deviation (SD) of the reference. For each gene a , we calculated the gene expression level as follows: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document} $${{\rm{Z}}_{\rm{a}}}{\rm{ = }}{{{{\rm{g}}_{{{\rm{R}}_{\rm{a}}}}}{\rm{ - mean}}\left({{{\rm{g}}_{{\rm{nRef}}}}} \right)} \over {{\rm{SD}}\left({{{\rm{g}}_{{\rm{nRef}}}}} \right)}}$$ \end{document} where g Ra represented the expression value of gene a in an individual RA subject, mean(g nRef ) stood for the mean expression value of genes a in all nRef cases, SD(g nRef ) stood for SD of the reference. To evaluate the iPAS by nRef, Average Z algorithm was employed, which presented well in highlighting pathway aberrance and in revealing clinical importance. 21 A vector Z = ( z 1 , z 2 , …, z n ) represented the expression status of a pathway, where z a stood for the standardized expression value of the a -th gene, and n represented the gene number in the specific pathway. The iPAS value of a pathway was calculated as follows: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document} $${\rm{iPAS = Average}}\,\,{\rm{Z = }}{{\sum\nolimits_{\rm{a}}^{\rm{n}} {{{\rm{Z}}_{\rm{a}}}} } \over {\rm{n}}}$$ \end{document} Then, the expression matrix (281 pathways × 18 RA subjects) was obtained for each pathway in each individualized RA subject from the nRef. The pathway statistics for each individual RA subject was calculated by wilcoxon-test and false discovery rate (FDR) was used to adjust the p-value. The pathways with p-value <0.05 were defined as differential pathways. The CMap contains more than 7000 gene expression profiles for 1309 compounds. Each drug presents a specific drug-induced gene expression changes of human cells, enabling us to identify the pathway aberrance. After data preprocessing, a total of 6919 genes, 888 drugs, and 281 pathways were selected for subsequent analysis. To identify the drug-induced pathway aberrance, iPAS algorithm was utilized to estimate the pathway levels. For each drug, we calculated the specific drug-induced iPAS status of each pathway using Average Z algorithm. Then, the expression matrix (281 pathways × 888 drugs) was obtained for each pathway in each drug. After the above treatment, we obtained RA-induced pathway aberrances and drug-induced pathway aberrances, respectively. Then, we systematically estimate the similarities between RA-induced pathway aberrances and drug-induced pathway aberrances using a to select drugs that might mimic or suppress RA. Prior to similarity analysis, we firstly built Prototype Ranked Lists (PRLs) 25 by merging all the samples corresponding to the same drug, after converting iPAS values to ranks (the iPAS value was used as a ranking procedure in our analysis). The expression matrix (281 pathways × 888 drugs) of PRLs was obtained for further analysis. Next, the pathway-drug PRL matrix was converted to a subject-oriented matrix. Here, a rank-based pattern-matching Enrichment Score (ES) strategy that was based on the weighted Kolmogorov–Smirnov (KS) statistic in Drug Set Enrichment Analysis (DSEA) 26 was employed to perform the converted procedure. Given a PRL x and an RA subject y, the ES xy was calculated through DSEA approach. The KS-weighted ES could quantitatively measure the enrichment of signatures in the top/bottom ranked region. The ES value gives a range from 0 to 1. ES value tending to 1 indicates complete similarity and the value tending to 0 indicates the complete opposite. Finally, we built an ES matrix (888 drugs × 18 RA samples), the row corresponded to drug and the column represented RA subjects. By calculating the ES values, the drugs tending to mimic or suppress RA were quantified. Based on the drug-subject matrix, we sorted each row x according to the ES xy values of the drug x across the y = 1, …, Y RA subjects, and obtained a rank-based drug matrix R. Given an element R xy in R, it represented the rank of drug x according to its effect on RA subject y. In this case, the ES could sign whether a drug was a mimic or inhibitor in the development of RA. The significance of a drug for RA subjects was assessed by applying a nonparametric, rank-based procedure. For each disease subject y, the rank value of drug x was defined as: \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document} $${\rm{Ran}}{{\rm{k}}_{{\rm{xy}}}}{\rm{ = }}{{{\rm{sum}}\left({\left| {{\rm{E}}\left. {{{\rm{S}}_{{\rm{total}}}}} \right|} \right.{\rm{ \gt }}\left| {{\rm{E}}\left. {{{\rm{S}}_{{\rm{xy}}}}} \right|} \right.} \right)} \over {{\rm{leangt}}{{\rm{h}}_{{\rm{total}}}}}}$$ \end{document} The larger the rank value, the greater the likelihood of a suppressant; the smaller the rank value, the greater the likelihood of a good mimic. In the present study, the top 1% drugs with larger rank value were predicted as therapeutic drugs, and the top 1% drugs with smaller rank value could be considered as good mimics.

In response to the high attrition rates in the traditional drug development process, drug repositioning which recaptures marketed drugs for new indications has attracted the attention from pharmaceutical companies and medical researchers. Ashburn and Thor indicated that drug repositioning might shorten the time of drug development from 10–17 years to 3–12 years. 27 Previous studies have proposed numerous methods to build predictive models and some have shown promising results. 10 , 12 , 28 , 29 These emerging technologies enable investigators to identify candidate drugs that will prevent this disease and its complications. For example, Zhang et al 30 presented a drug repositioning strategy based on “omics” data mining to screen candidates for new indications in diabetes treatment, and successfully identified 9 drugs that might have the potential to treat diabetes. Existing drug repositioning methods mainly focus on discovering candidate drugs for a kind of disease, and are not suitable for predicting candidate drugs for an individual sample. The drug response heterogeneity makes the raise of new strategies that target genotypically well-characterized subpopulations of patients, instead of targeting huge populations. This conversion drives the researchers focus to settle personalized pathway profile by which the disease works and how to intervene on it. In this study, we proposed a computational method to predict candidate drugs from CMap database for RA, in a personalized way, contributing to revealing the molecular mechanisms and the future application of custom therapeutic decisions. In the present study, we first identify the pathway aberrance of individual RA subjects using iPAS algorithm, and a total of 4 differential pathways were identified. It is well known that RA is a kind of autoimmune disorder. Our study identified 4 differential pathways (cardiac muscle contraction, amoebiasis, amino sugar and nucleotide sugar metabolism, and protein processing in endoplasmic reticulum) in RA. Among them, cardiac muscle contraction pathway showed the highest altered frequency, which altered in 14 of 18 RA subjects. Previous studies revealed that RA patients had an increased risk of premature death compared with the general population, mainly due to cardiovascular disease, and the two diseases shared genetic and environmental risk factors. 31 , 32 Our result that identified cardiac muscle contraction as a differential pathway might affirm this declaration. Although, a rapid evolution of the care and treatment of patients with RA has contributed to a decrease in disease and disease-related complications. Efficient and safe drugs are still pressing problem. By drug repositioning, several small molecules that might suppress RA were identified, which might be critically important for disease prevention and early treatment. Among these candidate drugs, celecoxib had a current drug indication of RA, which is a COX-2 selective nonsteroidal anti-inflammatory drug used to treat various forms of arthritis. 33 The candidate parthenolide showed a current drug indication of anti-inflammatory. López-Franco et al 34 indicated that parthenolide could modulate the NF-κB-mediated inflammatory responses during vascular damage. Moreover, parthenolide has a variety of reported in vitro biological activities, such as blocking lipopolysaccharide-induced osteolysis, 35 inducing apoptosis of leukemia cells, 36 and activity against a parasite Leishmania amazonensis . 37 The candidate drug harmine was reported to promote differentiation of osteoblasts 38 and chondrocytes, 39 and inhibit osteoclastogenesis. 40 Suloctidil and prenylamine both showed a current drug indication of vasodilator. Unfortunately, they were withdrawn from market due to liver toxicity and cardiac arrhythmias, respectively. The other candidate drugs predicted by our method had current indications to treat other diseases and might complement the drug group for RA treatment. Given these screened potential drug for RA treatment, previous literature have reported that some of these drugs could exert certain relieved effect on RA. For instance, Danazol has been identified that treatment of it attenuates refractory autoimmune thrombocytopenia in rheumatic diseases successfully. 41 , 42 Some studies have demonstrated that piperlongumine relieves RA through expansion of myeloid-derived suppressor cells (MDSCs) and the inhibition of the Th17 response and activation of fibroblast-like synoviocytes (FLS) 43 represses dendritic cell maturation by decreasing production of reactive oxygen species, 44 and suppresses proliferation, migration and invasion of FLS 45 via animal model and patient samples. Celecoxib has been approved to treat RA and osteoarthritis (OA); 46 , 47 it, combined with Cilostazol, impedes proinflammatory factors in FLS of RA patients. 48 These results may provide more favorable evidence for the clinical use of these drugs. Meanwhile, we predicted nine mimics based on the similarity of pathway aberrances between disease and drugs. In this study, the drugs with higher ES values showed more similar to disease in the matter of pathway aberrances. Theoretically, it happened that the drugs with high ES values might mimic the disease at the molecular level. While it may not make sense to consider that the pathway profiles of one drug is positively related to those of one disease from the biological view. Thus, this study focused on the candidate drugs for RA treatment in a personalized way.

In the present study, we predicted several candidate drugs for RA treatment, and conducted drug discovery in a personalized way contributing to the future application of custom therapeutic decisions.