VIGLA-M: visual gene expression data analytics

Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 https://doi.org/10.1186/s12859-019-2695-7 SOFTWARE Open Access VIGLA-M: visual gene expression data analytics Ismael Navas-Delgado1*, José García-Nieto1, Esteban López-Camacho1, Maciej Rybinski1, Rocio Lavado2, Miguel Ángel Berciano Guerrero2and José F. Aldana-Montes1 From The 2017 Network Tools and Applications in Biology (NETTAB) Workshop Palermo, Italy. 16-18 October 2017 Abstract Background: The analysis of gene expression levels is used in many clinical studies to know how patients evolve or to find new genetic biomarkers that could help in clinical decision making. However, the techniques and software available for these analyses are not intended for physicians, but for geneticists. However, enabling physicians to make initial discoveries on these data would benefit in the clinical assay development. Results: Melanoma is a highly immunogenic tumor. Therefore, in recent years physicians have incorporated immune system altering drugs into their therapeutic arsenal against this disease, revolutionizing the treatment of patients with an advanced stage of the cancer. This has led us to explore and deepen our knowledge of the immunology surrounding melanoma, in order to optimize the approach. Within this project we have developed a database for collecting relevant clinical information for melanoma patients, including the storage of patient gene expression levels obtained from the NanoString platform (several samples are taken from each patient). The Immune Profiling Panel is used in this case. This database is being exploited through the analysis of the different expression profiles of the patients. This analysis is being done with Python, and a parallel version of the algorithms is available with Apache Spark to provide scalability as needed. Conclusions: VIGLA-M, the visual analysis tool for gene expression levels in melanoma patients is available at http://khaos.uma.es/melanoma/. The platform with real clinical data can be accessed with a demo user account, physician, using password physician_test_7634 (if you encounter any problems, contact us at this email address: mailto: [email protected]). The initial results of the analysis of gene expression levels using these tools are providing first insights into the patients’ evolution. These results are promising, but larger scale tests must be developed once new patients have been sequenced, to discover new genetic biomarkers. Keywords: Gene Expression Level Analysis, Nanostring Immune Profiling Panel, Metastatic Melanoma Background Theanalysisofthegeneexpressionchangesandtheircorrelation with clinical changes can be developed using the tools provided by those companies that currently provide gene expression platforms, such as NanoString. The analysis of gene expression levels is a very important problem when studying how patients evolve, and looking for a way *Correspondence: [email protected] 1Khaos Research, Universidad de Málaga, Málaga, Spain, Arquitecto Francisco Peñalosa 18, 29071 Málaga, Spain Full list of author information is available at the end of the article to find new genetic biomarkers helping in clinical decisionmaking.Currentsolutionsarebasedondesktoptools designed for use by geneticists. However, it is the clinicians who are in direct contact with the patients, so their ability to access such data in prior analyses could prove highly beneficial. NanoStringDiff [1] is an R package for differential expression analysis of NanoString nCounter data, whose main function for differential analysis is glm.LRT.The use of this tool requires knowledge of R programming language. © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 2 of 11 NanoStringNorm [2] provides R functions for preprocessing methods proposed by NanoString in addition to an extensible plug and play framework to incorporate other methods. This tool also requires R knowledge from its users. nCounter[3]istheadvancedanalyticaltoolprovidedby NanoString Technologies, Inc. This tool provides different analytical tools for gene expression data. It gives high quality results, but it has been designed for use by geneticists and so, it has a steep learning curve for physicians and clinical researchers. These users hope to make initial discoveries based on the experimental results, which normally requires help of an experienced nCounter user. nSolver [4] is the simplest analytical tool provided by NanoString Technologies, Inc. This tool provides a desktop user interface where users can import reporter code count (RCC) files for different analyses. Each RCC file contains results of a single flow cell of the cartridge. Each RCC file contains the number of counts for each target mRNA molecule (measured within one sample). A separateRCCfileiscreatedforeachflowcell(orsample)inthe cartridge. Although easier to use than nCounter, the tool still requires some learning effort, even for simple data analyses. In addition, the lack of some parametrization options makes it difficult to extract relevant information. As an example, a simple heat map analysis does not provide a way of filtering out those genes with lower variability. NanoStriDE [5] allowed biologists to take raw data produced by a NanoString nCounter analysis system and easily interpret differential expression analysis of this data represented through a heatmap. However, the system was not easy to use, and in some cases clinicians required biostatisticians to help them with these tasks. The analysis of gene expression levels is a complex problem that is being studied from different perspectives. Thus, [6] presents an approach to parallelize the TriCluster algorithm to be used in this context. In [7], another parallelization approach is proposed for efficient large-scale transcriptome data analysis. The gene expression platform used in this work is NanoString, which provides nCounter. However, nCounter has some limitations as the type of analysis is restricted to a predefined set, and the introduction of clinical features is a complex task. This project presents an approach to collect the clinical information using a structured database and a web user interface to introduce this information, including the results of the gene expression measurements, as a move a step closer towards the clinicians, in comparison to the nCounter tool. This paper presents VIGLA-M, an analytics tool to provide clinical researchers with early insights into the gene expression samples of a patient (discovering changes in their gene expression levels) or sets of patients (discovering patient clusters with a possible clinical correlation). Once the physicians or clinical researchers find relevant information in these data, geneticists can build on the preliminary analysis using, for instance, nCounter. Thetoolproposedinthispaperhasbeentestedona small set of patients of metastatic melanoma patients, but it is intended to support similar cases in other diseases where gene expression variation analysis would help in discovering new biomarkers. At present, immunotherapy for metastatic melanoma is based on stimulating an individual’s own immune system, inmostcasesthroughtheuseofspecificmonoclonalantibodies. The use of immunotherapy has meant that many patients with melanoma have survived and, therefore, it constitutes a promising present and future treatment in this field [8–12]. At the same time, drugs have been developed targeting specific mutations, i.e. BRAF (BRAF is a human gene that encodes a protein called B-Raf), resulting in significant responses in tumor regression (set up in this clinical study for 18 months), as well as a higher percentage of long-term survivors [13–16]. These significant responses are triggered by a massive process of antigens creation, that is needed to reignite a state of immunocompetence in the patient. It also leads to the activation of biomarkers that could guide response prediction. Unfortunately, the recent interest in developing drugs against melanoma hasnotbeen accompaniedbytheidentification of biomarkers that are predictive of response or toxicity to them. Although preliminary data determined immunohistochemistry was useful in the tumor expression of the PDL-1 protein, subsequent work has failed to validate its predictive usefulness. In fact, one of the main problems in the study of the immune system is its plasticityandadaptivenature.Tumorinfiltratinglymphocytes (TIL), another possible postulated and non-predictive marker, validated that PDL-1 expression can fluctuate within shorttime frames anditsexpressionchanges when biopsies are performed serially over time on the same tumor [17]. Therefore, more biomarkers are needed to optimize the sequences of treatments against melanoma. For these reasons, our work aims to combine both concepts: the use of therapy directed against the BRAF mutation and the reactivation of acquired immunity against melanoma. Therefore, it would be useful to look for these changes in blood biomarkers in relation to immunological mechanisms to identify predictive markers of response to new treatments, in order to identify long-term survivors with targeted therapy. Unfortunately, many patients will not benefitfrom these treatments,as they are difficult tofund in terms of human and economic expenditure. However, this project seems to be an essential first step and the finding of these biomarkers is yet to be explored. Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 3 of 11 Additionally, we have pursued a first approach to parallelize the algorithms using Big Data technologies (Apache Spark), to support large scale clinical assays. At this stage, VIGLA-M has been tested for a clinical assay with metastatic melanoma patients. Implementation In this section we present the three main functionalities that have been developed in the VIGLA-M project: clinical data collection, gene expression data collection and gene expression level data analysis and visualization. To collect the required data, we have developed a structureddatabase(using Oracle11g ExpressEdition)tostore clinical and gene expression data. Figure 1shows part of thedatabaseschema,whichisdesignedtostore:aminimal set of personal patient data (gender, date of birth, medical records, medication, diagnosis date, metastasis diagnosis date, etc.); the response to the treatment; the cancer’s progress, including the different visits to the hospital and thebloodtestvaluesateachvisit;theBRAFmutation;primary location of the cancer, and the files associated with the gene expression analysis. The clinical information is collected by the patient’s physician (previously registered and accepted in the clinical assay by the clinical project coordinator). Thus, the effort of collecting the data is shared by all the physicians participating in the project (clinical assay). The files obtained from the gene expression assays are loaded by the biologists (previously registered in the system). These files are obtained from the NanoString platform. The NanoString platform can analyze 12 samples in each cartridge. As a result, we obtain 12 RCC files with the gene counts for each of the gene panels. Each file isassociatedwiththegeneexpressionofapatientataspecific stage of treatment. These files are locally stored with a patient code, in addition to being annotated with the patient code from the project database, the sample collection date and the Nanostring experiment date. Then, the biologist uploads the files attached to the patient profile, along with other metadata. We have used the Immune Profiling Panel NanostringTM (770 genes), as it has been specifically designed for cancer projects where immune aspects are studied. This panel includes 24 different immune cell types, common checkpoint inhibitors, CT antigens, and genes covering both, adaptive and innate, immune responses. Data normalization Using the gene expression files, a set of analytic functionalities has been developed to discover patterns in the change of the gene expression levels. However, these files need to be preprocessed, as NanoString returns the flat countsofthegeneexpressionlevels.Thepre-processingis done to normalize the counts according to a set of control measures. The preprocessing is described in NanoString guidelines[18],anditcanbesummarizedwiththefollowing steps: Fig. 1 Extract of the database. This figure has been produced using SQL Developer Data Modeler (Oracle Database Design Tool). Users are connected with one or more profiles (user types), meanwhile profiles can also be connected with several users. Samples cannot exist without a user, and a user can have several samples. Patients are linked to a user with a profile of type “physician” Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 4 of 11 •Step 1. Generation of quality control flags (binding density, positive control linearity, limit of detection). •Step 2. Background correction using the negative control samples. •Step 3. Calculation of lane-specific scaling factors, based on the housekeeping gene set. •Step 4. Adjusting the flat counts with the lane-specific scaling factor. Step 1 is to ensure that the end user is informed about thequalityofthesamples.Thebinding densityisincluded in the NanoString output, therefore it is simply read from the sample file. Positive control linearity is the process of evaluating whether the flat counts obtained for synthetic-positive-control samples maintain the expected linear relationship resulting from their known amounts. Therefore, linearity is calculated as the r2value between the NanoString values and known values. The limit of detection is evaluated by comparing one of the positive controls and its distance from the background. The flags are directly passed on to the visualization module, so that theendusergetsnotifiedofanypossiblequalityproblems. Step 2 is done through background thresholding – setting the values of flat counts that fall below a certain threshold value to this threshold value. The threshold valueiscalculatedasatwostandarddeviationsofthenegative control sample values over the geometric mean of negative control sample values. Step 3 is carried out along the same lines as those presentedin[19].Asubsetofatleast3most stably-expressed housekeepers is found and then used to calculate the scaling factors. Step 4 is done simply by multiplying the flat counts with the scaling factors calculated in Step 3. The method for finding the most stable genes is illustrated with the following pseudocode snippet: Code Snippet 1 Skeleton of a jMetalSP application //expression_values is a genes x samples matrix; //number_of_genes_to_find is an integer findStableGenes(expression_values, number_of_genes_to_find): totalGenes=expression_values.rows V=matrix(size:totalGenes,totalGenes) M=vector(size:totalGenes) for i in [0 to totalGenes]: for j in [0 to totalGenes]: if i!=j: val=[] for column in expression_values: val.append(log2(colVect[i]/colVect[j])) V[i,j]=std_dev(val) M[i]=np.mean(V[i,:]) return get_indexes_of_lowest_K_values(M, K=number_of_genes_to_find) Subsequently, the lane-specific normalization factors are calculated as the inverse of the geometric mean of the in-lane housekeeper counts multiplied by the average value of these geometric means (across all lanes). Gene filtering The Immune Profiling Panel NanostringTM contains 770 genes,so anyvisualizationincludingallthese geneswould be difficult to explore for the end-users. For this reason, thedevelopedtoolincludesafilteringmoduleabletofilter out those genes that do not change (or change less than the rest) over the analyzed samples. The filter is also based on the geNorm method. Specifically, geNorm postulates finding n(>3) most stable housekeeping genes. For filtering, we use the same model to find k(parameter specified by the user, e.g. 100) least stable genes, i.e., the ones that display most variation (see thecodesnippet). Results To enable users to exploit the processed data we have developed VIGLA-M, a web service for accessing and analyzing these data. Registered users can log on to the tool and browse through the data of their patients. The tool available at http://khaos.uma.es/ melanoma has a demo user with demo patients with real gene-expression data. User and password for demo access are: physician and physician_test_7634 respectively. Physicians can only access their patients’ data, biologists can only access gene expression data that they have uploaded, and clinical assay administrators can access all the data (see Fig. 2for the demo user). Users can select different gene expression samples from one or more patients to start the gene expression analysis (see Fig. 3for the demo user). Therefore, using the normalized gene expression values of the most variable genes, the proposed tool is designed to perform two main functionalities for data analysis and visualization: • Heat Maps . A heat map provides a graphical representation of data stored as a matrix. The values included in the matrix are represented as colors. The tool provides a sliding element to set a parameter for extracting only the most variable gene expression levels (as a percentage of the total number of genes). It is worth mentioning that gene expression levels are standardized by calculating z-scores, hence providing balanced values before the visualization step. In addition to this, hierarchical clustering is computed to generate groups of both, genes and experiments. This will provide additional information to physicians about how samples can be organized. Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 5 of 11 Fig. 2 List of patients for a demo user. For each patient the user can click on Sample to add one or more samples for them • Inference of Gene Regulatory Network . Gene Regulatory Networks (GRNs) can be reconstructed from a time series of gene expression data. This information is of particular clinical interest because it can serve to detect genes that are interacting and so, detect new pharmacological targets. The tool provides a sliding element to set the parameter for representing only the most relevant ( n ) links of the network (as large networks will be difficult to interpret). This filter keeps the nodes involved in the n relationships with highest weights (see Fig. 4). To generate the network, two well-known algorithms [20] can be used. Both implementations are provided in the arboretum Python package (https://arboreto. readthedocs.io/en/latest/, namely: GRNBoost2,a fast algorithm based on stochastic Gradient Boosting Machine regression with early-stopping regularization; and GENIE3, the popular classic GRN inference algorithm based on Random Forest (RF) or ExtraTrees (ET) regression. This choice option will enable physicians to inspect different learning models for inferring the underlying GRNs. Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 6 of 11 Fig. 3 Sample selection from the available patient data. The user can select one or more samples from the same patient Data analysis through heat map visualization Theuseof heatmaps in VIGLA-Mprovides a flexible way of producing different types of analysis with respect to the selected gene expression data and their order. Thus, depending on these features, it is possible to produce four different representations to help clinicians analyze the expression results: •Type I. Comparison of several samples for a single patient (patient evolution). The visualization of the expression change in a patient enables the individual analysis of the patient data with a high degree of detail. •Type II. Comparison of the basal status of the patients (classification of the patients depending on the initial expression profile). This visualization enables the analysis of the gene expression profiles to allow users to analyze those patients with a specific pattern. For example, a user could select those patients with certain genes with a lower expression than the rest of the patients, to later explore the possible expression growth. Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 7 of 11 Fig. 4 Selection slide for controlling the size of the gene regulatory network. Once the user has select one or more samples from one or more patients, the user can select between the different analytical tools. The figure shows the two current options: Heat Maps and Gene Regulatory Networks. The slides enable users to determine which section of the data will be included in the graphical representation of the results •Type III. Comparison of the second sample of the patients (classification of the patients according to their expression after the evolution). This visualization also allows comparisons of patients after treatment. •Type IV. Comparison of the patients including both samples (profiling on the patient evolution in the gene expression). This visualization type is the most complete one provided by this tool. It enables the classification and visualization of differences between patients in the evolution of the gene expression. In this case, the patient samples are ordered showing first the initial sample of all the patients, and then, the results for the second sample. Figure 5shows a clinical case for Type I. This case can be replicated using the samples from Demo Patient 1. This patient suffered an allergic reaction to the treatment for the melanoma cancer. The first graphical representation is a heat map with hierarchical clustering according to both, genes and assays. As input this heat map only needs a set of ordered gene expression files. The heat map shows the most significant changes that have occurred in the expression of some genes of immune function, at different times of the treatment for metastatic melanoma. Prior to treatment, there was an activation of several genes, whose activity may be associated with increased tumor replication. Thus, there are complement (C3AR1), tyrosine kinase (FLT3 or CD135, therapeutic target of drugs such as sorafenib or sunitinib), interleukins (IL1R2) and adhesion molecules (CEACAM6) receptors or fibronectin (FN1), which are usually activated in many tumors. With BRAF inhibitors, these genes reduce their expression, showing the efficacy of this treatment in the subject under study. With the allergic reaction to the BRAF inhibitor other genes were activated, which were initially not active, such as IL3RA, CCL2, HLA-DRB4, IL8, FOS, PTGS2, CXCL1, and which were subsequently corrected with the use of drug sensitization with corticosteroids and gradual exposure to the drugs. The activation of these genes can improve the immune response, through the interleukins(IL3andIL8)orthebetterpresentationofantigens (HLADRB4), and so is able to produce a good response to the drugs. Finally, the heat map also shows how, over time, there are several genes that are activated (LILRA4, IDO1, MARCO, HLA-DQA1, MSR1, C1OB, IF27) and others that are deactivated (EGR2 and 1, IL32, EBP, TREM1, F2RL1, LCN2, CREB5, SLC11A1, IL1TRAP). Some of these genes are the key to immunotherapy resistance based on inhibitors of immunocheckpoints, so they can also be key in the immune response induced by targeted therapy against BRAF. Thus, IDO not only appears to be a key factor in the gene interactions tree (see Fig. 6), but it is a novel target. The inhibition of this enzyme is promising given the results of early phase drug development [21]. Data analysis through gene regulatory network reconstruction Inferring the interaction network of genes is a fundamental step towards understanding how a cell or an organism will respond to its environment. Its use in clinical assays Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 8 of 11 Fig. 5 Heatmap with hierarchical clustering generated for patient with allergic reaction. The X axis shows the four different samples for this patient over a stretch of time (t1: sample before the treatment; t2: sample after the treatment with BRAF inhibitors; t3: sample just after the allergic reaction; t4: sample after the allergic reaction disappears). The Y axis show the genes that had more changes in this time series Navas-Delgado et al. BMC Bioinformatics 2019, 20(Suppl 4):150 Page 9 of 11 Fig. 6 Gene regulatory network obtained from GRNBoost2 algorithm for patient allergic reaction. Thicker part of the lines shows the direction of the arc (that is, the regulatory effect). For example, on the left we can see that FOS regulates LBP. The longer the line, the longer the regulatory strength would be of interest as it provides the insight on how different genes affect the expression of other related (or not) genes. Following the example of the patient with an allergic reaction, we have used the tool to reconstruct the Gene Regulatory Network. The resulting network is shown in Fig. 6, which visualizes how activated and deactivated genes (according to the previous analysis) take part in the regulatory graph with strong interactions between them. In this figure, the thicker part of the lines in edges show the direction of the arcs (in form of arrows), thus indicating the regulatory effect. This differentiates GRN from gene coexpression networks, in which the direction of edges is overlooked. In Fig. 6,itcanbeclearlyseenthatgeneslike MARCO, C1OB, IF27, and CREB5 act as transcription factors that play important role in the regulation phase of the network. Strong interactions are registered for the following interactions: IL8→MSR1, MSR1→MARCO, MARCO→CREB5, and IF27→CREB5. This could shed light on the inmunotherapy response of the organism, so the clinician is now supported with additional arguments in adjusting the treatment process.