extfieldsofstudy:
Medicine
Biology
Computer Science

text:
Prediction of Drosophila melanogaster gene function using Support Vector Machines

Background While the genomes of hundreds of organisms have been sequenced and good approaches exist for finding protein encoding genes, an important remaining challenge is predicting the functions of the large fraction of genes for which there is no annotation. Large gene expression datasets from microarray experiments already exist and many of these can be used to help assign potential functions to these genes. We have applied Support Vector Machines (SVM), a sigmoid fitting function and a stratified cross‐validation approach to analyze a large microarray experiment dataset from Drosophila melanogaster in order to predict possible functions for previously un‐annotated genes. A total of approximately 5043 different genes, or about one‐third of the predicted genes in the D. melanogaster genome, are represented in the dataset and 1854 (or 37%) of these genes are un‐annotated. Results 39 Gene Ontology Biological Process (GO‐BP) categories were found with precision value equal or larger than 0.75, when recall was fixed at the 0.4 level. For two of those categories, we have provided additional support for assigning given genes to the category by showing that the majority of transcripts for the genes belonging in a given category have a similar localization pattern during embryogenesis. Additionally, by assessing the predictions using a confidence score, we have been able to provide a putative GO‐BP term for 1422 previously un‐annotated genes or about 77% of the un‐annotated genes represented on the microarray and about 19% of all of the un‐annotated genes in the D. melanogaster genome. Conclusions Our study successfully employs a number of SVM classifiers, accompanied by detailed calibration and validation techniques, to generate a number of predictions for new annotations for D. melanogaster genes. The applied probabilistic analysis to SVM output improves the interpretability of the prediction results and the objectivity of the validation procedure.

http://www.biodatamining.org/content/6/1/8 and only 58% have a GO-BP (Biological Process) term [1]. If one excludes genes that only have an electronic annotation, then only 41% have a GO-BP term [1]. Functional genomic data and in particular microarray mRNA expression data have been used by numerous researchers as a means to help predict the function of un-annotated genes. The analysis of such data is based on the premise that genes involved in a particular biological, molecular, and/or biochemical process are often co-expressed. This co-expression is dependent on the presence of common cis-regulatory elements of the co-regulated genes that bind one or more common transcription factors. A common approach to examine co-expression profiles from microarray experiment data is to use clustering analysis. In this type of analysis, genes are organized and grouped based on their expression profile, with genes having similar expression pattern being "clustered" (grouped) together. The results of clustering analysis depend (a) on the metric used for comparing the expression profiles of genes that are analyzed and (b) on the clustering algorithm used. Euclidean distance and Pearson's correlation coefficient are two simple commonly used metrics [2]. Popular clustering algorithms are hierarchical clustering, k-means clustering and Self Organizing Maps [3].
Although clustering provides an insightful way of exploring gene co-expression patterns, it does not directly convey any information regarding potential functions of un-annotated genes. For that purpose, researchers have applied a further layer of analysis to the results of clustering, making use of additional data containing annotations of genes being analyzed. According to this type of analysis, clusters that are "enriched" for a particular gene function (i.e. they contain a larger than by chance number of genes having that function) are "assigned" that function. Subsequently, and following the "guiltby-association" principle, all un-annotated genes in that cluster are predicted to have the same function [4][5][6][7].
In addition to the clustering analyses of genes in an unsupervised fashion and without annotation information, a number of supervised methods have been proposed for predicting gene functions. These methods make use of available information about the annotations of genes in order to discover gene expression patterns that characterize those annotations. Then, functional predictions for un-annotated genes are made based on how well their expression profiles are matched with "annotated" patterns. Perhaps the most popular class of supervised methods are the binary classifiers, where first, expression data are separated into two classes either having or not having a particular annotation, and then un-annotated genes are predicted of having or not having the annotation, based on which class their expression profile falls under. Examples of applications of these methods include Support Vector Machines (SVM) [8,9], random forest [10], neural networks [11], factor analysis [12], logistic regression, linear discriminant and quadratic discriminant analysis [13]. The application of the aforementioned approaches for gene function prediction benefits from the use of large data sets where many different experimental treatments or conditions make up the microarray expression dataset (e.g. developmental time points, mutations, specific tissues, environmental conditions, drug treatments, etc.). For the most part, the studies that have tried to do large-scale gene prediction assignment have used well-known model organisms where large microarray datasets were available. These include studies in Saccharomyces cerevisiae [7,14], Arabidiopsis thaliana [13], D. melanogaster [15], and Mus musculus [9]. Other approaches include assembling and integrating a large number of datasets from different experimental approaches (e.g. microarray expression, genetic interaction, and protein-protein interaction) to create a http://www.biodatamining.org/content/6/1/8 network which in turn improves the robustness of the gene function predictions. These include studies in Saccharomyces cerevisiae [16], C. elegans [17], D. melanogaster [10,18] and Mus musculus [19].
In this study we focus on the use of SVM for the prediction of D. melanogaster gene functions. SVM is a popular machine learning method for classification and regression. Its proven high performance as well as its solid theoretical basis justify its frequent use in many fields, including bioinformatics and predictions of gene functions. As a two-class classification tool, SVM attempts to separate the data points not in the original feature space but in an "enlarged" higher-dimensional space instead. The seemingly highly computationally expensive data transformation is not performed but instead, ingeniously, the separation is performed "implicitly" based on their distances measured with the use of a kernel function. The SVM algorithm uses optimization techniques to find the surface that gives the optimal margin between the points of the two classes [20,21].
Despite its popularity in bioinformatics along with other research areas, SVM has been used for the prediction of GO-BP annotations for genes of various organisms (such as S. cerevisiae [8] and M. musculus [9]) but not for D. melanogaster. Our study fills this gap, investigating how well this very popular method works for D. melanogaster gene expression data, in particular when the dataset is of a specific structure imposed by the nature of time-course experiments, as the one we use in our study.
Using microarray data from the life cycle of D. melanogaster [22], and a controlled vocabulary for annotation of biological processes associated with D. melanogaster genes from the Gene Ontology Consortium (GO-BP) [23], in this study we propose a method of predicting gene function of un-annotated genes in the D. melanogaster genome by using Support Vector Machines and a two-level data splitting rotation scheme for validation (double cross-validation). Our prediction method was evaluated also externally with the use of an independent dataset.
Using this approach we have been able to provide a putative GO-BP term for about 77% of the un-annotated genes represented in the dataset and about 19% of all of the unannotated genes in the D. melanogaster genome helping to bridge the gap for the large number of genes that have little or no annotation. In addition, this SVM approach provides a precision and probability estimate that can help guide users as to the likelihood a given gene belongs to GO annotation class.

Microarray data and annotation sources
The microarray data used in this study was obtained from the series of 138 cDNA microarrays spanning the life cycle of Drosophila melanogaster [22]. The microarrays contained probes from 6765 cDNA clones and represented 5043 genes, roughly one third of the total number of genes in the D. melanogaster genome. From those genes, 1854 were not annotated with any GO-BP term at the time of the analysis. Data were obtained from the Stanford Microarray Database [24] and normalized using a ratiobased method according to the original publication. The dataset can also be obtained from the Gene Expression Omnibus [25], GEO accession number GSE4347. cDNA clone names were converted to primary Flybase Gene Identifiers (FBgn ids) from release 4.2 of the D. melanogaster genome using annotation available at Flybase [26]. Biological process annotation was downloaded from the Gene Ontology Consortium [23] in February http://www.biodatamining.org/content/6/1/8 2006. Only GO-BP categories containing a minimum of 10 and a maximum of 999 genes in the dataset were included, and 788 categories met these criteria. Clones with duplicate computer gene (CG) numbers were purposely not removed as we were interested in investigating the consistency of the predictions across the duplicates.

Support Vector Machines
A Support Vector Machine (SVM) is a classification and regression method originally developed by Vapnik [21]. Given a set of p-dimensional vector data x i and their labels y i taking the values {−1, +1}, a linear Support Vector Machine finds the optimal hyperplane that separates the "positive" from the "negative" class. This plane is maximizing the margin between the two classes. According to the mathematical formulation of the problem, the solution refers to a weight vector w and a scalar b that satisfy the optimization problem where N is the number of samples.
When classes are overlapping and misclassifications are allowed, the above constraints become where ξ i are slack variables and C a constant. Solutions can be generated by using the "Lagrange" formulation of the problem and its "Wolfe dual" problem. Given the solution, the weight vector w, for a new point x, the function f (x) = w T x + b is calculating the discriminant value for x, which can be used for its classification to the positive or negative class. It turns out that only a number of training points are important for the determination of the solution, the support vectors. This mechanism can be applied to problems of non-linear separation, after mapping the data onto a higher dimensional space, with the use of a mapping . is determined by a kernel K(., .) such that K(x, y) =< (x), (y) >. The optimization problem is similar with the linear separation case, and following the Lagrange formulation and given a solution of Lagrange multipliers α i , the discriminant function is given by . . , N s are the support vectors. It is important to note here that the solution and discriminant function depend on the data only through the kernel, and also that the explicit expression of the mapping is not needed for the solutions. This is one of the reasons that make SVM computationally attractive. More information on SVM can be found in [27].
It is evident that SVM depends on the choice of the kernel. After performing a small set of experiments with pilot runs evaluating a number of kernel choices, we decided to use a radial basis kernel, x T x, over a linear or polynomial kernel. Our decision was also supported by the large popularity of this kernel (previously used, for instance, in [8,9]). We used as σ a heuristic value equal to the median value of the distances between positive training points and their nearest negative training points. To combat asymmetric classes, a suggestion from [28] is implemented, where K(x, x) is augmented by λ · |C|/N http://www.biodatamining.org/content/6/1/8 where |C| is the size of the class C that x belongs to, and λ a tuning parameter. For our calculations we use λ = m/2, where m is the median value of the diagonal of the kernel matrix. For the implementation of the SVM algorithm we use the publicly available Gist package [29].

