Personalized analysis of pathway aberrance induced by sevoflurane and propofol

Anesthetic agents are used in surgical operations to reversibly reduce consciousness and pain. Sevoflurane is an inhalational anesthetic. Propofol is a short-acting intravenous general anesthetic. The mechanism of anesthetic agents at pathway level on individual patients has not been reported to date. In the present study, pathway aberrance in the human atrial tissue in response to anesthetics was examined. Microarray data of anesthesia-treated samples were downloaded from the Array Express database. Pathway information was obtained from the Reactome Pathway Database. The individual pathway aberrance score (iPAS) was introduced to identify dysregulated pathways in individual patients. The present data demonstrated 157 dysregulated pathways in the sevoflurane group, and 44 pathways were identified with the least P-values. A subset of 49 differentially expressed genes (DEGs) that were shared between the expression profiling results and the dysregulated pathways results were constructed into a co-expression network. The top 5 ranked DEGs, nuclear receptor subfamily 4 group A member 3 (NR4A3), JUNB proto-oncogene, MYC proto-oncogene, tachykinin precursor 1 and nicotinamide phosphoribosyltransferase, were identified as important in the topology analysis. In the propofol group, 87 dysregulated pathways were identified and 44 pathways had the least P-values. In total 28 DEGs were constructed into a co-expression network, of which 5 DEGs were important in the topology analysis, NR4A3, suppressor of cytokine signaling 3, cyclin dependent kinase inhibitor 1A, C-C motif chemokine ligand 2 and C-X-C motif chemokine ligand 1. A total of 72 dysregulated pathways were identified in common in the two groups. In conclusion, the two types of anesthetics induced partially similar mechanisms. The pathways enriched by DEGs, particularly those that were unique to sevoflurane and propofol, may affect surgical outcomes and aid the prevention of complications from anesthetics.


Introduction
Anesthetic reagents are used in surgical operations to reduce consciousness and pain reversibly. They have been used for ~150 years due to their efficacy and safety (1). Sevoflurane is an inhalational anesthetic used for general anesthesia. It dissolves poorly in blood, and can be breathed out quickly by the lungs, providing rapid recovery. Propofol is a short-acting, most commonly used, intravenous general anesthetic and it is administered as an alkylphenol formulated in a lipid emulsion (2). It dissolves from the blood to the surrounding tissues and central nerve systems quickly. In addition, propofol exhibits rapid clearance from the body (3). Since these two types of anesthetic are widely used in surgery, a better understanding of their effect in biological processes is warranted.
Although the study of Lucchinetti et al (4) provided insights into the gene expression effects mediated by anesthetics administration inhuman atrial tissues, pathway aberrance in individual patients remains unknown.
Gene-based pathway analysis is becoming an important method for understanding disease mechanisms. Especially pathway analysis based on differentially expressed genes (DEGs) has become an important approach to study pathogenesis mechanisms in an individual sample (5). The analysis of altered pathways in an individual patient may help to understand the disease status and suggest customized therapies. Although the personalized interpretation of pathways can be demanding, most current pathway analyses have been developed to investigate deregulated pathways between two phenotype groups. Three methods were described by Khatri et al (6), overrepresentation analysis, functional class scoring and pathway topology-based approach. They were mainly focused on pathway aberrance between normal samples and cancer samples, but not suitable for an individual sample. The method of individual pathway aberrance score (iPAS) proposed by Ahn et al (5), which is based on overrepresentation analysis and functional class scoring, provides a series of analyses steps, containing four parts: Data processing, gene-level statistics, iPAS and a significance test, with the clinical importance of this method being to provide pathway interpretation of a single cancer sample.
With the development of molecular biology and bioinformatics, network analysis has become a more informative and powerful approach to study disease mechanisms (7). The relationship of disease to a specific network, especially a co-expression network, provides valuable information to analyze gene regulation (8). The topological characteristics may reflect the properties of the network and may be used to investigate the underlying implications of the network.
To test the hypothesis that sevoflurane and propofol may induce pathway aberrance, microarray data of human atrial tissue treated with sevoflurane and propofol were examined separately. DEGs and individual pathway aberrance were identified using individual pathway analysis and co-expression network analysis.

Datasets
Gene expression data. Microarray data of E-GEOD-4386 (4), and its annotation file were downloaded from the Array Express database (http://www.ebi.ac.uk/arrayexpress/). The dataset includes 10 sevoflurane-treated atrial tissues, 10 propofol-treated atrial tissues and 10 matched non-anesthetic-treated atrial tissues. The tissues were obtained from patients who underwent coronary artery bypass graft (CABG) surgery and were given the corresponding anesthetics. The sample platform that was used in the study was the A-AFFY-44-Affymetrix Gene Chip Human Genome U133 Plus 2.0 [HG-U133_Plus_2]. According to platform annotation files, probes were mapped to gene symbols by filtering with the function of the feature Filter (9). The model description was as follows: We used (w, b) as the coefficients and intercept of L1 logistic regression. (1) where L was the loss function, λ was a regularization parameter which had the ability to dispose high dimensional data, xn ϵRN was the n-th feature and yn ϵ RN was the label of the n-th sample. After performing ID transformation with the get SYMBOL function (10), 20,102 genes were obtained.
Pathway data. Information from gene sets representing biological pathways was obtained from the Reactome Pathway Database (http://www.reactome.org) (11). Pathways with large number of genes are more difficult to understand for human experts. Therefore, pathways with gene set size >100 were filtered out (5). Considering that information from different platform may vary, we filtered out the pathways that had no transformational gene symbols. Eventually 1,005 pathways were obtained, including 4,852 genes.
DEG analysis DEG screening. At present, many methods are routinely used to screen for DEGs, including Linear Models for Microarray Data (Limma) (12), significant analysis of microarrays (13), statistical methods for ranking (14), and R package siggenes (15) In the present study, the Limma package was selected to screen for DEGs. With this method, we obtained a P-value for all the genes. The genes with log |fold change (FC)|≥2 and P≤0.01 were selected for further analysis. In total, 269 DEGs were obtained from the sevoflurane samples, and 129 DEGs were obtained from the propofol samples.

Co-expression gene analysis.
To explore the interactive relationship between genes under different conditions, differentially co-expressed genes were constructed using the DEGs, with the following functions in the weighted gene co-expression network analysis (WGCNA) package (16).
where sij is the absolute value of the correlation coefficient between the profiles of nodes i and j, aij encodes the network connection strength between nodes i and j, and τ is the 'hard' threshold parameter. A total of 180 co-expression relationships (813 gene pairs) were identified when samples were treated with sevoflurane, and 119 co-expression relationships (1,216 gene pairs) with propofol. Then a network was constructed by linking atrial genes with the selected gene signatures using Cytoscape (17), a free software package for visualizing, modeling and analyzing the integration of biomolecular interaction networks with high-throughput expression data and other molecular states for further analysis in the topological characteristics.
Pathway analysis. Pathway analysis has become important in capturing clinical information. An iPAS model was introduced to identify altered pathways in an individual sample (18). Normal tissue data were represented as nRef. The process of this analysis was described as follows. Gene level statistics. Gene level statistics was described by Irizarry et al (19,20) by making use of multichip averages. Genes from normal samples were normalized by the function: (5) where represents the unit diagonal, and qk=(qk1,...,qkn ) for k=1,..., p represents the vector of the k-th quantiles for all n arrays For normal controls, the mean and standard deviation of gene expression level were calculated. For individual tumor cases, quantile normalization was performed after combining the single tumor microarray with all nRef samples. The gene-level statistics of individual tumor samples were standardized using the mean and standard deviation of the reference. The formula was: (6) where Zi indicates the standardized expression value of i-th gene and n represents the number of genes belonging to the pathway (5).
Pathway statistics. For each pathway, the pathway statistics is the average of all the gene level statistics which was extracted from the same pathway, calculated as: (7) Where zi symbolizes the standardized expression value of i-th gene belonging to the pathway, and n symbolizes the number of genes belonging to the pathway. The clinical characteristics of a single anesthesia sample can be obtained by measuring its difference against common characteristic of normal samples, which was represented as nRef in our study.
Dysregulated pathways. Significance test was performed to assess altered pathways associated with anesthetics administration. To assess the statistical significance, pathway statistics, were tested by the two-group comparison method Wilcoxon (21), with the function of: where n is the number of samples. The significance levels were revised by false discovery rate (FDR) (22). Finally, each pathway yielded a P-value. A pathway with P<0.01 was considered to be significantly different between the anesthetics group and the normal group.
Individual altered pathway. Significance can be obtained against the null distribution generated from normal samples. The statistic of the null distribution was obtained by comparing every normal sample with all nRef samples. The P-value of each pathway was obtained by comparing each anesthetic sample with the null distribution. P<0.05 was considered to be significant difference.
Analysis of DEGs in dysregulated pathways. DEGs in dysregulated pathways were mapped to co-expression networks to analyze the co-expression relationships among these DEGs. In the sevoflurane group, 49 DEGs were constructed into a network. In the propofol group, 28 DEGs were constructed into a network. The procedure for analysis of topological characteristics of the networks was performed with the following steps.
The degree for counting the number of interactions of a given node was measured as described (23). Stress centrality (24) for the number of nodes in the shortest path between two other nodes was calculated as follows: (11) Closeness centrality (24) for measuring the average length of the shortest paths to access all other proteins in the network was calculated with the function of (12) Betweenness centrality (24) for determining how the neighbors of a node are inter-connected was calculated with the function of (13) Where σst is the total number of shortest paths from node s to node t and σst (v) is the number of those paths that pass through v. dG (s, t) represents the length of the shortest path between two nodes s and t in graph G, which is the sum of the weights of all edges on this shortest path.
Statistical analysis. All data were expressed as mean ± standard deviation. The Student's t-test was used for comparisons between groups unless otherwise noted. P<0.05 was considered to indicate a statistically significant difference.

DEG data analysis.
In the sevoflurane group, 269 DEGs were identified with logFC≥2 and P≤0.01, and they were mapped to 189 pathways. In the propofol group, 129 DEGs were identified, and they were overlapped to 140 pathways.
Dysregulated pathway analysis. Pathway statistics of anesthetic-treated samples and normal samples were obtained by using the method of iPAS. The pathway statistics of the pathways in which most DEGs were participating were displayed in Figs. 1 and 2 to illustrate the difference between the two groups. Fig. 1 presents pathway statistics for each of the samples and controls for the 'chemokine receptors bind chemokines' pathway, as this was the pathway that most DEGs were associated with in the sevoflurane group, the P-value for this pathway was P=0.002838. For the propofol group, the pathway that the most DEGs were associated with was the 'signaling by TGF-β receptor complex', and the P-value for this pathway was P=0.01213. However, although these two pathways had the highest number of DEGs that were associated with them, they did not exhibit the lowest P-values in their respective treatment groups. In total, 157 dysregulated pathways with P<0.01 were identified in the sevoflurane group. Following ranking of the pathways by P-value, the top 44 pathways with the lowest P-value (all P=0.000247) were listed in Table I. In the propofol group, 82 dysregulated pathways were obtained. Following ranking of the pathways by P-value, the top 29 pathways with the lowest P-value (all P=0.000375) were listed in Table II. Among these dysregulated pathways, 72 dysregulated pathways were common between the two groups.
The dysregulated pathways were converted into a heatmap to help understand pathway aberrance in different samples. Heatmaps are commonly used in pathway analysis and are constructed as a visual aid to analyze dysregulated pathways. The heatmap of the dysregulated pathways in the sevoflurane group is illustrated in Fig. 3 and the heatmap of the dysregulated pathways in the propofol group is illustrated in Fig. 4, and two

Analysis of DEGs in the dysregulated pathways.
To analyze the DEGs involved in the dysregulated pathways, the genes were mapped into the corresponding co-expression network, and then the interactive relationships and the topological characteristics of the DEGs were obtained. The interactive relationship networks are illustrated in Figs. 5 and 6 for the sevoflurane and the propofol groups, respectively.
For the sevoflurane group, 49 DEGs were constructed into a network (Fig. 5). In the network, only two genes were downregulated, Fanconi Anemia complementation group F and cingulin, but the two genes didn't relate with other DEGs (Fig. 5). All the other genes were upregulated at different extents. Amongst these, six genes, nuclear receptor subfamily 4 group A member 3 (NR4A3), MYC proto-oncogene (MYC), tachykinin precursor 1 (TAC1), JUNB proto-oncogene (JUNB), nicotinamide phosphoribosyltransferase (NAMPT) and interleukin 1 receptor antagonist (IL1RN), were central nodes of interactions, displaying the most complex relationships with the other DEGs.    The topology properties of all the DEGs in the network were analyzed, including betweenness centrality, closeness centrality, degree and stress. Following ranking DEGs by their topology characteristics, the top five ranked genes were NR4A3, JUNB, MYC, TAC1 and NAMPT (Table III).
For the propofol group, 28 DEGs were constructed into the network, all of which were upregulated (Fig. 6). The genes encoding interleukin (IL) 6 and IL8 were the most highly upregulated compared with the other DEGs. The top five ranked genes by topology characteristics were NR4A3, suppressor of cytokine signaling 3 (SOCS3), cyclin dependent kinase inhibitor 1A (CDKN1A), C-C motif chemokine ligand 2 (CCL2) and C-X-C motif chemokine ligand 1 (CXCL1) ( Table IV).

Discussion
In the present study, pathway aberrance induced by anesthetic treatment in patients who underwent cardiac surgery was analyzed. The microarray data used in the present study were based on a previously published study by Lucchinetti et al (4), and the aim was to identify altered pathways in individuals.
Sevoflurane and propofol were the experimental anesthetics involved in the data, as these are the most commonly used anesthetic agents in everyday clinical practice in several countries. Sevoflurane is a volatile anesthetic agent that is widely used by mask induction (25). Propofol is the intravenous anesthetic agent of choice in surgery, because of its favorable operating conditions and association with rapid recovery (26).
Pathway regulation patterns of samples treated with general anesthetic were distinct from those in the normal groups. With the iPAS method, anesthetic-treated samples were analyzed for pathway aberrance. A total of 157 dysregulated pathways were obtained from the sevoflurane samples. There were 44 pathways significantly influenced (all P= 0.000247) by sevoflurane treatment, such as the 'activation of C3 and C5' pathway, which participates in regulation of the immune system. Three pathways involved in vitamin metabolism were also significantly altered, which was consistent with the fact that vitamin C supplementation offers significant protection against anesthetic-induced neurotoxicity and behavioral alterations, as proposed by Xu et al (27). Seven pathways involving the Toll like receptor (TLR) were significantly altered, which is consistent with the report by Sun et al which demonstrated that sevoflurane exerts direct relaxant and anti-inflammatory effects by inhibiting the TLR4/ nuclear factor-κB pathway (28). The present results indicate that these identified pathways were closely related to sevoflurane treatment, and suggest that special attention should be given to patients with correlated diseases when administering sevoflurane anesthesia. In the propofol samples, there were 82 dysregulated pathways, of which 29 pathways were significantly influenced (all P=0.000375), such as the 'nicotinamide salvaging' pathway, which is involved in nicotinamide utilization. Propofol administration reduces the heart rate and the baseline cerebral blood flow. These 29 identified pathways in the present study were demonstrated to be closely related to propofol treatment. In the original study by Lucchinetti et al (4), the peroxisome proliferator-activated receptor γ coactivator 1-α pathway and the granulocyte-colony stimulator factor survival pathway, integral components of the protective program in the heart, were proposed as potential therapeutic targets in perioperative cardioprotection. A total of 72 pathways were common in both groups, indicating that, in some aspects, sevoflurane and propofol treatments caused similar effects. Previous studies have provided evidence that both sevoflurane and propofol enhance the function of the γ-aminobutyric acid type A receptor (GABA A R) in neurons and in recombinant α 1 β 2 γ 2 GABA A R systems (29,30). The difference between sevoflurane and propofol is that they bind to separate sites and modulate receptor function via distinctly different mechanisms (31).
To visualize dysregulated pathways and provide a sensitive measurement of a patient's clinical features, heatmaps of the two anesthetic groups were constructed. In the heatmap, intuitive information of the comparison of pathways between anesthetic-treated samples and normal samples was  obtained. The differences between the normal samples and anesthetic-treated samples are evident, and dysregulated pathways varied among the different anesthetic agents.
Networks were constructed for the DEGs involved in the dysregulated pathways. In the sevoflurane group, two genes were downregulated and 47 genes were upregulated. Six DEGs in the center of the interaction network were closely linked with other genes, so the corresponding pathways, nuclear receptor transcription and circadian clock pathways, may have been influenced by sevoflurane. In the propofol group, all the 28 DEGs were upregulated, and several genes in the center of the interaction network were significantly altered, including IL6 and IL8, which were also identified by Lucchinetti et al (4). As IL6 and IL8 participate in immune response, it is probable that administration of the anesthetic propofol may influence the immune system.
Following ranking genes by their topological characteristics, the top ranked gene in the sevoflurane group was NR4A3. In the interaction network, NR4A3 was a central node with a complex relationship with the other genes. NR4A3 is important in the regulation of cell proliferation, differentiation, metabolism (32) and apoptosis (33). It is involved in the 'nuclear receptor transcription' pathway. These genes may potentially be targets for therapies to prevent anesthetic-associated side effects in patients undergoing CABG surgery, and may contribute to the clinical action of sevoflurane and propofol. In the propofol network, the top ranked gene was SOCS3, which functions in the negative regulation of fetal liver hematopoiesis and placental development. SOCS3was involved in the pathway of 'growth hormone receptor signaling', 'interleukin-6 signaling', and 'regulation of interferon γ signaling'.
In summary, the present study identified 157 dysregulated pathways in the sevoflurane group and 87 dysregulated pathways in the propofol group. Of these, 72 dysregulated pathways were common in both groups. Following ranking DEGs by topology analysis, the top 5 ranked DEGs in the sevoflurane group were NR4A3, JUNB, MYC, TAC1 and NAMPT, and in the propofol group, NR4A3, SOCS3, CDKN1A, CCL2 and CXCL1. These genes and pathways may be important in regulating the effects of sevoflurane and propofol administration. Among them, several have been reported previously, while others were reported in the present study for the first time. In conclusion, both types of anesthetic caused gene and pathway aberrance, and the present results may provide indicators of physical condition in patients that may aid in preventing complications from anesthetics. Further prospective studies will be required to validate the present results.