Estimation of class membership probabilities
For the translation of the discriminant values f (x) to posterior probabilities of class membership, p(y = 1|x), we use the method proposed in [30]. The author proposes a sigmoid function model, where with parameters a, b. The estimation of the parameters a, b is done by solving a maximum likelihood problem and , where N + , N − are the sizes of the positive and negative class respectively. For the sigmoid model we use the model fitting algorithm proposed by [31], which is implemented in the Gist package.

Cross-validation and performance evaluation
For each selected GO-BP category a set of annotated genes consisting from the positively annotated and a subset of the negatively annotated of size equal to up to four times the size of the positively annotated is participating in the evaluation of the SVM prediction algorithm. If I is the index set of all these genes participating in the cross-validation, a partition {I 1 , I 2 , I 3 , I 4 } is generated, where each I j has equal size and equal number of positively labelled genes (to the degree that this is achievable). Two of those sets are used as "training set" for the SVM, one set is used for calibration and sigmoid fitting, and the last set is used as "test" set, where probability estimate values are output from SVM and compared with the known annotations of the genes for the evaluation of the method. For each permutation (α, β, γ , δ) of {1, 2, 3, 4}, and if gene g i belongs to set I α , we denote with p β|γ ,δ (i) the probability estimate for gene g i from an SVM classifier that was trained with training set I γ ∪ I δ and calibrated with the set I β . Keeping the test set fixed we can have three different ways of constructing the training and calibration sets, and therefore each gene of the test set has three different probability estimate values. For example, for gene g i in I 1 , we have p 2|3,4 (i), p 4|2,3 (i) and p 3|4,2 (i). Those three values are averaged to produce a unique probability estimate for every gene. For each test set (or fold) the probability estimates and true annotations of the genes are used for the calculation of performance measures such as precision at 20, 30 and 40. The values of those measures are then averaged to produced one final measure for the whole SVM procedure. Note here http://www.biodatamining.org/content/6/1/8 the distinction between the probability estimate generation for each gene that involves averaging over the three possible arrangements of training and calibration sets, and the calculation of the performance measure that is done independently for each one of the four folds (test sets), using each time only the genes of the test set.
For the calculation of precision and recall and the generation of the precision-recall plots the R package ROCR [32] was used. For the precision-recall plots, the curves corresponding to the different folds were averaged vertically (i.e. precision values were averaged for the same recall value).

Gene function prediction using Support Vector Machines
The prediction of gene functions was performed with the use of a number of SVM classifiers. Annotation predictions were made independently for each of 788 GO-BP categories selected based on their size (see Microarray data and annotation sources in Methods section), based on the results of an SVM classifier. Each SVM was trained with the use of a training set of data points of known labels (in this case, genes with Gene Ontology Biological Processes (GO-BP) annotation). To remedy any model fitting problem caused by the high imbalance between the positive and negative classes, for every GO-BP category and SVM classifier we limited the negative examples by randomly selecting a subset of them with a ratio of positive to negative examples being 1:4. The trained system was used for the classification of new, unlabelled data (i.e. prediction of functions for un-annotated genes). The classification was based on the discriminant value of the data point, which is output by the SVM. As the discriminant value measures the "distance" of the data point from the hyperplane that separates the two classes, traditionally points with discriminant value larger than some threshold (typically zero) are classified as positives and the rest as negatives. Here, we refrain from translating the discriminant values to binary classifications, but, instead, we used a previously published algorithm to estimate the posterior probability of a gene membership to each GO-BP category. The algorithm fits a sigmoid function to the discriminant values of labeled data and it uses this function to approximate posterior class probabilities for unlabelled data [31]. A flowchart describing the probability estimation procedure using SVM and the sigmoid fitting function is shown in Figure 1.
We validated the performance of SVM and sigmoid fitting function with the use of stratified cross-validation, where the initial set of known data points (positively and selected negatively annotated genes) was randomly partitioned into four equal size sets, with equal number of positively annotated genes [33]. For every possible set combination, two fourths were used as "SVM-training" set (used to train the SVM classifier), one fourth as a "sigmoid-training" set (used for fitting the sigmoid function to the discriminant values) and one fourth as an overall "test" set, for which the membership probabilities were estimated. Predictions results for the latter set can then be validated, since the true labels (gene functions) are known. For more details in the cross-validation procedure used see Methods section.
For any chosen probability threshold value p t we can calculate the corresponding precision value, defined as the ratio of the number of genes that have average probability estimate equal or larger than p t and belong to the class, divided by the number of genes that have a probability estimate equal or larger than p t (i.e. True Positives / Predicted Positives). Similarly, recall is defined as the proportion of the genes annotated with the http://www.biodatamining.org/content/6/1/8 term that also have average probability estimate equal or larger than the threshold (i.e. True Positives / Total Positives). While precision gives a measure of the accuracy of the prediction, recall measures the coverage or completeness. The prediction performance of SVM for a specific GO-BP category was evaluated with the use of a measure equal to the precision value corresponding to a specific recall value (e.g. 20%, 30%, 40%). Precision at 20, 30 and 40 were calculated for each fold and averaged over all four folds, and the number of categories with high performance measure values were calculated. More specifically, the method identified 39 high-precision categories reaching precision at 40 values equal or larger than 0.75. Table 1 presents these categories. Some of these categories are related through ancestor-predecessor relationships and therefore have a certain amount of redundancy. (The relationship between the categories is shown in a GO graph presented in Additional file 1: Figure S1.) Taking into account these relationships as well as the size of the categories, we came up with a subset of 24 categories of minimum redundancy.
Two of the high-precision GO-BP categories, DNA-dependent DNA replication and oxidative phosphorylation, were chosen for further analysis. Precision-recall plots were generated for both categories and presented in Figures 2A and 2B.
Developmental gene expression profiles for the genes in these two categories are shown in Figure 3. For all of the expression profiles, the gene expression ratio for a given gene at a given time point is calculated by dividing the intensity signal for that gene at that time point by the average intensity for that gene throughout development. For the genes belonging to the DNA-dependent DNA replication category, there is a strong upregulation of these transcripts from the earliest time points in development (0.5-1.5h embryos) through approximately mid-embryogenesis (11-12h embryos) ( Figure 3A). By http://www.biodatamining.org/content/6/1/8  pupal stages (M00h-M48h) these genes are down-regulated and then are dramatically upregulated in late metamorphosis (M60h) and stay high until early adulthood ( Figure 3B). It should be noted that during embryogenesis and metamorphosis the organisms are not feeding or moving but during the larval and adult stages they are. Interestingly, the oxidative phorphorylation genes are up-regulated well in advance of the organisms moving or feeding, suggesting that the organism is preparing itself for the next stage of the development that will require the enzyme activities to carry out this process.

Gene function predictions for un-annotated genes
One approach of declaring an un-annotated gene as predicted to be annotated with some GO-BP term, is to ascertain if its probability estimate output from the SVM specific for http://www.biodatamining.org/content/6/1/8 this term was larger than some threshold value (e.g. 0.5) and then assign the un-annotated gene with that term. However, drawbacks to this approach include the arbitrary choice of the threshold and the fact that it ignores the evaluation of the performance of the classifier, as measured by cross-validation. In this approach, any gene that has probability estimate score larger than 0.5 (or any other pre-specified threshold value) for two categories A and B, will be predicted to belong to both A and B, regardless how well SVM is performing for these categories. A third problem with this method is that it assumes that probability estimate values are comparable between different GO-BP categories. http://www.biodatamining.org/content/6/1/8 Although the purpose of the sigmoid function fit method is to provide with comparable estimates, one cannot be certain that this is necessarily obtained. We have chosen an approach that takes into account the predictive behaviour of SVM for a specific GO-BP term, when a prediction is made for an un-annotated gene, as well as it addresses the other issues mentioned above. For each gene-GO category pair a "gene precision" score is assigned to be equal to the maximum precision that can be achieved using thresholds smaller or equal to the probability estimate of the gene (output from SVM). If the gene is to be included in the "positive" predictions from SVM (and therefore to be assigned the GO-BP function), the threshold has to be such that it will give a precision value for the classifier not larger than the "gene precision" score for the gene-GO category pair. This score measures how probable is for the gene to truly belong to GO-BP, given its probability estimate output from SVM, and also given the precision value of SVM for the particular category. This gene precision score is calculated separately for each one of the four folds (since for each fold probability estimates are generated for all un-annotated genes) and the four fold-specific scores are averaged.
We can then declare as gene function predictions all the gene-category pairs that have "gene precision" score equal or larger than some threshold value. More specifically, 9887 new predictions representing 1422 genes (CG numbers) were generated corresponding to a gene precision score larger than 0.75. A total of approximately 5043 different genes, or about one-third of the predicted genes in the D. melanogaster genome, are represented in the dataset. 1854 (or 37%) of these genes are un-annotated and therefore this method is providing gene function predictions for about 77% (1422/1854) of the un-annotated genes represented on the microarray and about 19% of all of the un-annotated genes in the D. melanogaster genome. A graphical representation of these results is shown in Figure 4. The complete list of the these predictions can be found in Additional file 2: Table S1.

Experimental support of predicted gene functions
The un-annotated genes contained in two of the high-precision GO-BP categories, DNAdependent DNA replication and oxidative phosphorylation are listed in Tables 2 and 3. For each gene, the gene precision score is indicated. The transcripts for genes involved in the same function are often co-localized within cells. Each of the genes (i.e. CG numbers) on each list was searched in a D. melanogaster FISH database [34,35] to visually examine the spatial and temporal mRNA expression for that gene during early embryogenesis. Figure 4 The Proportion of genes that are Annotated, Un-annotated, and Predicted by the SVM to Belong to a GO-BP Category. Genes from the developmental time course data set that had no GO-BP annotation were broken down into genes that were predicted with high confidence to belong to a GO-BP category (predicted), and those that had low prediction values (not predicted). http://www.biodatamining.org/content/6/1/8 0.798 0.566 * CG number not known DNA-dependent DNA replication GO-BP category (along with oxidative phosphorylation) was selected for further biological validation using fluorescent in situ hybridization (FISH) analysis of early D. melanogaster embryos. In this table, genes that are predicted to belong to this category are shown along with information regarding their presence in the FISH database. Each of the genes (i.e. CG numbers) in the list was searched in a D. melanogaster FISH database to visually examine the spatial and temporal mRNA expression for that gene during early embryogenesis. Genes that are marked with (+) or (−) sign had images present in the database. Those that had patterns that largely matched that of known genes in the category are marked with a (+). If either their temporal or their spatial pattern did not match the known gene pattern, they are marked with (−). For each gene the "gene-precision" score and the average probability estimate output from SVM are reported.

Table 2 Gene Lists and Fluorescent In Situ Hybridization (FISH) analysis for DNA-dependent DNA replication GO-BP category
Genes that are marked with a (+) or (−) had images present in the database. For the genes with images present in the FISH database, their spatial and temporal pattern were compared with those of known genes annotated with the category. Genes that had patterns that largely matched that of known genes in the category are marked with a (+). If either their temporal or their spatial pattern did not match the known gene pattern, they are marked with (−). http://www.biodatamining.org/content/6/1/8 One annotated gene from the dataset known to belong to the GO-BP and one unannotated gene from each of the two GO-BP categories listed in each of Tables 2 and 3 were chosen and FISH images for those genes was retrieved from the FlyFISH database. Four different stage categories of early embryogenesis are shown in Figure 5. Green fluorescence represents the mRNA localization pattern for that transcript/gene and red fluorescence is showing the position of nuclei within the organism at that stage. For the known gene belonging to the GO-BP DNA-dependent DNA replication, the annotation http://www.biodatamining.org/content/6/1/8

Figure 5 Fluorescent In Situ Hybridization (FISH) Images of Annotated and Un-annotated Gene mRNAs.
One annotated gene and one un-annotated gene from each of the two GO-BP categories shown in Tables 2 and 3 were chosen and FISH images for those genes was retrieved from the FlyFISH database. Four different stage categories of early embryogenesis are shown. Green fluorescence represents the mRNA localization pattern for that transcript/gene and red fluorescence is showing the position of nuclei within the organism at that stage. CG1584 = Orc6, Origin recognition complex subunit 6. Its molecular function is described as DNA binding and it is involved in the biological processes: DNA replication initiation; DNA-dependent DNA replication; chromatin silencing. CG1970-NADH dehydrogenase (ubiquinone) activity. It is involved in the biological process mitochondrial electron transport.
for gene CG1584 (Gene name = Orc6, i.e. Origin recognition complex subunit 6) has a molecular function described as DNA binding. It is involved in the following biological processes: DNA replication initiation; DNA-dependent DNA replication; chromatin silencing. The un-annotated gene, CG1578 (from Table 2), had a very similar in situ hybridization pattern ( Figure 5A). For the oxidative phophorylation GO-BP category, the known gene is CG1970, annotated as having NADH dehydrogenase (ubiquinone) activity. It is involved in the biological process mitochondrial electron transport. The un-annotated gene, CG1715 (from Table 3), has an in situ hybridization pattern closely resembling the pattern seen for CG1970.
Of the un-annotated genes that are predicted to belong to the DNA-dependent DNAreplication GO-BP genes, 22 had images in the FlyFISH database and 18 out of 22 (82%) had FISH images that were the same or very similar to the known genes belonging to the category (Table 2). Similarly, of the un-annotated genes that are predicted to belong to the oxidative phosphorylation GO-BP category, 19 had images in the FlyFISH database http://www.biodatamining.org/content/6/1/8 and 16 out of 19 (84%) had FISH images that were the same or very similar to the known genes belonging to the category (Table 3).

Discussion
We have applied Support Vector Machines (SVM) and a stratified cross-validation approach to analyze a large microarray experiment dataset from D. melanogaster in order to predict possible functions for previously un-annotated genes. This approach successfully generated a preliminary GO Biological Process annotation for a large number of these genes.
The SVM analysis, employing a sigmoid fitting function, generated annotation probability values for all the genes in the dataset. Using a recall value of 40%, 39 high-precision GO-BP categories were identified and 77% of the un-annotated genes in the dataset have been assigned to one or more of these GO-BP categories. Assigning new annotations based on gene-precision scores takes into account both the overall precision of the GO-BP category, as well as the individual probability of the gene, outcome of SVM and the sigmoid function. In that way, gene predictions of high confidence can be considered even for GO-BP categories that have an overall poor precision performance.
We chose two of the high-precision GO-BP categories to validate using an independent data set (the FlyFISH database). These two categories, DNA-dependent DNA replication and oxidative phosphorylation, were selected for a number of reasons but chiefly because each contained a fairly large number of previously unannotated genes allowing for a larger number of genes that could be cross-checked in the FlyFISH database and that the GO-BP description was fairly precise allowing for the potential to examine whether the unannotated genes are actually involved in the biological process they had been assigned to. While it could be argued that the use of FlyFISH data is perhaps too similar of a data type to be useful (i.e. is also gene transcription data), our main goal was to use data that: was from a completely independent source; had at least some data that differed in type from the original microarray data (in this case spatial data in addition to temporal data); where the raw data could be inspected independently at the gene level so that any computer or user assigned functional classification of the genes would not bias the validation. Given the high degree of precision that some of the unannotated genes have for belonging to a specific GO-BP category, especially for the oxidative phosphorylation category, we feel quite confident that a large number of these genes will be shown in the future to be directly or indirectly involved in these biological processes. Gene expression pattern data has been used in the past to uncover the function of a gene. In a study by Hughes and coworkers [14], a large compendium of yeast gene expression data was examined for co-expression patterns using hierarchical clustering and among several findings, they discovered that a previously unannotated gene, YER044c (now ERG28), had a gene expression pattern that was highly similar to genes that were known to be involved in ergosterol biosynthesis. Further biochemical and genetic experiments supported that the gene was indeed involved in this process [14].
In the original Arbeitman and coworkers study [22], the authors performed different hierarchical clustering analyses to identify groups of genes that were temporally coregulated and/or restricted to specific biological sample type (e.g. mutant or sex). The authors also identified several GO classes that could be assigned to subsets of the coregulated genes. While many of the GO classes reported were fairly low (i.e. general) in http://www.biodatamining.org/content/6/1/8 the GO annotation hierarchy (e.g. enzyme, metabolism) others were more detailed. Five of the more detailed GO categories reported by Arbeitman and coworkers [22] overlap with the 39 GO biological categories identified in this study including: mitochondria/electron transport; ribosome; cell cycle; cell adhesion; and phototransduction.
We attempted to compare our results with those from the study of Yan and co-workers [10], who also assigned GO terms to groups of similarly behaving genes to a dataset which in part included the data we used in our study. The authors of that study apply a different computational method (Random Forest) to a number of data sources (including gene expression data) for the prediction of GO-BP annotations. For the comparison with our study, first, predictions based on expression data only have to be extracted and, subsequently, precision and recall measures need to be calculated using methods similar to those used in our study. There were a few issues that made the comparison between the two studies difficult. First, Yan and co-workers do not apply stratification in their 10fold cross-validation. As a result, for some of the GO-BP categories (and the associated classifiers) there are a number of folds that contain only negative classes, and therefore classification (and calculation of precision) cannot be performed. Estimates of measures such as precision at 40 could be generated based on the remaining folds, but this practice is generally not recommended as it may introduce bias [33]. Second, the number of categories that are represented in both studies (273) is much smaller than those that are represented in each study individually. If we exclude those that give folds without positively labeled representatives, we are left out with a very small number of GO-BP categories for comparison. Due to these issues a strict comparison of the results from the two studies is not particularly valuable. For reference purposes we mention here that in Yan et al. study, 8 GO-BP categories gave a precision at 40 value equal or larger than 0.75 after excluding the GO terms with less than 5 folds representing both classes, while the number of selected categories dropped to 1 if we exclude those with less than 8 folds representing both classes. Also, the reported new predictions in Yan et al. involve 2062 genes at the confidence level of 0.2, 213 genes at the confidence level of 0.5 and 11 genes at the confidence level of 0.75. It should be noted, though, that the confidence level is Yan et al. study was computed differently from the gene precision score in our study and the two are not necessarily comparable.
We believe that the methods used in this study can be applied to other microarray and RNA-seq datasets. Our own experience trying this approach on other data suggests that the robustness of the results is largely dependent on both the number of array experiments that are in the dataset as well as on the type of experiment(s) the data represent. In this case, a developmental time course with numerous time points provided the ability to identify a large number of co-expressed genes.

Conclusions
In this study we have applied Support Vector Machines for the prediction of GO-BP annotation terms for D. melanogaster genes, using published time-course gene expression data. We have assessed the predictive ability of the method using an elaborate stratified double cross-validation procedure, involving the fitting of a sigmoid function on the raw output from SVM. Precision and recall values for each GO-BP function were calculated, describing in detail the performance of the prediction. On a second level, our results were validated externally using independent data from FlyFISH database. Finally, we applied http://www.biodatamining.org/content/6/1/8 SVM for the prediction of GO-BP annotations of un-annotated genes, providing certainty estimates for those predictions, based on a method that uses the results from the cross-validation and the predictive performance of SVM for the specific GO-BP category.
We believe that this is the first study that investigates the performance of SVM for GO-BP prediction for D. melanogaster genes and provides a comprehensive list of new annotation predictions that can be further researched experimentally. We believe that our study is of great value to researchers interested in computational annotation of genes, as well as to the community of D. melanogaster researchers.

Additional files
Additional file 1: The Parent-Child Relationship Between Ontology Biological Process Categories. Figure showing a joined Gene Ontology graph providing an overview of the relationship between the 39 Gene Ontology Biological Process (GO-BP) categories that were identified from the output of the SVM as having a precision-at-40 value equal or larger than 0.75, as indicated in the highlighted polygonal nodes. Table showing all new predictions with average gene-precision score equal or larger than 0.75. Detailed information on gene-precision scores and probability estimates for each one of the four folds is shown.

Additional file 2: Selected GO-BP Predictions for Previously Unannotated Genes.
Abbreviations GO-BP: Gene Ontology Biological Processes; SVM: Support Vector Machine; CG number: computer gene number.


extfieldsofstudy: Biology Medicine

text:
Effect of FK506 in reducing scar formation by inducing fibroblast apoptosis after sciatic nerve injury in rats

We previously demonstrated that FK506, a generally applied immunosuppressant in organ transplantation, could promote peripheral nerve regeneration through reducing scar formation. However, little is known about how FK506 reduces scar formation. Herein we investigated the influence of FK506 on fibroblast proliferation and its correlation with scar formation after sciatic nerve injury in rats, and further explored the effect of FK506 on fibroblast proliferation and apoptosis in vitro. Masson staining and immunohistochemistry revealed that scar area and fibroblast number in the nerve anastomosis of sciatic nerve-injured rats were significantly reduced after FK506 administration. The scar area had a significant positive correlation with the fibroblast number, as detected by linear correlation analysis. CCK-8 assay and flow cytometry indicated that FK506 also inhibited proliferation and induced apoptosis of fibroblasts in vitro. It was primarily phosphorylation of JNK and ERK that were activated during the apoptosis of fibroblast. Pretreatment of cells with JNK inhibitor, SP600125, or ERK inhibitor, PD98059, could inhibit FK506-induced fibroblast apoptosis, respectively. Moreover, simultaneous application of both inhibitors had additive roles in cell protection from apoptosis. These results suggest that FK506-induced fibroblast apoptosis contributes to the suppression of fibroblast proliferation and then results in the reduction of scar formation in sciatic nerve-injured rat, and that JNK and ERK are involved in FK506-induced fibroblast apoptosis.

Following traumatic peripheral nerve injury, proliferation of scar tissue from the epineurium will result in impediment to the regenerating axons that need to traverse the injury site. [1][2][3][4][5] A reduction in scar formation at a site of nerve repair is associated with better recovery of neural function. 6,7 It has been demonstrated that FK506, a generally applied immunosuppressant in organ transplantation, has powerful effect of promoting axon regeneration [8][9][10][11][12] through its immunosuppressive action 13 and neurotrophic action. [14][15][16][17][18] Moreover, our previous study 19 found that FK506 could promote peripheral nerve regeneration and accelerate the recovery of neural function by reducing scar formation. It is believed that increased collagen biosynthesis from the proliferative fibroblasts contributes to scar formation at the repair site, and preventing fibroblast proliferation can reduce scar area. 20,21 In order to gain further insight into the mechanism of action of FK506, this study investigated the effect of FK506 on fibroblast proliferation and its correlation with scar formation in sciatic nerve-injured rat model, and further explored the influence of FK506 on fibroblast proliferation and apoptosis in vitro and its possible cell signaling. The results from this study may provide further explanation for the effect of FK506 on promoting peripheral nerve regeneration.

Results
Effect of FK506 on scar formation and fibroblast proliferation in histological analysis. At 4 weeks after establishing the nerve injury model, Masson staining showed extensive collagen tissue hyperplasia in nerve anastomosis of the model group ( Figure 1A), as seen by the blue staining. Immunohistochemistry showed that a large amount of fibroblasts mixed with fibrocytes appeared in the model group ( Figure 1a). In the FK506 group ( Figures 1B and b), collagen tissue hyperplasia and fibroblasts mixed with fibrocytes obviously decreased in nerve anastomosis compared with those of the model group. Meanwhile, collagen tissue and fibroblasts were rarely visible in nerve anastomosis of the normal control group (Figures 1C and c ).
Effect of FK506 on fibroblast number in sciatic nerve anastomosis of rats and its correlation with scar area. As shown in Figure 2a, the scar area of nerve anastomosis was significantly increased in the model group (4.36±0.68 mm 2 ) compared with the normal control group (0.37 ± 0.14 mm 2 , Po0.01), and significantly decreased in the FK506 group (1.61 ± 0.42 mm 2 ) compared with the model group (Po0.05). As shown in Figure 2b, the fibroblast number of nerve anastomosis was significantly increased in the model group (37.3 ± 4.85) compared with the normal control group (3.8±1.03, Po0.01), and significantly decreased in the FK506 group (16.5 ± 3.24) compared with the model group (Po0.05). Linear correlation analysis results showed that the scar area of nerve anastomosis augmented with increased fibroblast number. Thus, the scar area had a significant positive correlation with the fibroblast number (r ¼ 0.9915, Po0.0001), suggesting that FK506 reduces scar formation by inhibiting fibroblast proliferation (Figure 2c).
FK506 inhibited fibroblast proliferation in vitro in a dose-dependent manner. Rat skin fibroblasts were treated with FK506 at increasing concentrations for 8 h. The Cell Counting Kit-8 (CCK-8) assay demonstrated that FK506 could induce a dramatic loss in the viability of fibroblasts. Cell viability reached a relatively minimal level at 75 mM ( Figure 3). FK506 is an effective inhibitor of fibroblast proliferation in vitro and this inhibitory effect is dose-dependent.
FK506 induced significant apoptosis of fibroblast. As the CCK-8 assay demonstrated that FK506 could inhibit the proliferation of fibroblast in a dose-dependent manner, flow cytometry was further used to identify the types of cell death. As shown in Figure 4, fibroblasts showed a dosedependent apoptosis, including early as well as late apoptotic cell death. The mean percentage of apoptotic cell death was about (14.7 ± 6.89), (45.7 ± 7.05) and (86.8 ± 6.83)% in cells treated with 12.5, 25 and 50 mM FK506, respectively. In contrast, the cells that had not been treated with FK506 demonstrated normal cell viability without abnormal cell death.
Morphological changes by fluorescence microscopy. In order to characterize the pattern of apoptotic status mediated by FK506 in fibroblasts, morphological examinations were performed. As shown in Figure 5, normal fibroblasts exhibited intact nuclei and adqulis chromatin. In comparison, after FK506 treatment and fluorescence microscopy using Hoechst 33342 staining, shrunken cells with condensed or fragmented nuclei were observed. In cells treated with 50 mM  FK506, significant chromatin agglutination and nuclear fragmentation, both characteristic of apoptosis, were observed.
FK506 activated JNK and ERK, and enhanced expressions of cytosolic cytochrome c and cleaved caspase-3. JNK, ERK, cytochrome c and cleaved caspase-3 are all known to be involved in the regulation of apoptosis. Their expression levels after FK506 treatment for 8 h were detected by western blotting. As shown in Figure 6, the levels of GAPDH expression were comparable among the negative control group, the dimethyl sulfoxide (DMSO) group and the three FK506 treatment groups. In contrast, phosphorylation of JNK (p-JNK) could not be activated in either the negative control group or the DMSO group. Increasing expression of p-JNK was observed in fibroblasts after FK506 treatment at increasing concentrations; this expression peaked at the concentration of 50 mM. This suggests that JNK can be activated by FK506. The kinetics and the extent of phosphorylation of ERK (p-ERK), cytosolic cytochrome c and cleaved caspase-3 had a similar pattern of increase   JNK inhibitor, SP600125, or/and ERK inhibitor, PD98059, reduced FK506-induced apoptosis of fibroblasts. The effects of JNK inhibitor, SP600125, and ERK inhibitor, PD98059, were investigated to demonstrate the role of JNK and ERK in FK506-induced apoptosis. As shown in Figure 7, preincubation with SP600125 or PD98059 for 30 min before 50 mM FK506 treatment significantly reduced the apoptotic rate. The mean percentage of apoptotic cell death decreased to (61.1 ± 7.18) and (67.6±7.59)% in cells pretreated with SP600125 and PD98059, respectively, compared with that of cells treated with 50 mM FK506 alone (87.9 ± 6.69)% (Po0.05). In addition, preincubation with the two inhibitors together further decreased the apoptotic percentage to (39.87± 4.57)%, compared with SP600125 or PD98059 pretreatment, respectively (Po0.05).
JNK inhibitor, SP600125, or/and ERK inhibitor, PD98059, weakened expressions of p-JNK or p-ERK and cleaved caspase-3. As shown in Figure 8, SP600125 prevented both the phosphorylation of JNK and the expressions of cytosolic cytochrome c and cleaved-caspase-3. However, PD98059 only prevented the phosphorylation of ERK and the expression of cleaved caspase-3, and had no effect on the expression of cytosolic cytochrome c. Moreover, simultaneous application of both inhibitors further decreased the expression of cleaved caspase-3 (Po0.05), and almost had no influence on the expressions of p-JNK, p-ERK and cytosolic cytochrome c, compared with SP600125 or PD98059 pretreatment, respectively.

Discussion
Following traumatic peripheral nerve injury, fibroblasts were activated by inflammatory cytokines and growth factors, such as transforming growth factor-b (TGF-b) and basic fibroblast growth factor. Proliferative fibroblasts produce a great amount of collagen fibers that contribute to scar formation at the nerve anastomotic sites, and preventing fibroblast proliferation can reduce scar area. 20,21 In this study, multiple parameters including the histological analysis, the area of scar tissue and the number of fibroblasts were used to evaluate the effect of FK506 on fibroblast proliferation and scar formation in rats' model after sciatic nerve injury. Furthermore, we performed successive transversal sections through the anastomotic sites to make the scar area analysis and fibroblast counting in histological evaluation. Therefore, the mean values could reflect the amount of scar area and fibroblast number, and it was a valid method to the quantitative evaluation of fibroblast proliferation and scar formation.
With FK506 gavage (4 mg/kg per day) for 4 weeks after sciatic nerve injury, the scar area and fibroblast number were significantly decreased in the FK506 group compared with the model group. These results showed that FK506 could inhibit fibroblast proliferation and reduce scar formation at the nerve anastomotic sites. Moreover, correlation analysis indicated that the scar area was positively correlated with the fibroblast number, which further confirmed that FK506 can reduce scar formation in some way by inhibiting fibroblast proliferation after sciatic nerve injury in rats.
The study further investigated the influence of FK506 on fibroblast proliferation in vitro. The CCK-8 assay demonstrated that FK506 could inhibit the proliferation of cultured rat skin fibroblasts in a dose-dependent manner. Then the effect of FK506 on fibroblast apoptosis in vitro was examined. Hoechst 33342 staining showed that the fibroblasts subjected to FK506 exhibited condensed or fragmented nuclei. Furthermore, flow cytometric analysis showed a significant increase, in a dose-dependent manner, in the percentage of apoptotic cells. The more the apoptotic cells, the less the cells available for proliferation. Thus, these results suggest that FK506-induced fibroblast apoptosis contributes to the suppression of fibroblast proliferation and then results in the reduction of scar formation in sciatic nerve-injured rat.
Furthermore, we explored the cell signaling that may be involved in FK506-induced fibroblast apoptosis. We found that fibroblast apoptosis induced by FK506 could be inhibited by JNK and ERK are members of the mitogen-activated protein kinase (MAPK) family. JNK can be rapidly activated by various stress stimuli and mostly mediates apoptosis, 22,23 whereas the ERK cascade is activated by mitogenic stimuli, and is critical for proliferation and survival. 24 However, JNK and ERK are controversial in their role in cell death or proliferation. JNK has been reported to induce pro-survival responses in several types of cells [25][26][27] and ERK also has been suggested to mediate apoptosis. [28][29][30][31] What determines whether JNK and ERK will act in a pro-apoptotic or pro-survival fashion remains an important unanswered question, but the different stimuli as  well as the cell type may be determining factors. The role of MAPKs on apoptotic pathways may be unique to the specific model tested, and to the apoptotic trigger assessed. 32 This study found that FK506 induced sustained activation (phosphorylation) of JNK and ERK in fibroblasts in a dosedependent manner. JNK mostly mediates apoptosis, whereas the accumulation of p-ERK in the cytoplasm raises the possibility that kinase may not be available to act on potential downstream pro-survival targets and may contribute to apoptosis through the activation of cytoplasmic cell death mediators. 33 SP600125 or PD98059 effectively blocked FK506-induced JNK or ERK activation and attenuated FK506-induced apoptosis, suggesting that the sustained activation of JNK and ERK signaling pathways may be involved in FK506-induced fibroblast apoptosis. Western blot analysis of cytosolic cytochrome c protein was performed to examine whether FK506 induced fibroblast apoptosis involving the mitochondria. Mitochondria act as an important apparatus for signals during apoptosis, and the loss of mitochondrial membrane integrity can induce the release of cytochrome c from mitochondria into cytosol. 34 The release of cytochrome c is important in this apoptotic pathway, as cytochrome c binds to apoptotic proteaseactivating factor-1 (Apaf-1) forming a complex that in turn cleaves caspase-3. 35 In this study, the level of cytosolic cytochrome c significantly increased in a dose-dependent manner after FK506 treatment. After preincubation with JNK inhibitor SP600125, the p-JNK and cytosolic cytochrome c were markedly decreased. Preincubation with ERK inhibitor PD98059, however, resulted in a significant decrease of p-ERK, leaving the level of cytosolic cytochrome c almost unchanged. These results indicated that mitochondria have a critical role and that p-JNK, but not p-ERK, is involved in mitochondrial pathway of FK506-induced fibroblast apoptosis.
Caspase-3 is one of the key agents of apoptosis, as it is either partially or totally responsible for the proteolytic cleavage of many key proteins. 36 Although p-ERK had no effect on the release of cytochrome c, it still contributed to the activation of caspase-3. Preincubation with SP600125 or PD98059 both led to a significant decrease of cleaved caspase-3. Moreover, simultaneous application of the two inhibitors had additive roles that further decreased the expression of cleaved caspase-3 and the apoptotic percentage, and almost had no influence on the expressions of p-JNK, p-ERK and cytosolic cytochrome c, compared with SP600125 or PD98059 pretreatment, respectively. These observations further demonstrated that JNK and ERK pathways are independently involved, but have additive roles in FK506-induced fibroblast apoptosis.
In conclusion, FK506 activates both JNK and ERK with or without mitochondrial cytochrome c release, followed by the cleavage of caspase-3, subsequently leading to the apoptosis of fibroblasts, and FK506 has a valid effect on reducing scar formation in sciatic nerve-injured rat by inducing fibroblast apoptosis. Further research is needed to study at the effect of FK506 on other extracellular matrix components in the process of scar formation. Experiments will be continuously performed in vitro to further elucidate the signaling pathways of FK506-induced apoptosis in fibroblasts.

Materials and methods
Animals and establishment of rat sciatic nerve transection anastomosis models. Healthy male Sprague-Dawley rats, weighing 250-300 g were purchased from Shanghai Laboratory Animal Center, the Chinese Academy of Sciences (Shanghai, China). All procedures were approved by the Animal Ethics Committee of Nanjing Medical University in China and were in compliance with the Guide for the Care and Use of Laboratory Animals published by the Chinese National Institute of Health. Rats were anesthetized by 30 mg/kg sodium pentobarbital via intraperitoneal injection and fixed in an abduction position at room temperature. After the right thigh and buttock were shaved, a longitudinal incision was made along the median femur to expose B2 cm of the sciatic nerve. The greater sciatic foramen was cut at the tendon level of the obturator internus muscle. The incision was sutured with three sutures under Â 10 magnification using a surgical microscope (Latric Optical Instrument Factory, Shanghai, China) and 10/0 nylon monofilament (Latric Suture Needle Company Ltd., Ningbo, China). Rats were postoperatively injected with penicillin 1 600 000 U/kg at the non-operated side to prevent infection.
Drug intervention in vivo. FK506 was purchased from Sigma (St. Louis, MO, USA), containing 1 mg per capsule. After the capsule dressing was removed, the powder was dissolved in saline and prepared into a homogeneous suspension with a final concentration of 0.5 mg/ml by microwaving for 5 min. Rats were intragastrically injected with the suspension at a dose of 4 mg/kg per day after model establishment. Thirty rats were randomly divided into three groups (ten rats in each group): model group, saline gavage for 4 weeks after sciatic nerve injury; FK506 group, FK506 gavage for 4 weeks after sciatic nerve injury; normal control group, saline gavage for 4 weeks without sciatic nerve injury.
Harvesting the sciatic nerve. At 4 weeks after model establishment, all the rats of each group were sacrificed. The operative sciatic nerve was completely dissected, a 5-mm-long nerve trunk was cut from 2.5 mm lateral to the distal and proximal ends of the anastomosis, followed by labeling, conventional fixation, dehydration, transparency, permeation, embedding, sectioning, coating, paraffin mounting and slicing into 4-mm thick transversal sections. Twelve successive transversal sections of each rat were made through the anastomosis from the proximal to the distal end.
Masson staining. Six odd sections of each rat were stained with Masson's trichrome. Sections were dewaxed and mordanted (10% potassium dichromate þ 10% trichloroacetic acid) for 30 min, nuclei were stained with hematoxylin for 20 min, differentiated with hydrochloric acid and ethanol for 15 s, returned to blue with weak ammonia for 15 s, stained with Masson solution (Cell Signaling Technology, Irvine, CA, USA) for 1 min, rinsed with 1% acetic acid, dehydrated with an increasing ethanol series, cleared with xylene I and II for 10 min to make sections transparent and finally sealed in resin. Then collagen fiber proliferation was observed under the light microscope (Olympus, Tokyo, Japan) with the magnification of Â 400.
Immunohistochemical analysis. Six even sections of each rat were dewaxed and rehydrated through gradient ethanol solutions. The sections were incubated in citrate buffer to activate the antigenicity and exposed to 3% H 2 O 2 for 10 min to block endogenous peroxidase. All subsequent incubations were performed in a humidified chamber at 37 1C unless otherwise stated. The sections were blocked in common goat serum solution for 15 min, followed by incubation overnight with rabbit polyclonal anti-TGF-b (1 : 50) at 4 1C. Sections were then washed in phosphate buffered saline (PBS) and incubated with horseradish peroxidase-conjugated goat anti-rabbit IgG antibody for 15 min. The cells were visualized by using 3,3 0 -diaminobenzidine for 5 min at room temperature, and the sections were counterstained with hematoxylin for 2 min. The sections were analyzed under light microscope with the magnification of Â 400.
Histological image analysis. Six sections of each rat stained with Masson's trichrome were selected and the histological images at Â 12.5 magnification were obtained. The areas of scar tissue were automatically calculated by Image-pro plus 5.0 image analysis software (Media Cybernetics, Shanghai, China). The area of scar tissue for each rat was defined as the mean scar area of six sections and expressed in mm 2 . Another six sections stained with TGF-b antibody of each rat were selected, and the histological images at Â 400 magnification were also obtained. The counting area was selected and each was about 100 mm Â 100 mm from three fields on each section. The number of FK506 induces apoptosis of fibroblast J Que et al fibroblast was calculated, and the number of fibroblast for each section was defined as the mean number of three fields, and the number for each rat was defined as the mean number of six sections.
Culture of rat skin fibroblasts and drug intervention in vitro. Primary cultures of rat skin fibroblasts were obtained from male Sprague-Dawley newborn rats (Shanghai Laboratory Animal Center, the Chinese Academy of Sciences) 1-2 days after birth as reported previously, 37 with minor modifications. In brief, after thiopental anesthesia, the skin of dorsal part of the lumbar region was shaved and disinfected with 70% ethanol. Skin samples were aseptically removed from the rats and then minced and incubated for 3 h in 0.1% collagenase solution (37 1C, 5% CO 2 ). After incubation, the cell suspension was centrifuged at 250 Â g for 5 min. The supernatant was discarded, and the isolated cells were cultured in DMEM (HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (JRH Biosciences, Lenexa, KS, USA), 100 U/ml penicillin G and 100 U/ml streptomycin. After 2 h, the nonadherent cells were washed out. Cells were cultured in an atmosphere of 5% CO 2 at 37 1C. After confluence was attained, cells were subcultivated by a 5-min incubation in 0.1% trypsin and 5 mM ethylenediaminetetraacetic acid (EDTA). The free cells were added to an equal volume of culture medium. Following centrifugation (1000 Â g for 5 min), cells were resuspended in culture medium. The thin and triangular cells with light cytoplasm present in the culture were identified as fibroblasts. Western blotting. Cells were harvested and washed three times with cold PBS. Then cells were lysed with lysis buffer (1% Nonidet P-40, 50 mM Tris, pH 7.5, 100 mM NaCl, 5 mM EDTA, 10 mg/ml aprotinin, 10 mg/ml leupeptin and 1 mM PMSF). Lysates was centrifuged at 12 000 Â g for 30 min at 4 1C. The supernatant was collected and protein concentrations were determined using BCA assay kit (Pierce, Rockford, IL, USA) according to the manufacturer's instructions. Cellular extract containing 50 mg of total protein from each sample were separated by 12% SDS-PAGE and transferred to PVDF membranes using a standard protocol. The membranes were blocked with TBS-0.05% Tween-20 containing 5% (v/v) bovine serum albumin (Sigma) for 90 min at 37 1C, then incubated with primary antibodies overnight at 4 1C. The following antibodies were used: GAPDH, JNK, ERK, p-JNK, p-ERK (Cell Signaling Technology, 1 : 1000 dilution) and cleaved caspase-3 (Abcam, Cambridge, UK, 1 : 200 dilution). The second day, membranes were washed three times and hybridized with horseradish peroxidase-conjugated secondary antibody (Sigma, 1 : 10 000 dilution) for 2 h at 37 1C. After washing three times, protein bands specific for antibody were visualized by fluorography using an enhanced chemiluminescence kit (Amersham Pharmacia Biotech, Piscataway, NJ, USA). Band densities were quantified using Image-J, a public domain Java image-processing program inspired by NIH and were normalized to GAPDH to correct for variations in protein loading.
Measurement of cytosolic cytochrome c. Cells were resuspended in buffer (250 mM sucrose, 20 mM HEPES, 10 mM KCl, 1.5 mM MgCl 2 , 1 mM EDTA, 1 mM EGTA, 1 mM DTT, 17 mg/ml of PMSF, 8 mg/ml of aprotinin and 2 mg/ml of leupeptin, pH 7.4) on ice for 30 min and then homogenized gently with a glass Dounce homogenizer for 20 strokes. The homogenates were centrifuged at 750 Â g for 10 min at 4 1C to remove nuclei, and the supernatants were centrifuged at 20 000 Â g for 30 min at 4 1C to remove mitochondria and other insoluble fragments. Protein (50 mg) was subjected to 12% SDS-PAGE, and then western blot analysis was performed using anticytochrome c antibody (Abcam, 1 : 500 dilution).
Statistical analysis. Data are presented as mean ± S.D. Statistical significance was determined using analysis of variance (ANOVA). A P-value of o0.05 was considered statistically significant. The correlation between scar area and fibroblast number was analyzed using linear correlation analysis, with Po0.05 indicating correlation.

Conflict of Interest
The authors declare no conflict of interest.


s2fieldsofstudy: Biology Medicine

text:
Complement alternative pathway genetic variation and Dengue infection in the Thai population

Dengue disease is a mosquito-borne infection caused by Dengue virus. Infection may be asymptomatic or variably manifest as mild Dengue fever (DF) to the most severe form, Dengue haemorrhagic fever (DHF). Mechanisms that influence disease severity are not understood. Complement, an integral component of the immune system, is activated during Dengue infection and the degree of activation increases with disease severity. Activation of the complement alternative pathway is influenced by polymorphisms within activation (factor B rs12614/rs641153, C3 rs2230199) and regulatory [complement factor H (CFH) rs800292] proteins, collectively termed a complotype. Here, we tested the hypothesis that the complotype influences disease severity during secondary Dengue infection. In addition to the complotype, we also assessed two other disease-associated CFH polymorphisms (rs1061170, rs3753394) and a structural polymorphism within the CFH protein family. We did not detect any significant association between the examined polymorphisms and Dengue infection severity in the Thai population. However, the minor allele frequencies of the factor B and C3 polymorphisms were less than 10%, so our study was not sufficiently powered to detect an association at these loci. We were also unable to detect a direct interaction between CFH and Dengue NS1 using both recombinant NS1 and DV2-infected culture supernatants. We conclude that the complotype does not influence secondary Dengue infection severity in the Thai population.

Introduction
Dengue disease, a mosquito-borne infection, is a public health concern in many subtropical and tropical countries. The causative agent is Dengue virus (DV), a positive-sense single-stranded RNA virus and a member of the Flaviviridae family. This virus can be classified into four serotypes (DV 1-4) whose genomes are closely related but distinct. The genome expresses three structural proteins (Capsid, PrM and Envelop protein), and seven nonstructural proteins (NS1, NS2a, NS2b, NS3, NS4a, NS4b and NS5). DV infection may be asymptomatic or manifest as undifferentiated fever, Dengue fever (DF) or Dengue haemorrhagic fever (DHF) with plasma leakage. DHF can be classified on severity into four subtypes (DHF 1-4). The most severe disease, which has been termed Dengue shock syndrome (DSS), comprises a syndrome of plasma leakage with shock and severe bleeding [1]. Both the host immune response and the DV strain contribute to the development of DHF [2][3][4]. Non-neutralizing antibodies enhance infectivity in an Fc-dependent manner [5]. A number of host susceptibility genes have been identified and include genes encoding: Fc gamma receptor II (FcγRIIA), human leucocyte antigens (HLA-A1, and HLA-A2), the vitamin D receptor, tumour necrosis factor (TNF)-α and dendritic cell-specific intercellular adhesion molecule-3-grabbing non-integrin (DC-SIGN) [6][7][8].
The complement system is an integral component of host defence and appears to influence the host response to flavivirus infection in murine models. Mice lacking complement components [C3, C1q, C4 or factor B (FB)] were more susceptible to West Nile virus (WNV) infection [9]. Activation of complement increased antibody-mediated neutralization of flavivirus, as C1q reduced the antibody-dependent enhancement of both WNV and DV infection [10]. Complement activation occurs during the course of DV infection. An increase in C3 activation, together with reduced levels of C4 and FB, has been shown in DHF [11]. C3 levels were reduced while the levels of C3a and C5a (anaphylatoxins) were increased in DHF compared to DF [12]. C5b-9 complexes, along with high levels of DV NS1 protein, correlated with DV infection severity [13].
The molecular mechanisms through which complement could influence flavivirus infection are incompletely understood. DV NS1 can bind complement C4, a classical and lectin pathway activation protein and trigger C4 degradation [14]. WNV NS1 binds the complement regulatory protein, complement factor H (CFH), an interaction that prevented the lysis of WNV-infected cells [15]. In contrast, Japanese encephalitis virus (JEV) NS1 did not bind CFH and did not modulate the host innate immune response [16].
Recently, it has been shown that polymorphisms (collectively termed 'complotype') in complement proteins C3, FB and CFH can be associated with an up to sixfold variation in complement haemolytic activity in vitro [17]. C3 and FB are activation proteins of the complement alternative pathway, while CFH is the major regulator of this pathway. Therefore, the description of the complotype demonstrated that alternative pathway activation is influenced by genetic factors. These polymorphisms have been associated with susceptibility to age-related macular degeneration (AMD [18][19][20]), suggesting that the degree of haemolytic activation driven by these polymorphisms in vitro is biologically relevant. Individuals with the 'high activation' complotype have increased susceptibility to AMD [FB-32 R/R: rs12614 (C) and rs641153 (G); CFH-62 V/V: rs800292 (G); and C3-102 F/F: rs2230199 (G)], whereas risk is reduced in those with the 'low activation' complotypes [FB-32 Q/Q rs12614 (C) and rs641153 (A); CFH-62 I/I: rs800292 (A); and C3-102 S/S: rs2230199 (C)]. Consistent with these data was the finding that serum complement activation appears to be enhanced in AMD [21].
In addition, genetic variation across the CFH-CFHR locus (the region encoding CFH and CFH-related proteins) is linked to AMD susceptibility. This variation includes a polymorphism within CFH (CFH-402 Y/H, rs1061170, [22]) and a copy number polymorphism that results in deletion of the CFHR1 and CFHR3 genes (ΔCFHR3-1, [23]). Recently, CFH promoter (-332 C/T, rs3753394) and CFH (V62I, rs800292) polymorphisms have been associated with Dengue disease [24]. The functional consequences of these polymorphisms are not understood.
Although studies have implicated a role for complement in DV infection severity, only one study to date has explored the relationship between complement genetic variation and DV infection [24]. In this study we specifically analysed the relationship between DV infection severity and the complotype (CFH rs800292, FB rs12614/rs641153 and C3 rs2230199), as this variation has been associated with differences in complement activation. We also included the AMD-associated CFH rs1061170 variant, the CFH promoter rs3753394 variant and the ΔCFHR3-1 deletion polymorphism. We did not find any strong correlation between DV infection severity and these specific polymorphisms. Our data demonstrate that variation in these polymorphisms is not a major host factor for determining the phenotypic response to this infection.

Patients and controls
Genomic DNA samples were extracted from patients attending either Khon-Kean or Song-Khla Hospitals, Thailand. Viral diagnosis was established using reverse transcription-polymerase chain reaction (RT-PCR) and severity classified using World Health Organization (WHO) guidelines: DF and DHF grades 1-4. Dengue immunoglobulin (Ig)G/IgM capture enzyme-linked immunosorbent assay (ELISA) was used to classify primary and secondary infection. An IgM : IgG ratio equal to or greater than 1·8 was used to define primary infection [25]. Patients with other febrile illnesses (OFI) (n = 63) were those presenting with fever and having undetectable Dengue viral genome and levels of anti-Dengue antibodies lower than the cut-off value for an acute Dengue infection. Healthy volunteers were recruited from Siriraj Hospital, Thailand. Written informed consent was obtained and sample collection approved by the Institutional review board.

Genotyping
Genomic DNA was amplified by PCR using high-fidelity polymerase (Certamp; Biotools, Madrid, Spain) and primers listed in Table S1. Following purification (ChargeSwitch PCR clean-up kit; Invitrogen, Carlsbad, CA, USA) amplicons were sequenced and single nucleotide polymorphism (SNP) genotype-recorded. The CFH rs800292 variant was determined using the Taqman SNP genotyping assay (no. 4351379; Applied Bioscience, Island, NY, USA). ΔCFHR3-1 copy number was detected using a combination of (i) Taqman copy number assay probes (Applied Bioscience), (ii) genotyping of rs6677604 (a SNP within intron 12 of the CFH gene that has been shown to tag the ΔCFHR3-1 deletion) [26] and (iii) multiplex ligation-dependent probe amplification assay (MLPA; MRC Holland, the Netherlands). MLPA was performed according to the manufacturer's instructions. PCR products were mixed with GeneScan 500 ROX Size Standard (Applied Bioscience) in Hi-Di formamide (Applied Bioscience) and separated using capillary electrophoresis (3730xl DNA Analyser; Applied Biosystems, Grand Island, NY, USA).

Dengue virus infection in vitro
C6/36 cell lines were cultured in 10% fetal bovine serum (FBS), 10% tryptose-phosphate-broth (TPB) Leibovitz L-15 with streptomycin and penicillin at 28°C. Vero cell lines were maintained with modified Eagle's medium (MEM) supplemented with 10% FBS at 37°C in a CO2 incubator. Dengue virus serotype 2 strain 16681 (DV2) was propagated in C6/36 cells. Infected cells were incubated further at 28°C until a cytopathic effect (CPE) could be seen. Mock was used as a negative control. To titrate the virus in the supernatant, Vero cells were plated on each well in a 96-well plate and incubated for 2 days in a CO2 incubator. On the day of infection, the supernatant were serially diluted 10-fold with 3% MEM and added to the cells for 2 h at 37°C. The infected cells were then overlaid with 1·5% carboxymethylcellulose and incubated further for 2 days. Infected foci were observed by immunoperoxidase staining. Briefly, the infected cells were washed with PBS, fixed and permeabilized with 3·7% formaldehyde and 2% Triton-X 100, respectively. Monoclonal antibody against Dengue E protein (4G2) and rabbit anti-mouse-horseradish peroxidase (HRP) diluted in 1/1000 (Dako, Glostrup, Denmark) were added sequentially before adding 3,3′diaminobenzidine tetrahydrochloride (DAB) (Sigma-Aldrich, St Louis, MO, USA). The number of foci were counted and reported as foci-forming units per ml.

CFH and NS1 interaction assays
To determine if CFH interacted with NS1 we utilized ELISA, far Western and immunopreciptiation techniques.
For far Western blotting, mock and DV2 supernatant was separated using sodium dodecyl sulphate-polyacrylamide gel electrophoresis (SDS-PAGE). Following membrane transfer, membranes were incubated with 5 μg/ml CFH. Following washing, either goat anti-human factor CFH (Quidel) or 1F11 was used to detect CFH and NS1, respectively. CFH (0·5 μg) was also applied to non-reducing gels, and following transfer, membrane-incubated with either mock or DV2 supernatant. Following washing, either goat anti-human factor CFH (Quidel) or 1F11 was used to detect CFH and NS1, respectively. Recombinant NS1 (rNS1; Abcam, Cambridge, UK) was also tested for interaction with CFH.
For co-immunoprecipitation, mock or DV2 supernatant was incubated with 5 μg/ml CFH for 3 h at 4°C followed by 1F11 for 1 h. Protein A (Roche, San Francisco, CA, USA) was then added and the mixture was rotated overnight. After centrifuging, beads were washed three times using 0·1% PBS-T, bound protein eluted using non-reducing loading buffer and eluate analysed by SDS-PAGE. NS1 and CFH were detected using 1F11 and goat anti-human CFH antibody, respectively (Quidel).

Statistical analyses
The analyses were performed on SPSS version 19·0 (SPSS, Inc., Chicago, IL, USA). Hardy-Weinberg equilibrium and association tests were carried out by using the χ 2 test. Odds ratio (OR) was the odds of the genotype comparing either DHF to DF or DHF with shock to DHF without shock. Logistic regression model was used to adjust OR with sex and age. Bonferroni's correction was performed by dividing α-level with the total number of polymorphisms compared (n = 6) to obtain an adjusted P-value of 0·05/6 = 0·008.

Study population
The study population demographics are depicted in Table 1. The age of patients with either Dengue infection or non-Dengue febrile infection was 2-15 years. The range for healthy controls was higher, at 19-47 years. The number of males and females within each of the three groups was comparable with the exception of the DHF4 subgroup, which was comprised of two females. We included only patients with secondary infection.

Complement gene polymorphisms did not differ in frequency between patients with severe and non-severe Dengue infection in the Thai population
The individual SNP frequencies in FB (rs12614, rs641153), C3 (rs2230199) and CFH (rs3753394, rs1061170 and rs800292) did not differ between patients with Dengue infection and either healthy controls or individuals with non-Dengue febrile illness (Table S2). We next examined whether the SNP frequencies differed among Dengue infection subgroups, defined as those with DF (n = 121) and those with DHF (n = 187). We also compared genotypes between DHF patients with (n = 43, defined as DHF3 and DHF4) and without shock (n = 144 defined as DHF1 and DHF2, Table 2). No differences were noted for any of the tested polymorphisms. Notably, the minor allele frequency of functional C3 rs2230199 polymorphism was low (2%), as were the minor allele frequencies of the adjacent rs12614 and rs641153 FB polymorphisms (<10%). In contrast, the CFH rs800292 polymorphism had a minor allele frequency of 42% in this population, but nevertheless its frequency did not differ between Dengue infection subgroups. Due to the rarity of variation at C3 rs2230199, FB rs12614 and FB rs641153FB we did not detect any individual in our Dengue infection cohort with either the most [C3F/F, FB-32R/R (C/G at the adjacent rs12614/rs641153) and CFH 62V/V] or least [C3S/S, FB-32Q/Q (C/A at the adjacent rs12614/ rs641153) and CFH 62I/I] active complotype described in vitro (Table 3) [17]. One high (defined as FB-32R/R and CFH 62V/V) and one low (defined as CFH 62I/I and C3S/S) activation complotype were present in our Dengue infection cohort, but neither was associated with infection severity ( Table 3).

The CFHR3-1 deletion polymorphism is not associated with Dengue infection severity
The ΔCFHR3-1 deletion allele frequency in our controls (healthy and non-Dengue febrile illness) and Dengue infection was 6·79 and 6·94%, respectively, and did not differ significantly between these two groups or between subgroupings of Dengue infection (Table 4). In this analysis, we noted one individual with three copies of CFHR3-1, suggesting duplication of these genes on one allele. In addition, nine individuals had heterozygous deletion of the CFHR1 gene in the presence of two copies of the CFHR3 genes. This is most likely to indicate the presence an allele containing deletion of both the CFHR1 and CFHR4 genes (ΔCFHR4-1), a variant reported in healthy controls, but at a lower frequency than that of the ΔCFHR3-1 deletion allele [27].
Dengue NS1 protein did not bind to the complement factor H protein NS1 protein of WNV virus has an immune evasion activity by binding to CFH protein, resulting in decreased comple-ment recognition of infected cells [15]. We next determined if Dengue NS1 interacts with CFH using either recombinant NS1 or NS1 from DV2-infected cell supernatants (Fig. 1). Using a combination of ELISA, far Western blotting and immunoprecipitation assays we could not detect any evidence of an interaction between CFH and Dengue NS1 (Fig. 1). As a positive control for the far Western blotting, we utilized the known interaction between CFH and the meningococcal factor H-binding protein (Fig. S1).

Discussion
In this study we assessed the relationship of selected complement gene polymorphisms and Dengue infection severity. We included variation within four nucleotides, two adjacent in FB and one in each of C3 and CFH, that have been shown to determine a sixfold variation in complement haemolytic activity in vitro [17]. Our data showed that the minor allele frequency in the C3 and FB polymorphisms was very low, and consequently any population effect of variation at this locus would be minimal and difficult to detect without large numbers. However, although the minor allele frequency at the CFH rs800292 polymorphism was 42% (Table 2), we did not detect any significant variation in the distribution of genotypes among the subgroups of Dengue infection. A recent study of 121 Brazilian patients with DENV-3 (genotype III Sri Lankan-Indian strain) assessed the relationship between four CFH polymorphisms (rs3753394, rs800292, rs3753396, rs1065489) and infection severity [24]. The minor alleles at rs3753394 and rs800292 were associated with reduced infection severity. Our analysis included these loci, but did not replicate the association. Hence, if the association is robust in the Brazilian population, it does not apply in the Thai population.
Retrospective studies of Dengue infection have shown a correlation between the disease severity and the complement activation [11,12]. Plasma CFH was lower in patients with DHF during the acute phase of the illness, returning to normal during convalescence. In contrast CFH levels did not alter during DF infection [12]. However, in this study Healthy and other febrile illness (OFI) subgroups were utilized as non-Dengue controls. DF: Dengue fever; DHF1-4: Dengue haemorrhagic fever 1-4.  we did not find a significant correlation with severity and the examined CFH polymorphisms. The absence of the genes encoding the CFH-related proteins, CFHR1 and CFHR3, is a frequent polymorphism (ΔCFHR3-1), which is thought to have developed through non-allelic homologous recombination within the CFHR locus [22,23]. This allele confers protection against AMD and IgA nephropathy [23,28]. The frequency of the ΔCFHR3-1 allele in our Dengue population was approximately 7% ( Table 4). The reported allele frequencies in other populations include 4·7% in European Caucasians, 16% in African populations and 2% in Asians [22]. In our analysis there was no evidence that this allele conferred protection against Dengue infection or its severity. During copy number analysis across the CFHR locus we noted uncommon but previously reported variation, including one individual with three copies of the CFHR3 and CFHR1 genes, most probably a result of duplication of theses genes on one allele [29,30]. We also detected nine individuals with deletion of the CFHR1 gene but normal copies of the CFHR3 genes. One explanation for these findings is the presence of a less common allele in which there is deletion of the CFHR4 and CFHR1 genes (ΔCFHR4-1). Alternatively, this could represent deletion of the CFHR1 genes in isolation [27,31]. Our copy number analysis did not include probes for CFHR4 gene, so this remains to be shown.
It is important to note that in our study we excluded individuals with primary infection, because in many cases infection is asymptomatic. Secondly, our sample size for severe disease (DHF3 and DHF4) was relatively small, and it remains possible that we could lack the power to detect small effects. Thirdly, the low frequency of the minor allele of the studied C3 and FB polymorphisms means that we cannot conclude what effects variation could have on infection severity in non-Thai populations, where there is greater variation at these loci.
Levels of Dengue soluble NS1 protein correlates with disease severity [32] and Dengue soluble NS1 protein can directly activate the complement cascade, a phenomenon that is enhanced by the presence of anti-Dengue antibodies [13]. Dengue NS1 has been shown to interact with Clusterin, a complement regulatory protein [33]. Recently, some studies have shown that WNV NS1, but not JEV, could bind to CFH, and thereby reduce complement activation on the virus in vitro [15,16]. We were unable to detect a direct interaction between CFH and Dengue NS1 protein, so we currently have no evidence to suggest that NS1 interferes with CFH-mediated regulation of complement during infection, at least in a direct manner.
In summary, our data demonstrate that the severity of Dengue infection is not influenced strongly by the CFH rs3753394 promoter polymorphism, the AMD-associated CFH rs1061170 polymorphism, the AMD-protective ΔCFHR3-1 allele and the complotype-associated CFH rs800292 polymorphism in Thai population. We also show     Dengue non-structural protein 1 (NS1). (a) and (b) enzyme-linked immunosorbemt assays (ELISA) to detect interaction between CFH and NS1. (a) Supernatants were added to plates coated with human CFH (5 μg/ml) and binding detected with anti-NS1 antibody. (b) Anti-NS1 antibody was coated on the plate and NS1 captured from Dengue virus (DV)2-infected cell supernatant. Purified CFH was then added and binding detected using an anti-CFH antibody. The experiments were performed three times independently. (c,d,e) Far Western blotting assays to detect interaction between CFH and NS1. CFH (c), supernatant (d) and recombinant NS1 (e) were immobilized using gel electrophoresis. Following addition of supernatant (c), recombinant NS1 (c) and CFH (d,e) interaction was tested by addition of the appropriate antibodies. (f) Immunoprecipitation assay to detect interaction between CFH and NS1. Immunoprecipitation using anti-NS1 antibody (lanes 2 and 3) demonstrated NS1 in the absence of CFH in DV2-infected cell supernatants only. In lane 1 purified CFH was loaded on to the gel as a positive control for CFH detection by the anti-CFH antibody.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Fig. S1. Far Western blot between Neisseria complement factor H binding protein (fHbp) and complement factor H (CFH). fHbp protein (a gift from Professor Chris Tang, University of Oxford, UK) was loaded into the gel, transferred to a membrane which was then incubated with human CFH (lane 1, 5 μg/ml; Quidel, San Diego, CA, USA) or buffer (lanes 2 and 3). The interaction between gelimmobilized fhbp and the added CFH was detected using an anti-human CFH antibody (lane 1, antibody obtained from Quidel). Gel-immobilized fHbp was visualized using an anti-fHbp (lane 2, antibody a gift from Professor Chris Tang, University of Oxford, UK). The anti-human CFH antibody did not recognize gel-immobilized fhbp in the absence of CFH (lane 3). Table S1. Primer sequences used in custom genotyping assays. Table S2. Relationship of complement gene polymorphisms between Dengue infection and controls.