# Prediction and Characterization of Missing Proteomic Data inDesulfovibrio vulgaris

Authors: Feng Li; Lei Nie; Gang Wu; Jianjun Qiao; Weiwen Zhang
Journal: Comparative and Functional Genomics
 (2011)

---

## Abstract

Proteomic datasets are often incomplete due to identification range and sensitivity issues. It becomes important to develop methodologies to estimate missing proteomic data, allowing better interpretation of proteomic datasets and metabolic mechanisms underlying complex biological systems. In this study, we applied an artificial neural network to approximate the relationships between cognate transcriptomic and proteomic datasets ofDesulfovibrio vulgaris, and to predict protein abundance for the proteins not experimentally detected, based on several relevant predictors, such as mRNA abundance, cellular role and triple codon counts. The results showed that the coefficients of determination for the trained neural network models ranged from 0.47 to 0.68, providing better modeling than several previous regression models. The validity of the trained neural network model was evaluated using biological information (i.e. operons). To seek understanding of mechanisms causing missing proteomic data, we used a multivariate logistic regression analysis and the result suggested that some key factors, such as protein instability index, aliphatic index, mRNA abundance, effective number of codons (Nc) and codon adaptation index (CAI) values may be ascribed to whether a given expressed protein can be detected. In addition, we demonstrated that biological interpretation can be improved by use of imputed proteomic datasets.


---

## Body

## 1. Introduction

Application of various experimental “omics” tools has enhanced our understanding of complex biological systems in the past decade [1–8]. However, due to technical limitations of these high throughput technologies and constraints on experimental design, most of these high-throughput “omics” datasets still suffer from significant missing values. The incomplete data available has impeded scientists from extracting correct information regarding cell metabolism. To address the issue, attempts to apply computational tools to impute missing values in various high throughput “omics” datasets have been made [9–11]. The most successful examples of such efforts were for missing transcriptomic data [12].Compared with transcriptomic dataset, proteomic datasets suffer even more serious missing data since its identification range and sensitivity are still not fully comparable with transcriptomic measurements [13]. In particular, when partial proteomics dataset was used in various integrated “omics” studies, the undetected proteins were simply assigned a “zero’’ value and were excluded from relationship modeling, which could bias any conclusion resulted from the integrated studies [14]. To overcome this problem, several methods have been adapted from the estimation of missing values in transcriptomic data to estimate the missing proteomics values by using the available measurements from other proteins, such as the k nearest neighbor and Bayesian Principal Component Analysis (BPCA) methods for imputing missing proteomic values in gel-based proteomics dataset [15, 16], and by integrating the GO (Gene Ontology) information into the proteomic data imputation [17]. Based on the assumption that there exists meaningful correlation between two types of datasets [14], in recent years we have developed the Zero-inflated Poisson (ZIP) linear regression model [18] and a stochastic Gradient Boosted Trees (GBT) nonlinear model [19] to uncover possible relationships between transcriptomic and proteomic data and to predict protein abundance for the proteins not experimentally detected.As part of a continuous effort to develop statistical tools for missing value imputation, we in the study applied an artificial neural network (ANN) approach to predict abundance for experimentally undetected proteins for several cognate transcriptomic and proteomic datasets of sulfate reducing bacteriumDesulfovibrio vulgaris [20–23]. The advantages of artificial neural networks include ability to tackle complex relationships and to deal with noisy data with missing values, and its capability of generalization [24, 25]. In the approach, a multilayer perceptron (MLP) with sigmoid activation function in the single hidden layer was used to capture the nonlinear relationship between transcriptomic data and protein abundance. The top ten relevant predictors for protein abundance identified in a previous study [19] were used as inputs to the MLP for training the output of the MLP to approximate the corresponding protein abundance of each gene. The trained MLP was validated using biological information. In addition, by using a multiple logistic regression analysis, we quantified the contribution of some sequence-based factors, such as mRNA abundance, protein instability index, gene length, “effective number of codons’’ (Nc), and CAI (Codon Adaptation Index) values to missing proteomic values. We also demonstrated that biological interpretation can be improved by using the imputed proteomic datasets.


## 2. Materials and Methods

### 2.1. Datasets

Two datasets fromDesulfovibrio vulgaris were analyzed in this study. The datasets covers the same species, although the experimental conditions are different and they were obtained by two independent studies [19–23]. The raw intensity values from both datasets were normalized using a quantile normalization in an R package (caret) available through the R project (http://www.r-project.org/). Absolute fluorescence intensity and peptide hits were used as abundance for transcriptomic and proteomic data, respectively [26]. These datasets, labeled Dataset 1 and Dataset 2, will be used throughout the article.


### 2.2. Genome Information

Gene annotated attributes, such as sequence length, protein length, molecular weight, GC content, triple codon counts, and cellular functional categories of all genes in theD. vulgaris genome were downloaded from the comprehensive microbial resource (CMR) of TIGR (http://cmr.tigr.org/) [22]. To define operons, we used the distance-only approach, a relatively low threshold, in this study to cover more of the possible operons. The genes that are transcribed in the same direction and are with intergenic region less than 15 base pairs were defined as one operon [19].


### 2.3. Construction of Relationship between Proteomics and mRNA through Artificial Neural Network

We trained an artificial neural network to learn a particular function by adjusting the values of the weights between neurons. It has been shown that a multilayer perceptron with a single hidden layer can approximate any continuous function to arbitrary precision by increasing the number of neurons in the hidden layer [27–29]. The weights connecting the input signals as well as the number of neurons can be determined by minimizing the sum of the squared difference between the network outputs and target outputs for each set of input signals. Nonlinear optimization methods in combination with a backpropagation algorithm [30] were implemented to find the optimal weights. The inputs to the network and the targets were usually rescaled to [−1, 1]. The outputs of the network were then rescaled back to approximate the targets in original scale. In this work, we trained a multilayer perceptron to learn the nonlinear relationship between transcriptomic data and proteomic data. The inputs were the transcriptomic data and the target was the protein abundance.A feed-forward network with a hidden layer can approximate any continuous function with arbitrary precision, which provides ANN the power to deal with nonlinear relationship, also open a door to overfitting if too many neurons are used in the hidden layer. To address the possible overfitting, we implemented a cross-validation to prevent overfitting. Through this method the data was partitioned inK equal parts where K-1 sets were used to train the model and the other set was used to calculate prediction errors [31]. This was repeated K times, hence K prediction errors values were computed one at every fold. An average and standard deviation can be extracted to select the most representative model for future prediction. Cross-validation was generally effective to prevent overfitting. Beside cross-validation to prevent overfitting, we also used some statistical tools to diagnose the selected model. Once the best model has been selected, based on cross-validation, the model was evaluated based on its coefficient of determination (R2). This coefficient represents the variation explained by the model. Furthermore, to alternatively assess the goodness of the model we studied the predictions of small sets of genes which are grouped based on operon information. In order to describe the variation within a data set, such as “molar abundance” of proteins within one operon, we computed the coefficient of variation (CV) for each set of proteins. The CV is defined as the ratio of the standard deviation and the mean of the “molar abundance” for a set of proteins [18, 19, 32] where the calculation of CV score is independent of the sample size. These coefficients of variations were computed for all operons and compared to a distribution of permuted CVs where permutation of genes was performed.


### 2.4. Analysis of the Factors Ascribed to Missing Protein Measurements

For DNA and protein sequence features, theBioPerl module was used to batch query and parse results from the ExPASy ProtParam program, returning “AliphaticIndex” and “InstabilityIndex”. CodonW was used to perform correspondence analysis of both codon and amino acid usages, returning the first four axes representing the major trends (AA_axis1–4 for amino acid usage, CR_axis1–4 for Relative Synonymous Codon Usage). The frequency of each base at the third codon position (synonymous site, designated as A3s, T3s, C3s, and G3s), as well as two additional measures on protein property such as Grand Average of Hydropathy (GRAVY) and AROMO, were also computed with CodonW. The Codon Adaptation Index (CAI) was computed with an in-house perl script using D. vulgaris ribosomal protein genes as reference for highly expressed genes. The calculation and analysis of these features were described in several of our previous works [33–35].A logistic regression model was used to assess how each of the features attributes to whether a given protein will be detected [36, 37]. Assume the features for a gene i are x1i,…,xki. In logistic regression, the logit, natural log of the odds, of the unknown probability, pi, that the protein is not detected is modeled as a linear function of x1i,…,xki:(1)log⁡it(pi)=ln⁡(pi1-pi)=β0+β1x1i+⋯+βkxki.
The model is equivalent to (2)pi=exp⁡(β0+β1x1i+⋯+βkxki)exp⁡(β0+β1x1i+⋯+βkxki)+1.
The unknown coefficients are usually estimated by maximum likelihood using a method common to all generalized linear models [35]. In the model, the parameter βj is the additive effect on the log of the odds for a unit change in the jth explanatory feature variable. That is, one unit change of feature j, that is, xji indicates the odds changes eβj times. Positive βj indicated that the probability that the protein will not be detected is increased with increasing the jth feature. By using a logistic regression analysis, the effects of relevant features on whether a protein will be detected were quantified and can be compared. And the fitted model can also be used to predict the chance that any other proteins will not be detected.


## 2.1. Datasets

Two datasets fromDesulfovibrio vulgaris were analyzed in this study. The datasets covers the same species, although the experimental conditions are different and they were obtained by two independent studies [19–23]. The raw intensity values from both datasets were normalized using a quantile normalization in an R package (caret) available through the R project (http://www.r-project.org/). Absolute fluorescence intensity and peptide hits were used as abundance for transcriptomic and proteomic data, respectively [26]. These datasets, labeled Dataset 1 and Dataset 2, will be used throughout the article.


## 2.2. Genome Information

Gene annotated attributes, such as sequence length, protein length, molecular weight, GC content, triple codon counts, and cellular functional categories of all genes in theD. vulgaris genome were downloaded from the comprehensive microbial resource (CMR) of TIGR (http://cmr.tigr.org/) [22]. To define operons, we used the distance-only approach, a relatively low threshold, in this study to cover more of the possible operons. The genes that are transcribed in the same direction and are with intergenic region less than 15 base pairs were defined as one operon [19].


## 2.3. Construction of Relationship between Proteomics and mRNA through Artificial Neural Network

We trained an artificial neural network to learn a particular function by adjusting the values of the weights between neurons. It has been shown that a multilayer perceptron with a single hidden layer can approximate any continuous function to arbitrary precision by increasing the number of neurons in the hidden layer [27–29]. The weights connecting the input signals as well as the number of neurons can be determined by minimizing the sum of the squared difference between the network outputs and target outputs for each set of input signals. Nonlinear optimization methods in combination with a backpropagation algorithm [30] were implemented to find the optimal weights. The inputs to the network and the targets were usually rescaled to [−1, 1]. The outputs of the network were then rescaled back to approximate the targets in original scale. In this work, we trained a multilayer perceptron to learn the nonlinear relationship between transcriptomic data and proteomic data. The inputs were the transcriptomic data and the target was the protein abundance.A feed-forward network with a hidden layer can approximate any continuous function with arbitrary precision, which provides ANN the power to deal with nonlinear relationship, also open a door to overfitting if too many neurons are used in the hidden layer. To address the possible overfitting, we implemented a cross-validation to prevent overfitting. Through this method the data was partitioned inK equal parts where K-1 sets were used to train the model and the other set was used to calculate prediction errors [31]. This was repeated K times, hence K prediction errors values were computed one at every fold. An average and standard deviation can be extracted to select the most representative model for future prediction. Cross-validation was generally effective to prevent overfitting. Beside cross-validation to prevent overfitting, we also used some statistical tools to diagnose the selected model. Once the best model has been selected, based on cross-validation, the model was evaluated based on its coefficient of determination (R2). This coefficient represents the variation explained by the model. Furthermore, to alternatively assess the goodness of the model we studied the predictions of small sets of genes which are grouped based on operon information. In order to describe the variation within a data set, such as “molar abundance” of proteins within one operon, we computed the coefficient of variation (CV) for each set of proteins. The CV is defined as the ratio of the standard deviation and the mean of the “molar abundance” for a set of proteins [18, 19, 32] where the calculation of CV score is independent of the sample size. These coefficients of variations were computed for all operons and compared to a distribution of permuted CVs where permutation of genes was performed.


## 2.4. Analysis of the Factors Ascribed to Missing Protein Measurements

For DNA and protein sequence features, theBioPerl module was used to batch query and parse results from the ExPASy ProtParam program, returning “AliphaticIndex” and “InstabilityIndex”. CodonW was used to perform correspondence analysis of both codon and amino acid usages, returning the first four axes representing the major trends (AA_axis1–4 for amino acid usage, CR_axis1–4 for Relative Synonymous Codon Usage). The frequency of each base at the third codon position (synonymous site, designated as A3s, T3s, C3s, and G3s), as well as two additional measures on protein property such as Grand Average of Hydropathy (GRAVY) and AROMO, were also computed with CodonW. The Codon Adaptation Index (CAI) was computed with an in-house perl script using D. vulgaris ribosomal protein genes as reference for highly expressed genes. The calculation and analysis of these features were described in several of our previous works [33–35].A logistic regression model was used to assess how each of the features attributes to whether a given protein will be detected [36, 37]. Assume the features for a gene i are x1i,…,xki. In logistic regression, the logit, natural log of the odds, of the unknown probability, pi, that the protein is not detected is modeled as a linear function of x1i,…,xki:(1)log⁡it(pi)=ln⁡(pi1-pi)=β0+β1x1i+⋯+βkxki.
The model is equivalent to (2)pi=exp⁡(β0+β1x1i+⋯+βkxki)exp⁡(β0+β1x1i+⋯+βkxki)+1.
The unknown coefficients are usually estimated by maximum likelihood using a method common to all generalized linear models [35]. In the model, the parameter βj is the additive effect on the log of the odds for a unit change in the jth explanatory feature variable. That is, one unit change of feature j, that is, xji indicates the odds changes eβj times. Positive βj indicated that the probability that the protein will not be detected is increased with increasing the jth feature. By using a logistic regression analysis, the effects of relevant features on whether a protein will be detected were quantified and can be compared. And the fitted model can also be used to predict the chance that any other proteins will not be detected.


## 3. Results and Discussion

### 3.1. Determine the Variable Importance and Partial Variable Dependence

In a previous study, we used the gradient boosted trees (GBT) model to measure the contribution of each of the 70 variables input variables to the RNA-protein correlation, and the 10 top-ranked variables were identified to be mRNA expression level, cellular role, and several triple codon usages for theD. vulgaris datasets we used (See Table  1 of Supplmentary Material available online at doi: 10.1155/2011/780973) [19]. In this study, we used the same top ten predictors as inputs to the artificial neural network approach. The datasets used are consistent with previous work for comparison. The results showed that for all datasets, we observe a smooth relationship between protein abundance and the abundance of corresponding mRNA controlling the other predictors at their mean values for each growth condition. The partial dependence between protein and mRNA abundance tends to be smoother when the abundance of mRNAs are in the high range (Figure 1). Figure 2 shows the scatter plots of predicated protein abundance versus mRNA for proteins detected and undetected for datasets 1 and dataset 2, respectively. For both datasets, under each growth condition, there are much fewer genes with very high mRNA or protein abundance. Those proteins undetected are mostly with very low predicted protein abundance and observed mRNA abundance as shown in Figures 2(b) and 2(d) for dataset 1 and dataset 2, respectively. The proteins undetected have much smaller variance in predicted values especially for dataset 2 as shown in Figure 2(d). In comparison, those detected proteins have a larger range of protein abundance and are of very high variability.Partial dependency plots. Predicted protein abundance for given values of mRNA controlling the other predictors at their mean values. (a) Dataset 1; (b) dataset 2.
(a)(b)Prediction plot for proteins. Scatter plots of predicted protein abundance versus mRNA for proteins detected and undetected. (a) Plot for genes with detected proteins for dataset 1; (b) plot for genes with undetected proteins for dataset 1; (c) plot for genes with detected proteins for dataset 2; (d) plot for genes with undetected proteins for dataset 2.
(a)(b)(c)(d)In the previous work of using GBT model to uncover the nonlinear relationship between transcriptomic and proteomic datasets, an interesting pattern, a plateau effect, was observed in the partial dependence plot of mRNA abundance versus protein abundance for high values of mRNA. This is potentially due to the fact that protein abundance data is sparse, with high variance, where the tree models do not generate splits among the predictors [19]. In the neural network model, the output of an MLP is an inherent smooth function of the network inputs, hence a relatively smooth prediction is always guaranteed. However, we can still see the similar pattern that the curve becomes flat when the variance for the mRNA abundance is high. For example, in the region of high mRNA values there are a small number of genes/peptides whose protein values ranges from (0, 40) for Dataset 1 and (0, 500) for Dataset 2. This could reflect problems with the sensitivity of the current proteomic technologies for high abundant proteins.


### 3.2. Construction of the Nonlinear Artificial Neural Network Model

A feed-forward network with a single hidden layer was trained for each condition for both datasets with mRNA abundance and other predictors as inputs and protein abundance as target. The neural networks were trained using the functions in the neural network toolbox provided in standard software Matlab (7.0.1). The activation function for the neurons in the hidden layer was the hyperbolic tangent function, which is given by(3)f(v)=1-exp⁡⁡(-v)1+exp⁡⁡(-v).
The activation function of the output neuron is a pure linear function, which is the weighted sum of the outputs from the neurons of the hidden layer. The optimal number of neurons in the hidden layer and the optimal weights were determined using 10-fold cross-validation for training a neural network. For each dataset under each experimental condition, the networks with the same inputs but with different number of neurons (from 5 to 12) in the hidden layer were each trained using cross-validation to learn the relationship between the mRNA features and protein abundance. The optimal number of neurons in the layer was determined based on the one generating the minimum averaged validation error.Both transcriptomic and proteomic data were log-transformed. For dataset 1, the genes with no missing value for three out of four replicates were selected for model training. The numbers of genes used for training were 354, 312, and 308 for the FL, LL, and LS conditions, respectively. For each growth condition, the inputs to the multilayer perceptron were the top ten predictors identified previously [19] (See Table  1 at Supplementary Material). The optimal network had 6, 6, and 9 neurons in the hidden layer for FL, LL, and LS conditions, respectively. Coefficients of determinations R2 were 0.47, 0.52, and 0.51 for FL, LL, LS conditions, respectively. For dataset 2, a total of 986, 1001, and 986 genes with no missing protein measurements were used for model training for CT0, CT120, and ST120 conditions, respectively. For each of the conditions, the inputs to the multilayer perceptron were the top ten predictors identified previously [19] (See Table  1 at Supplementary material). The optimal network had 8, 7, and 6 neurons in the hidden layer for CT0, CT120, and ST120 conditions, respectively. Coefficients of determination R2 were 0.65, 0.62, and 0.68 for CT0, CT120, and ST120 conditions, respectively.


### 3.3. Validation of the Prediction through the Nonlinear Artificial Neural Network Model

The validation of the model prediction was conducted through the calculation of coefficient of variation (CV) across conditions for every operon ofD. vulgaris as described before [18, 19]. The genes in the operons are considered to have less dispersion than any random set of genes in terms of their gene expression and/or protein abundance, because of their intrinsic biological relationship. With relatively low threshold for operon identification as described in Section 2, a total of 609 operons were identified in D. vulgaris, which ranged from 2 to 13 genes per operon. To compare these CV values, we have used a permutation test in the following manner: a CV is computed from the predicted protein values for a permuted number of genes. This step is repeated 2000 times through resampling of genes without replacement. The number of genes randomly selected is equal to the size of the gene group to be compared to. For dataset 1, 82% to 89% of the operon groups had smaller CV values than the permuted groups, and for dataset 2, 79% to 89% of the operon groups had smaller CV values than the permuted, demonstrating in general good performance of the model prediction. The validation results are summarized in Table 1.Table 1
Model validation: correlated expressions of proteins in selective operons.
Dataset 1Dataset 2FLLLLSCT0CT120ST120OperonCVPCVCVPCVCVPCVCVPCVCVPCVCVPCV10.360.7070.60.7970.3010.9011.63*1.2480.8141.1861.2691.28621.6*0.6381.564*0.7251.48*0.6310.7731.1020.9681.0440.7021.12930.4990.6591.04*0.7660.5630.8661.295*1.1341.568*1.1061.324*1.16940.2160.6480.2390.7260.3230.8130.3331.1060.4011.050.4531.13550.4180.6460.5590.7310.450.7940.4811.1050.3981.050.4661.12260.180.6470.4540.7450.3530.8090.4951.1160.4031.0420.4731.12470.2970.6510.1620.7290.3120.80.4081.1010.3431.0420.691.11880.3070.6930.2920.7720.280.9080.7671.2040.5021.1470.571.23690.3890.6760.3520.7640.8510.860.8841.1710.6311.0980.5581.197100.6210.6240.4450.7650.740.8590.8891.1660.671.0961.0871.19110.3040.5950.2310.7040.30.7730.781.1530.8161.0820.5671.182120.5990.8080.8310.9080.7291.0761.2011.5351.1741.3661.563*1.528130.3530.6460.4350.7230.5310.8260.3661.1180.5321.060.3691.124140.3590.6650.4070.730.2670.7520.6531.1570.8251.0870.5971.18150.3810.7170.7070.8010.780.9220.6911.2460.8381.1660.6741.293160.2990.7490.5380.8330.3070.9890.6131.3420.4561.2521.0591.372170.3210.6490.2190.7260.2220.8060.6571.0590.4920.9430.5291.046180.4730.6450.4860.7230.6280.7670.4741.1640.5391.0910.5781.185190.5920.7410.4790.8240.6310.9821.68*1.3171.0941.2481.403*1.342200.767*0.6210.5860.6950.5610.7520.2671.0370.2880.9640.4921.053210.3890.5480.4160.6930.4720.7571.0341.0370.8130.9780.8231.061220.2560.6540.4330.7210.490.80.5321.0950.6021.0530.5691.123230.6310.6940.7070.7750.6320.8770.861.1960.6771.1250.9251.242240.3760.6550.6480.6550.4060.6991.23*1.0990.9271.0311.138*1.123250.869*0.6741*0.7541.194*0.852.337*1.1562.119*1.0942.303*1.185260.657*0.6180.6830.7360.7470.8560.6961.1510.521.1140.71.176270.3780.850.4520.9380.4321.1040.5881.6310.6671.4450.7711.629281.013*0.670.963*0.7641.526*0.9131.2551.2591.375*1.1811.332*1.282CV is computed by dividing SD by the mean of the prediction values for protein abundance for a specific set of genes (group). The protein prediction values were normalized by molecular weight before CV calculation. PCV is the mean of CV values computed through permutation test for selected operons.*CV values greater than PCV.To further demonstrate how the current measurement of dispersion (CV) compares with the distribution of the permuted coefficients of variation, we also calculated the percentile score for all coefficients of variations, based on the biological knowledge that genes from the same operon should be less dispersed than permuted sets of genes and the percentile scores are thus expected to be low. The results showed that even though the calculated CV for most groups was less than the mean CV value for permuted sets of genes, the percentile scores shows that only a small collection of these groups fall within a percentile less than 0.20. The percentage of groups with percentile scores less than 0.2 is 43% to 57%, and 57% to 61% for dataset 1 and dataset 2, respectively.


### 3.4. Comparison of Artificial Neural Network to Previous Methods (GBT and ZIP)

Similar to the GBT method, the artificial neural network method aims to identify an unknown and potentially nonlinear relationship between proteomics and mRNA abundance data, which distinguish them from the ZIP method, which carries zero-inflated linear assumptions [18]. Another key difference is that GBT and neural network methods utilized additional sequence based information, such as cellular role, gene length, amino acid usage, and codon usage known to be important for translation efficiency to model the relationship between transcriptomic and proteomic dataset, while the ZIP method establishes the relationship between proteomic and mRNA solely based on experimental transcriptomic and proteomic dataset.Some improvements were observed in terms of prediction validation by the neural network model in comparison to the GBT model. The percentage of operon groups with CV smaller than permuted CV was 82% to 89% for dataset 1 in comparison to 75% to 79% for the GBT model. The percentage of operon groups with CV smaller than permuted CV was 79% to 89% for dataset 2 by the neural network model in comparison to 75% to 82% by the GBT model. For dataset 1, the percentage of groups with percentile scores less than 0.2 is 43% to 57% in the neural network model in comparison to 36% to 50% in the GBT model. For dataset 2, the percentage of groups with percentile scores less than 0.2 is 57% to 61% in the neural network model in comparison to 32% to 54% in the GBT model. The coefficient of determination was 0.47 to 0.58 for dataset 1 and 0.62 to 0.68 for dataset 2 in the neural network model in comparison to 0.39 to 0.58 for both datasets in the GBT model, demonstrating the artificial neural network model could be a good alternative in analyzing the missing proteomic values.


### 3.5. Analysis of Factors Ascribed to Missing Protein Measurements

Through the neural network model, the majority of the proteins that were experimentally undetected were predicted to have protein abundance with greater than 1.0 peptide hit. This thus tempted for us to look into the possible factors which may be ascribed to the missing protein values. To do so, we collected some of the protein or mRNA/DNA sequence-based features which may affect protein detectability for use along with experimental transcriptomic data in a logistic regression procedure to determine their relative contribution to detectability of proteins in theD. vulagris genome. The features/measurements we used were from the following several categories: (i) features related to gene expressivity: experimental mRNA abundance data, Nc and CAI values [33, 34], codon usage such as CR_axis1-CR_axis4, T3s, C3s, A3s, G3s, GC3s, and GC content [35]; (ii) features related to protein and RNA stabilities: protein instability index and AliphaticIndex, and minimal free energy (MFE) of RNA [35]; (iii) features related to protein solubility: protein solubility indices such as Grand Average of Hydropathy (GRAVY) [33]; (iv) amino acid usage such as AA_axis1-AA_axis4, and others such as gene length [35]. The logistic regression procedure in the SAS software package was used to model the probability of whether a protein will be experimentally detected given the features and measurements for the gene/protein. The stepwise selection method was used to select those most influential factors for indicating missing proteomic value. Those highly correlated features were removed automatically. The 0.05 significance threshold was used for the stepwise selection. The mRNA abundance and gene length values were both log2 transformed before being used in the procedure. The protein stability index, aliphaticIndex, and Nc values were divided by 10 and the CAI, N_MFE, GRAVY, AA_axis1-AA_axis4, and CR_axis1-CR_axis4 were all multiplied by 10 so that they scaled to 1. This standardization will not change the results, but make illustration of the results easier.Tables2 and 3 present the odds ratio of each selected factor for dataset 1 and dataset 2, respectively. The results showed that some of the features are highly correlated with probability of missing proteomic values. For dataset 1, seven factors were indentified consistently by the regression model. Among the features that are positively correlated with probability of missing proteomic values (i.e., when the features/measurement increase), protein instability index which provides an estimate of the stability of a given protein in a test tube [38] was found as the most significant factor for all three growth conditions, and aliphaticIndex which is regarded as a positive factor for the increase of thermostability of globular proteins [39], was found significant for growth conditions LL and FL, but not for condition LS. In addition, AA_axis2 and CR_axis3 correlated most strongly with GRAVY and G3s, respectively, and they were also found significant in all three growth conditions, while the N_MFE was only significant for LS condition (Table 2). For dataset 2, the features that are kept in the regression model are more the same for the three growth conditions in dataset 1, suggesting that these factors could be universal across experimental platforms. Six features that are significant for dataset 1 are also found to be critical for all the conditions and affect the probability in the same direction for dataset 2 include mRNA abundance, protein instability index, gene length, Nc, CAI, and AA_axis2 corresponding to GRAVY (Table 3). In general, the analysis suggested that those proteins with smaller mRNA expression values, smaller Nc value, higher protein instability index, or high value of aliphaticIndex, have a greater probability to have missing proteomic values in the D. vulgaris datasets. Based on this preliminary analysis, it may be possible to apply a logistic regression model to predict whether a given protein will or will not be detected by choosing a proper threshold for the predicted probability.Table 2
Estimates of odds ratios of missing proteomic values for dataset 1.
Dataset 1LLFLLSeffectPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitslog2 of mRNA0.4580.4050.5170.4420.3920.4990.4050.3480.471Protein instability index/101.6281.3791.9221.4931.2821.7391.5391.2971.825log2 of gene length0.4770.3730.6080.5170.4150.6440.4190.3240.543Nc/100.4650.2710.7980.4460.2720.7320.3140.1740.566CAI*100.2950.2120.4100.2750.2040.3700.2220.1550.319AA_axis2*101.3681.1931.5691.3261.1701.5021.3231.1661.501CR_axis3*101.2641.0461.5281.3091.1121.5411.3271.0921.613N_MFE*10——————1.8801.1783.000AliphaticIndex/101.1401.0131.2821.1601.0401.294———Table 3
Estimates of odds ratio of missing proteomic values for dataset 2.
Dataset 2CT0CT120ST120effectPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitslog2 of mRNA0.4070.3700.4480.4510.4120.4940.4990.4590.542Protein instability index/101.3101.1611.4781.2831.1401.4441.2641.1241.421log2 of gene length0.2610.2160.3150.2870.2400.3450.2940.2470.352Nc/100.5120.3360.7800.5180.3430.7830.5520.3740.815GRAVY*100.9140.8620.9700.9140.8630.9680.9230.8730.977CAI*100.5990.4210.8510.5590.3980.7860.5360.3840.749AA_axis1*100.6860.6270.7500.6810.6240.7440.6860.6290.749AA_axis2*101.8111.5582.1061.7211.4871.9901.6741.4521.931AA_axis4*100.8270.7470.9160.8370.7570.9250.8460.7670.934CR_axis1*100.7710.6470.9190.7800.6570.9250.7810.6600.923CR_axis4*100.7160.6350.8080.7060.6270.7950.7250.6450.814


### 3.6. Data Interpretation Using Imputed Datasets

A total of 308–354 proteins were experimentally identified for dataset 1, and 986–1001 proteins were experimentally identified for dataset 2, which are slightly over 10–30% of the whole proteome inD. vulgaris. Through the ANN modeling and imputation by integrating transcriptomic and proteomic data, we were able to assign abundance values to approximately 33–3400 proteins (dependent on growth conditions), approximately 96–98% of the whole proteomes inD. vulgaris. Using the imputed datasets, we first analyzed the proteins involved in energy metabolism when D. vulgaris was grown in three cultivation conditions (LL: Lactate-based growth in Exponential phase; LS: Lactate-based growth at Stationary phase; FL: Formate-based growth at Exponential phase). As shown in Table 4, even for proteins involved in energy metabolism such as ATP synthase which are typically expressed at relatively high levels, significant missing proteomic data was still obvious in the experimental datasets, only two out of nine putative ATP synthase proteins were experimentally detected under three growth conditions. By analysis of the imputed protein abundance for all putative ATP synthase proteins, a clear trend of downregulation was observed for LS and FL condition when compared with LL condition, which is consistent with the previous data suggesting lactate is a favorable electron donor for D. vulgaris [22]. When comparing lactate-based growth with formate-based at the exponential phase, we found that enzymes involved for formate utilization were predicted to be upregulated during growth on formate (Table 4). It has been reported that c-type cytochromes are highly expressed in sulphate-reducing bacteria [20], however, only low amounts of c-type cytochrome proteins have been identified from D. vulgaris by the proteomic (Table 4), probably due to the fact that c-type cytochromes undergo a complex posttranslational maturation process involving covalent attachment of heme groups [40]. Through the ANN model, significant expression of various cytochrome c proteins was predicted in D. vulgaris grown under three conditions. In general, cytochrome c proteins were predicted to have higher expression at the exponential phase in both lactate- and formate-based media than the stationary phase, consistent with their important roles in electron transfer during fast growth [20]. Two cytoplasmic hydrogenases (Ech and Coo) were previously assigned the putative function of generating hydrogen in the cytoplasm [20]. In agreement with previous prediction from the ZIP model [18], the ANN model predicted that all subunits of Ech hydrogenase (EchACDE) were upregulated, while all subunits of Coo hydrogenase (CooHLUX) were down-regulated on formate-based growth, suggesting their differential roles in H2 generation when grown with different carbon sources. We also used the imputed data to analysis the stress responses when the D. vulgaris cells were treated under salt stress. The analysis allowed identification of some proteins which were not experimentally detected, however, were predicted to be upregulated under the salt treatment condition (ST120) as compared with control (CT120) (Table 5). Among them were transporters and various multidrug resistance proteins and some signal transduction proteins. This list of proteins can serve as putative candidates for further experimental validation. In general, our analysis indicates that biological interpretation may benefit from imputing missing data using computational methods by integrating temporal transcriptomic and proteomic data. The full list of all protein-abundance predictions for all proteins has been provided in the Table  2 of the supplementary material.Table 4
Regulation of proteins involved in energy metabolism for Dataset 1*.
Gene_IDProtein abundance-Predicted-LL conditionRNA abundance-Measured-LL conditionProtein abundance-Measured-LL conditionProtein abundance-Predicted-LS conditionRNA abundance-Measured-LS conditionProtein abundance-Measured-LS conditionProtein abundance-Predicted-FL conditionRNA abundance-Measured-FL conditionProtein abundance-Measured-FL conditionDescriptionDVU091810.082315.20x2.46416.75x5.711668.90xATP synthase F0, A subunit (atpB)DVU07796.552505.40x3.30497.75x5.081397.10xATP synthase F0, B subunitDVU07806.623156.20x3.881050.10x5.812352.70xATP synthase F0, B subunitDVU09177.606086.50x3.851234.90x10.855944.40xATP synthase F0, C subunit (atpE)DVU09203.801590.50x3.59947.88x5.751198.20xATP synthase protein I (atpI)DVU077739.242535.0020.008.40369.958.0017.76735.6817.33ATP synthase, F1 alpha subunit (atpA)DVU077550.112228.6049.3323.76711.5330.3339.801198.4038.33ATP synthase, F1 beta subunit (atpD)DVU07784.251974.50x2.71374.22x3.21714.57xATP synthase, F1 delta subunit (atpH)DVU07744.852043.70x2.32336.55x4.141113.40xATP synthase, F1 epsilon subunit (atpC)DVU07022.62288.72x2.20372.82x2.05417.681.00Cytochrome c family proteinDVU09221.6966.60x1.85149.38x1.49149.52xCytochrome c family proteinDVU12883.441172.20x2.06298.10x3.12826.90xCytochrome c family proteinDVU248310.24456.23x3.24284.22x7.36490.33xCytochrome c family proteinDVU24845.06575.50x2.39277.35x4.63665.95xCytochrome c family proteinDVU27914.77683.65x2.61530.00x3.32667.05xCytochrome c family proteinDVU310715.101009.302.332.98177.151.6713.041026.303.00Cytochrome c family proteinDVU31445.58297.85x3.58142.45x8.24342.52xCytochrome c family proteinDVU062530.201822.80x14.80764.85x16.612467.50xCytochrome c nitrite reductase, catalytic subunit NfrADVU18154.91439.87x2.98324.93x7.12833.88xCytochrome c oxidase, subunit IDVU18124.42304.57x2.72215.70x4.63335.70xCytochrome c oxidase, subunit IIDVU18142.53397.35x2.31283.90x3.80898.67xCytochrome c oxidase, subunit IIIDVU28093.11447.42x2.32395.15x2.10601.50xCytochrome c3DVU317113.243937.003.335.24966.25x11.825719.201.00Cytochrome c3DVU25241.77114.95x2.08135.07x1.61170.70xCytochrome c3,DVU30413.03894.60x2.62713.10x2.73684.27xCytochrome c553DVU18176.342280.0011.677.261858.203.339.615400.2012.00Cytochrome c-553 (cyf)DVU10482.31255.07x2.14243.58x1.53185.47xCytochrome c-type biogenesis protein CcmB (ccmB)DVU10472.87435.45x2.06284.87x2.23259.25xCytochrome c-type biogenesis protein CcmC (ccmC)DVU10513.30915.92x3.59769.45x4.19732.07xCytochrome c-type biogenesis protein CcmE (ccmE)DVU10507.42581.02x4.03408.20x6.27463.25xCytochrome c-type biogenesis protein CcmF (ccmF)DVU327135.798217.80x12.152295.10x11.386954.00xCytochrome d ubiquinol oxidase, subunit I (cydA)DVU327018.093877.70x5.501145.70x13.903897.10xCytochrome d ubiquinol oxidase, subunit II (cydB)DVU04341.28177.22x2.79168.72x5.13468.00xEch hydrogenase, subunit EchADVU04332.10217.82x2.34166.95x2.84422.28xEch hydrogenase, subunit EchBDVU04322.11222.90x2.24133.13x3.12595.60xEch hydrogenase, subunit EchCDVU04312.20303.27x2.03215.73x3.30778.70xEch hydrogenase, subunit EchDDVU04302.68259.00x2.52199.20x4.61685.78xEch hydrogenase, subunit EchEDVU04292.24239.03x2.07201.57x2.05405.30xEch hydrogenase, subunit EchFDVU28240.90124.20x3.54106.78x9.30102.10xFormate acetyltransferaseDVU22721.0597.08x3.69113.331.0016.21160.88xFormate acetyltransferaseDVU05782.11154.12x2.23217.48x2.03247.52xFormate dehydrogenase accessory protein FdhDDVU05772.56237.93x2.29148.02x1.97513.40xFormate dehydrogenase formation protein FdhEDVU28103.39423.00x2.27312.02x2.47510.12xFormate dehydrogenase formation protein FdhEDVU05880.9235.70x2.6535.25x1.2679.03xFormate dehydrogenase, beta subunitDVU24812.39227.43x2.37304.92x2.35399.82xFormate dehydrogenase, beta subunitDVU28113.22335.33x2.52181.45x3.13435.95xFormate dehydrogenase, beta subunitDVU22882.68781.45x2.24412.85x1.39139.17xHydrogenase, CooL subunitDVU228616.67308.33x3.29138.72x6.0881.55xHydrogenase, CooM subunitDVU22902.84703.223.002.10316.453.00xHydrogenase, CooU subunitDVU22893.92823.18x2.14553.33x2.36319.80xHydrogenase, CooX subunit*Protein abundance is in peptide counts; RNA abundance is in raw fluorescence measurements in microarrays. “x’’ indicated proteins are undetected experimentally.Table 5
Up-regulated proteins under stress condition for Dataset 2*.
Gene_IDProtein abundance-Predicted-CT120RNA abundance-Measured-CT120Protein abundance-Measured-CT120Protein abundance-Predicted-ST120RNA abundance-Measured-ST120Protein abundance-Measured-ST120DescriptionAmino acid biosynthesisDVU02855.003005.10x7.434051.50xImidazole glycerol phosphate synthase (hisH)DVU166521.8220264.00x50.5530766.00x3-dehydroquinate dehydratase, type II (aroQ)Biosynthesis of cofactors, prosthetic groups, and carriersDVU01575.457559.00x6.315092.30xThiamin-monophosphate kinase (thiL)DVU09315.411730.40x6.874401.40xPhosphomethylpyrimidine kinase (thiD)DVU12005.825439.20x6.755079.60xRiboflavin synthase, alpha subunit (ribE)DVU27492.851427.30x3.59856.35xPrecorrin-6Y C5,15-methyltransferase (decarboxylating) (cobL)DVU330714.6520224.00x26.4230616.00x3-octaprenyl-4-hydroxybenzoate carboxy-lyase (ubiX)Cell envelopeDVU02942.441034.50x3.15726.82xGlycosyl transferase, group 2 family proteinDVU03082.971552.00x7.936986.70xMembrane proteinDVU03288.595628.30x9.125042.40xGlycosyl transferase, group 1 family proteinDVU040726.6520237.00x38.7032568.00xRare lipoprotein A family proteinDVU15342.481792.60x3.562183.30xMembrane proteinDVU20532.982406.60x3.852507.60xMembrane proteinDVU21022.91628.30x3.88843.26xOuter membrane protein, OMP85 familyDVU235620.546021.30x36.267883.90xMembrane proteinDVU26723.303253.30x4.475875.90xMembrane proteinDVU27255.204925.80x13.978894.70xMembrane proteinCellular processesDVU030228.3522358.00x37.2121180.00xChemotaxis protein CheXDVU20692.18811.21x3.72658.60xDNA processing protein DprADVU28164.82539.18x8.651502.80xMultidrug resistance proteinDVU28934.599413.30x6.789241.60xFlagellar basal-body rod proteinDVU33263.13782.67x3.682732.60xMultidrug resistance protein, Smr familyDVUA00025.532257.50x12.966384.00xParA family proteinDNA metabolismDVU08783.201967.30x7.1311149.00xdnaK suppressor proteinDVU13532.261516.00x10.31868.27xDNA polymerase III, alpha subunit (dnaE)DVU17031.93906.31x2.69628.53xType I restriction-modification enzyme, R subunitEnergy metabolismDVU00804.422476.60x6.455175.70xFumarate hydratase, class II (fumC)DVU05316.57328.82x11.021157.90xHmc operon protein 6DVU21373.64389.78x4.45338.82xSuccinyl-CoA synthase, beta/alpha subunits (sucCD)Other categoriesDVU01973.46933.13x4.40285.54xPhage portal protein, lambda familyDVU14904.72351.64x5.70372.18xTail tape measure proteinRegulatory functionsDVU00633.743823.20x4.796404.70xTranscriptional regulator, MarR familyDVU13317.546812.50x8.3311619.00xTranscriptional regulator, LysR familyDVU16457.227075.60x14.0323482.00xTranscriptional regulator, ArsR familyDVU30622.341205.40x2.92570.26xSensor histidine kinase/response regulatorDVU309512.3017923.00x14.8623528.00xTranscriptional regulator, Fur familyDVU07214.201085.90x7.19513.28xSensory box histidine kinaseDVU08345.732795.50x12.1713860.00xRibonuclease HII (rnhB)DVU107513.6330772.00x15.4527576.00xRibonuclease P protein component (rnpA)Transport and binding proteinsDVU05932.58992.23x3.281151.30xL-lysine exporterDVU10683.92301.98x5.76270.52xBranched-chain amino acid ABC transporter, permease proteinsDVU23406.0011276.00x7.0011625.00xAmino acid ABC transporter, permease proteinDVU25723.192831.50x5.777937.50xFerrous iron transport protein ADVU25744.625005.50x8.858857.00xFerrous ion transport proteinDVU274417.961201.00x34.423053.50xHigh-affinity branched-chain amino acid ABC transporterDVUA01363.29673.57x3.89346.83xZinc transporter (zupT)Unknown functionDVU01558.3516519.00x12.1518298.00xType I phosphodiesterase/nucleotide pyrophosphatase family proteinsDVU093810.1215409.00x11.7613138.00xIsoamylase N-terminal domain proteinDVU10237.26123.66x8.15144.26xRhomboid family proteinDVU12025.393480.30x5.925395.40xCytidine/deoxycytidylate deaminase family proteinDVU15585.77189.07x6.47234.72xCysteine-rich domain/iron-sulfur cluster-binding domain proteinsDVU17303.161288.70x3.744503.10xDNA-binding proteinDVU17456.55319.33x7.58182.38xDNA-binding proteinDVU17475.08258.33x12.80230.72xATPase, histidine kinase-, DNA gyrase B-, and HSP90-like domainsproteinDVU20364.271086.70x5.07350.61xHelix-turn-helix protein, CopG familyDVUA00393.75269.47x4.27261.42xChain length determinant family protein*Protein abundance is in peptide counts; RNA abundance is in raw fluorescence measurements in microarrays. “x’’ indicated proteins are undetected experimentally.In conclusion, missing values in large-scale proteomic analysis is a frequent problem and in some cases it has created difficulty for accurate interpretation of proteomic data for complex biological systems. However, in spite of its obvious importance, the methodologies to deal with missing proteomic data are still lagging well behind the methods for transcriptomic analysis. In this study, we applied an artificial neural network approach to approximate the relationships between cognate transcriptomic and proteomic datasets ofD. vulgaris, and to predict protein abundance for the proteins not experimentally detected based on relevant predictors, such as mRNA abundance, cellular role, and triple codon counts. The results showed that the coefficients of determination for the trained neural network models ranged from 0.47 to 0.68, providing better description than several previous models [18, 19]. In an attempt to identify possible causes for missing proteomic values, a logistic regression analysis of various experimental measurements and sequence-based features was performed for the experimentally undetected proteins but with high predicted protein abundance values. Although the results are in general consistent with previous knowledge that proteins with low possible expressivity and high instability will have higher chance to be missed, in the study, we for the first time quantified the odds ratio of each of those possible factors in the D. vulgaris datasets. In addition, in a preliminary analysis we demonstrated that biological interpretation could be improved by using the imputed proteomic datasets. Finally, although the model was well validated statistically, caution needs to be exercised when interpreting experimental data based on predicted protein expression values because the predicted abundance values are constrained by the quality of experimental proteomic data used as input, and a dearth of large-scale quantitative predictors which can be included in the model. Nevertheless, the initial success in applying the ANN method is encouraging, and the model could serve as a basis for developing more sophisticated models to further improve the biology interpretation.


## 3.1. Determine the Variable Importance and Partial Variable Dependence

In a previous study, we used the gradient boosted trees (GBT) model to measure the contribution of each of the 70 variables input variables to the RNA-protein correlation, and the 10 top-ranked variables were identified to be mRNA expression level, cellular role, and several triple codon usages for theD. vulgaris datasets we used (See Table  1 of Supplmentary Material available online at doi: 10.1155/2011/780973) [19]. In this study, we used the same top ten predictors as inputs to the artificial neural network approach. The datasets used are consistent with previous work for comparison. The results showed that for all datasets, we observe a smooth relationship between protein abundance and the abundance of corresponding mRNA controlling the other predictors at their mean values for each growth condition. The partial dependence between protein and mRNA abundance tends to be smoother when the abundance of mRNAs are in the high range (Figure 1). Figure 2 shows the scatter plots of predicated protein abundance versus mRNA for proteins detected and undetected for datasets 1 and dataset 2, respectively. For both datasets, under each growth condition, there are much fewer genes with very high mRNA or protein abundance. Those proteins undetected are mostly with very low predicted protein abundance and observed mRNA abundance as shown in Figures 2(b) and 2(d) for dataset 1 and dataset 2, respectively. The proteins undetected have much smaller variance in predicted values especially for dataset 2 as shown in Figure 2(d). In comparison, those detected proteins have a larger range of protein abundance and are of very high variability.Partial dependency plots. Predicted protein abundance for given values of mRNA controlling the other predictors at their mean values. (a) Dataset 1; (b) dataset 2.
(a)(b)Prediction plot for proteins. Scatter plots of predicted protein abundance versus mRNA for proteins detected and undetected. (a) Plot for genes with detected proteins for dataset 1; (b) plot for genes with undetected proteins for dataset 1; (c) plot for genes with detected proteins for dataset 2; (d) plot for genes with undetected proteins for dataset 2.
(a)(b)(c)(d)In the previous work of using GBT model to uncover the nonlinear relationship between transcriptomic and proteomic datasets, an interesting pattern, a plateau effect, was observed in the partial dependence plot of mRNA abundance versus protein abundance for high values of mRNA. This is potentially due to the fact that protein abundance data is sparse, with high variance, where the tree models do not generate splits among the predictors [19]. In the neural network model, the output of an MLP is an inherent smooth function of the network inputs, hence a relatively smooth prediction is always guaranteed. However, we can still see the similar pattern that the curve becomes flat when the variance for the mRNA abundance is high. For example, in the region of high mRNA values there are a small number of genes/peptides whose protein values ranges from (0, 40) for Dataset 1 and (0, 500) for Dataset 2. This could reflect problems with the sensitivity of the current proteomic technologies for high abundant proteins.


## 3.2. Construction of the Nonlinear Artificial Neural Network Model

A feed-forward network with a single hidden layer was trained for each condition for both datasets with mRNA abundance and other predictors as inputs and protein abundance as target. The neural networks were trained using the functions in the neural network toolbox provided in standard software Matlab (7.0.1). The activation function for the neurons in the hidden layer was the hyperbolic tangent function, which is given by(3)f(v)=1-exp⁡⁡(-v)1+exp⁡⁡(-v).
The activation function of the output neuron is a pure linear function, which is the weighted sum of the outputs from the neurons of the hidden layer. The optimal number of neurons in the hidden layer and the optimal weights were determined using 10-fold cross-validation for training a neural network. For each dataset under each experimental condition, the networks with the same inputs but with different number of neurons (from 5 to 12) in the hidden layer were each trained using cross-validation to learn the relationship between the mRNA features and protein abundance. The optimal number of neurons in the layer was determined based on the one generating the minimum averaged validation error.Both transcriptomic and proteomic data were log-transformed. For dataset 1, the genes with no missing value for three out of four replicates were selected for model training. The numbers of genes used for training were 354, 312, and 308 for the FL, LL, and LS conditions, respectively. For each growth condition, the inputs to the multilayer perceptron were the top ten predictors identified previously [19] (See Table  1 at Supplementary Material). The optimal network had 6, 6, and 9 neurons in the hidden layer for FL, LL, and LS conditions, respectively. Coefficients of determinations R2 were 0.47, 0.52, and 0.51 for FL, LL, LS conditions, respectively. For dataset 2, a total of 986, 1001, and 986 genes with no missing protein measurements were used for model training for CT0, CT120, and ST120 conditions, respectively. For each of the conditions, the inputs to the multilayer perceptron were the top ten predictors identified previously [19] (See Table  1 at Supplementary material). The optimal network had 8, 7, and 6 neurons in the hidden layer for CT0, CT120, and ST120 conditions, respectively. Coefficients of determination R2 were 0.65, 0.62, and 0.68 for CT0, CT120, and ST120 conditions, respectively.


## 3.3. Validation of the Prediction through the Nonlinear Artificial Neural Network Model

The validation of the model prediction was conducted through the calculation of coefficient of variation (CV) across conditions for every operon ofD. vulgaris as described before [18, 19]. The genes in the operons are considered to have less dispersion than any random set of genes in terms of their gene expression and/or protein abundance, because of their intrinsic biological relationship. With relatively low threshold for operon identification as described in Section 2, a total of 609 operons were identified in D. vulgaris, which ranged from 2 to 13 genes per operon. To compare these CV values, we have used a permutation test in the following manner: a CV is computed from the predicted protein values for a permuted number of genes. This step is repeated 2000 times through resampling of genes without replacement. The number of genes randomly selected is equal to the size of the gene group to be compared to. For dataset 1, 82% to 89% of the operon groups had smaller CV values than the permuted groups, and for dataset 2, 79% to 89% of the operon groups had smaller CV values than the permuted, demonstrating in general good performance of the model prediction. The validation results are summarized in Table 1.Table 1
Model validation: correlated expressions of proteins in selective operons.
Dataset 1Dataset 2FLLLLSCT0CT120ST120OperonCVPCVCVPCVCVPCVCVPCVCVPCVCVPCV10.360.7070.60.7970.3010.9011.63*1.2480.8141.1861.2691.28621.6*0.6381.564*0.7251.48*0.6310.7731.1020.9681.0440.7021.12930.4990.6591.04*0.7660.5630.8661.295*1.1341.568*1.1061.324*1.16940.2160.6480.2390.7260.3230.8130.3331.1060.4011.050.4531.13550.4180.6460.5590.7310.450.7940.4811.1050.3981.050.4661.12260.180.6470.4540.7450.3530.8090.4951.1160.4031.0420.4731.12470.2970.6510.1620.7290.3120.80.4081.1010.3431.0420.691.11880.3070.6930.2920.7720.280.9080.7671.2040.5021.1470.571.23690.3890.6760.3520.7640.8510.860.8841.1710.6311.0980.5581.197100.6210.6240.4450.7650.740.8590.8891.1660.671.0961.0871.19110.3040.5950.2310.7040.30.7730.781.1530.8161.0820.5671.182120.5990.8080.8310.9080.7291.0761.2011.5351.1741.3661.563*1.528130.3530.6460.4350.7230.5310.8260.3661.1180.5321.060.3691.124140.3590.6650.4070.730.2670.7520.6531.1570.8251.0870.5971.18150.3810.7170.7070.8010.780.9220.6911.2460.8381.1660.6741.293160.2990.7490.5380.8330.3070.9890.6131.3420.4561.2521.0591.372170.3210.6490.2190.7260.2220.8060.6571.0590.4920.9430.5291.046180.4730.6450.4860.7230.6280.7670.4741.1640.5391.0910.5781.185190.5920.7410.4790.8240.6310.9821.68*1.3171.0941.2481.403*1.342200.767*0.6210.5860.6950.5610.7520.2671.0370.2880.9640.4921.053210.3890.5480.4160.6930.4720.7571.0341.0370.8130.9780.8231.061220.2560.6540.4330.7210.490.80.5321.0950.6021.0530.5691.123230.6310.6940.7070.7750.6320.8770.861.1960.6771.1250.9251.242240.3760.6550.6480.6550.4060.6991.23*1.0990.9271.0311.138*1.123250.869*0.6741*0.7541.194*0.852.337*1.1562.119*1.0942.303*1.185260.657*0.6180.6830.7360.7470.8560.6961.1510.521.1140.71.176270.3780.850.4520.9380.4321.1040.5881.6310.6671.4450.7711.629281.013*0.670.963*0.7641.526*0.9131.2551.2591.375*1.1811.332*1.282CV is computed by dividing SD by the mean of the prediction values for protein abundance for a specific set of genes (group). The protein prediction values were normalized by molecular weight before CV calculation. PCV is the mean of CV values computed through permutation test for selected operons.*CV values greater than PCV.To further demonstrate how the current measurement of dispersion (CV) compares with the distribution of the permuted coefficients of variation, we also calculated the percentile score for all coefficients of variations, based on the biological knowledge that genes from the same operon should be less dispersed than permuted sets of genes and the percentile scores are thus expected to be low. The results showed that even though the calculated CV for most groups was less than the mean CV value for permuted sets of genes, the percentile scores shows that only a small collection of these groups fall within a percentile less than 0.20. The percentage of groups with percentile scores less than 0.2 is 43% to 57%, and 57% to 61% for dataset 1 and dataset 2, respectively.


## 3.4. Comparison of Artificial Neural Network to Previous Methods (GBT and ZIP)

Similar to the GBT method, the artificial neural network method aims to identify an unknown and potentially nonlinear relationship between proteomics and mRNA abundance data, which distinguish them from the ZIP method, which carries zero-inflated linear assumptions [18]. Another key difference is that GBT and neural network methods utilized additional sequence based information, such as cellular role, gene length, amino acid usage, and codon usage known to be important for translation efficiency to model the relationship between transcriptomic and proteomic dataset, while the ZIP method establishes the relationship between proteomic and mRNA solely based on experimental transcriptomic and proteomic dataset.Some improvements were observed in terms of prediction validation by the neural network model in comparison to the GBT model. The percentage of operon groups with CV smaller than permuted CV was 82% to 89% for dataset 1 in comparison to 75% to 79% for the GBT model. The percentage of operon groups with CV smaller than permuted CV was 79% to 89% for dataset 2 by the neural network model in comparison to 75% to 82% by the GBT model. For dataset 1, the percentage of groups with percentile scores less than 0.2 is 43% to 57% in the neural network model in comparison to 36% to 50% in the GBT model. For dataset 2, the percentage of groups with percentile scores less than 0.2 is 57% to 61% in the neural network model in comparison to 32% to 54% in the GBT model. The coefficient of determination was 0.47 to 0.58 for dataset 1 and 0.62 to 0.68 for dataset 2 in the neural network model in comparison to 0.39 to 0.58 for both datasets in the GBT model, demonstrating the artificial neural network model could be a good alternative in analyzing the missing proteomic values.


## 3.5. Analysis of Factors Ascribed to Missing Protein Measurements

Through the neural network model, the majority of the proteins that were experimentally undetected were predicted to have protein abundance with greater than 1.0 peptide hit. This thus tempted for us to look into the possible factors which may be ascribed to the missing protein values. To do so, we collected some of the protein or mRNA/DNA sequence-based features which may affect protein detectability for use along with experimental transcriptomic data in a logistic regression procedure to determine their relative contribution to detectability of proteins in theD. vulagris genome. The features/measurements we used were from the following several categories: (i) features related to gene expressivity: experimental mRNA abundance data, Nc and CAI values [33, 34], codon usage such as CR_axis1-CR_axis4, T3s, C3s, A3s, G3s, GC3s, and GC content [35]; (ii) features related to protein and RNA stabilities: protein instability index and AliphaticIndex, and minimal free energy (MFE) of RNA [35]; (iii) features related to protein solubility: protein solubility indices such as Grand Average of Hydropathy (GRAVY) [33]; (iv) amino acid usage such as AA_axis1-AA_axis4, and others such as gene length [35]. The logistic regression procedure in the SAS software package was used to model the probability of whether a protein will be experimentally detected given the features and measurements for the gene/protein. The stepwise selection method was used to select those most influential factors for indicating missing proteomic value. Those highly correlated features were removed automatically. The 0.05 significance threshold was used for the stepwise selection. The mRNA abundance and gene length values were both log2 transformed before being used in the procedure. The protein stability index, aliphaticIndex, and Nc values were divided by 10 and the CAI, N_MFE, GRAVY, AA_axis1-AA_axis4, and CR_axis1-CR_axis4 were all multiplied by 10 so that they scaled to 1. This standardization will not change the results, but make illustration of the results easier.Tables2 and 3 present the odds ratio of each selected factor for dataset 1 and dataset 2, respectively. The results showed that some of the features are highly correlated with probability of missing proteomic values. For dataset 1, seven factors were indentified consistently by the regression model. Among the features that are positively correlated with probability of missing proteomic values (i.e., when the features/measurement increase), protein instability index which provides an estimate of the stability of a given protein in a test tube [38] was found as the most significant factor for all three growth conditions, and aliphaticIndex which is regarded as a positive factor for the increase of thermostability of globular proteins [39], was found significant for growth conditions LL and FL, but not for condition LS. In addition, AA_axis2 and CR_axis3 correlated most strongly with GRAVY and G3s, respectively, and they were also found significant in all three growth conditions, while the N_MFE was only significant for LS condition (Table 2). For dataset 2, the features that are kept in the regression model are more the same for the three growth conditions in dataset 1, suggesting that these factors could be universal across experimental platforms. Six features that are significant for dataset 1 are also found to be critical for all the conditions and affect the probability in the same direction for dataset 2 include mRNA abundance, protein instability index, gene length, Nc, CAI, and AA_axis2 corresponding to GRAVY (Table 3). In general, the analysis suggested that those proteins with smaller mRNA expression values, smaller Nc value, higher protein instability index, or high value of aliphaticIndex, have a greater probability to have missing proteomic values in the D. vulgaris datasets. Based on this preliminary analysis, it may be possible to apply a logistic regression model to predict whether a given protein will or will not be detected by choosing a proper threshold for the predicted probability.Table 2
Estimates of odds ratios of missing proteomic values for dataset 1.
Dataset 1LLFLLSeffectPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitslog2 of mRNA0.4580.4050.5170.4420.3920.4990.4050.3480.471Protein instability index/101.6281.3791.9221.4931.2821.7391.5391.2971.825log2 of gene length0.4770.3730.6080.5170.4150.6440.4190.3240.543Nc/100.4650.2710.7980.4460.2720.7320.3140.1740.566CAI*100.2950.2120.4100.2750.2040.3700.2220.1550.319AA_axis2*101.3681.1931.5691.3261.1701.5021.3231.1661.501CR_axis3*101.2641.0461.5281.3091.1121.5411.3271.0921.613N_MFE*10——————1.8801.1783.000AliphaticIndex/101.1401.0131.2821.1601.0401.294———Table 3
Estimates of odds ratio of missing proteomic values for dataset 2.
Dataset 2CT0CT120ST120effectPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitsPoint estimate95% wald confidence limitslog2 of mRNA0.4070.3700.4480.4510.4120.4940.4990.4590.542Protein instability index/101.3101.1611.4781.2831.1401.4441.2641.1241.421log2 of gene length0.2610.2160.3150.2870.2400.3450.2940.2470.352Nc/100.5120.3360.7800.5180.3430.7830.5520.3740.815GRAVY*100.9140.8620.9700.9140.8630.9680.9230.8730.977CAI*100.5990.4210.8510.5590.3980.7860.5360.3840.749AA_axis1*100.6860.6270.7500.6810.6240.7440.6860.6290.749AA_axis2*101.8111.5582.1061.7211.4871.9901.6741.4521.931AA_axis4*100.8270.7470.9160.8370.7570.9250.8460.7670.934CR_axis1*100.7710.6470.9190.7800.6570.9250.7810.6600.923CR_axis4*100.7160.6350.8080.7060.6270.7950.7250.6450.814


## 3.6. Data Interpretation Using Imputed Datasets

A total of 308–354 proteins were experimentally identified for dataset 1, and 986–1001 proteins were experimentally identified for dataset 2, which are slightly over 10–30% of the whole proteome inD. vulgaris. Through the ANN modeling and imputation by integrating transcriptomic and proteomic data, we were able to assign abundance values to approximately 33–3400 proteins (dependent on growth conditions), approximately 96–98% of the whole proteomes inD. vulgaris. Using the imputed datasets, we first analyzed the proteins involved in energy metabolism when D. vulgaris was grown in three cultivation conditions (LL: Lactate-based growth in Exponential phase; LS: Lactate-based growth at Stationary phase; FL: Formate-based growth at Exponential phase). As shown in Table 4, even for proteins involved in energy metabolism such as ATP synthase which are typically expressed at relatively high levels, significant missing proteomic data was still obvious in the experimental datasets, only two out of nine putative ATP synthase proteins were experimentally detected under three growth conditions. By analysis of the imputed protein abundance for all putative ATP synthase proteins, a clear trend of downregulation was observed for LS and FL condition when compared with LL condition, which is consistent with the previous data suggesting lactate is a favorable electron donor for D. vulgaris [22]. When comparing lactate-based growth with formate-based at the exponential phase, we found that enzymes involved for formate utilization were predicted to be upregulated during growth on formate (Table 4). It has been reported that c-type cytochromes are highly expressed in sulphate-reducing bacteria [20], however, only low amounts of c-type cytochrome proteins have been identified from D. vulgaris by the proteomic (Table 4), probably due to the fact that c-type cytochromes undergo a complex posttranslational maturation process involving covalent attachment of heme groups [40]. Through the ANN model, significant expression of various cytochrome c proteins was predicted in D. vulgaris grown under three conditions. In general, cytochrome c proteins were predicted to have higher expression at the exponential phase in both lactate- and formate-based media than the stationary phase, consistent with their important roles in electron transfer during fast growth [20]. Two cytoplasmic hydrogenases (Ech and Coo) were previously assigned the putative function of generating hydrogen in the cytoplasm [20]. In agreement with previous prediction from the ZIP model [18], the ANN model predicted that all subunits of Ech hydrogenase (EchACDE) were upregulated, while all subunits of Coo hydrogenase (CooHLUX) were down-regulated on formate-based growth, suggesting their differential roles in H2 generation when grown with different carbon sources. We also used the imputed data to analysis the stress responses when the D. vulgaris cells were treated under salt stress. The analysis allowed identification of some proteins which were not experimentally detected, however, were predicted to be upregulated under the salt treatment condition (ST120) as compared with control (CT120) (Table 5). Among them were transporters and various multidrug resistance proteins and some signal transduction proteins. This list of proteins can serve as putative candidates for further experimental validation. In general, our analysis indicates that biological interpretation may benefit from imputing missing data using computational methods by integrating temporal transcriptomic and proteomic data. The full list of all protein-abundance predictions for all proteins has been provided in the Table  2 of the supplementary material.Table 4
Regulation of proteins involved in energy metabolism for Dataset 1*.
Gene_IDProtein abundance-Predicted-LL conditionRNA abundance-Measured-LL conditionProtein abundance-Measured-LL conditionProtein abundance-Predicted-LS conditionRNA abundance-Measured-LS conditionProtein abundance-Measured-LS conditionProtein abundance-Predicted-FL conditionRNA abundance-Measured-FL conditionProtein abundance-Measured-FL conditionDescriptionDVU091810.082315.20x2.46416.75x5.711668.90xATP synthase F0, A subunit (atpB)DVU07796.552505.40x3.30497.75x5.081397.10xATP synthase F0, B subunitDVU07806.623156.20x3.881050.10x5.812352.70xATP synthase F0, B subunitDVU09177.606086.50x3.851234.90x10.855944.40xATP synthase F0, C subunit (atpE)DVU09203.801590.50x3.59947.88x5.751198.20xATP synthase protein I (atpI)DVU077739.242535.0020.008.40369.958.0017.76735.6817.33ATP synthase, F1 alpha subunit (atpA)DVU077550.112228.6049.3323.76711.5330.3339.801198.4038.33ATP synthase, F1 beta subunit (atpD)DVU07784.251974.50x2.71374.22x3.21714.57xATP synthase, F1 delta subunit (atpH)DVU07744.852043.70x2.32336.55x4.141113.40xATP synthase, F1 epsilon subunit (atpC)DVU07022.62288.72x2.20372.82x2.05417.681.00Cytochrome c family proteinDVU09221.6966.60x1.85149.38x1.49149.52xCytochrome c family proteinDVU12883.441172.20x2.06298.10x3.12826.90xCytochrome c family proteinDVU248310.24456.23x3.24284.22x7.36490.33xCytochrome c family proteinDVU24845.06575.50x2.39277.35x4.63665.95xCytochrome c family proteinDVU27914.77683.65x2.61530.00x3.32667.05xCytochrome c family proteinDVU310715.101009.302.332.98177.151.6713.041026.303.00Cytochrome c family proteinDVU31445.58297.85x3.58142.45x8.24342.52xCytochrome c family proteinDVU062530.201822.80x14.80764.85x16.612467.50xCytochrome c nitrite reductase, catalytic subunit NfrADVU18154.91439.87x2.98324.93x7.12833.88xCytochrome c oxidase, subunit IDVU18124.42304.57x2.72215.70x4.63335.70xCytochrome c oxidase, subunit IIDVU18142.53397.35x2.31283.90x3.80898.67xCytochrome c oxidase, subunit IIIDVU28093.11447.42x2.32395.15x2.10601.50xCytochrome c3DVU317113.243937.003.335.24966.25x11.825719.201.00Cytochrome c3DVU25241.77114.95x2.08135.07x1.61170.70xCytochrome c3,DVU30413.03894.60x2.62713.10x2.73684.27xCytochrome c553DVU18176.342280.0011.677.261858.203.339.615400.2012.00Cytochrome c-553 (cyf)DVU10482.31255.07x2.14243.58x1.53185.47xCytochrome c-type biogenesis protein CcmB (ccmB)DVU10472.87435.45x2.06284.87x2.23259.25xCytochrome c-type biogenesis protein CcmC (ccmC)DVU10513.30915.92x3.59769.45x4.19732.07xCytochrome c-type biogenesis protein CcmE (ccmE)DVU10507.42581.02x4.03408.20x6.27463.25xCytochrome c-type biogenesis protein CcmF (ccmF)DVU327135.798217.80x12.152295.10x11.386954.00xCytochrome d ubiquinol oxidase, subunit I (cydA)DVU327018.093877.70x5.501145.70x13.903897.10xCytochrome d ubiquinol oxidase, subunit II (cydB)DVU04341.28177.22x2.79168.72x5.13468.00xEch hydrogenase, subunit EchADVU04332.10217.82x2.34166.95x2.84422.28xEch hydrogenase, subunit EchBDVU04322.11222.90x2.24133.13x3.12595.60xEch hydrogenase, subunit EchCDVU04312.20303.27x2.03215.73x3.30778.70xEch hydrogenase, subunit EchDDVU04302.68259.00x2.52199.20x4.61685.78xEch hydrogenase, subunit EchEDVU04292.24239.03x2.07201.57x2.05405.30xEch hydrogenase, subunit EchFDVU28240.90124.20x3.54106.78x9.30102.10xFormate acetyltransferaseDVU22721.0597.08x3.69113.331.0016.21160.88xFormate acetyltransferaseDVU05782.11154.12x2.23217.48x2.03247.52xFormate dehydrogenase accessory protein FdhDDVU05772.56237.93x2.29148.02x1.97513.40xFormate dehydrogenase formation protein FdhEDVU28103.39423.00x2.27312.02x2.47510.12xFormate dehydrogenase formation protein FdhEDVU05880.9235.70x2.6535.25x1.2679.03xFormate dehydrogenase, beta subunitDVU24812.39227.43x2.37304.92x2.35399.82xFormate dehydrogenase, beta subunitDVU28113.22335.33x2.52181.45x3.13435.95xFormate dehydrogenase, beta subunitDVU22882.68781.45x2.24412.85x1.39139.17xHydrogenase, CooL subunitDVU228616.67308.33x3.29138.72x6.0881.55xHydrogenase, CooM subunitDVU22902.84703.223.002.10316.453.00xHydrogenase, CooU subunitDVU22893.92823.18x2.14553.33x2.36319.80xHydrogenase, CooX subunit*Protein abundance is in peptide counts; RNA abundance is in raw fluorescence measurements in microarrays. “x’’ indicated proteins are undetected experimentally.Table 5
Up-regulated proteins under stress condition for Dataset 2*.
Gene_IDProtein abundance-Predicted-CT120RNA abundance-Measured-CT120Protein abundance-Measured-CT120Protein abundance-Predicted-ST120RNA abundance-Measured-ST120Protein abundance-Measured-ST120DescriptionAmino acid biosynthesisDVU02855.003005.10x7.434051.50xImidazole glycerol phosphate synthase (hisH)DVU166521.8220264.00x50.5530766.00x3-dehydroquinate dehydratase, type II (aroQ)Biosynthesis of cofactors, prosthetic groups, and carriersDVU01575.457559.00x6.315092.30xThiamin-monophosphate kinase (thiL)DVU09315.411730.40x6.874401.40xPhosphomethylpyrimidine kinase (thiD)DVU12005.825439.20x6.755079.60xRiboflavin synthase, alpha subunit (ribE)DVU27492.851427.30x3.59856.35xPrecorrin-6Y C5,15-methyltransferase (decarboxylating) (cobL)DVU330714.6520224.00x26.4230616.00x3-octaprenyl-4-hydroxybenzoate carboxy-lyase (ubiX)Cell envelopeDVU02942.441034.50x3.15726.82xGlycosyl transferase, group 2 family proteinDVU03082.971552.00x7.936986.70xMembrane proteinDVU03288.595628.30x9.125042.40xGlycosyl transferase, group 1 family proteinDVU040726.6520237.00x38.7032568.00xRare lipoprotein A family proteinDVU15342.481792.60x3.562183.30xMembrane proteinDVU20532.982406.60x3.852507.60xMembrane proteinDVU21022.91628.30x3.88843.26xOuter membrane protein, OMP85 familyDVU235620.546021.30x36.267883.90xMembrane proteinDVU26723.303253.30x4.475875.90xMembrane proteinDVU27255.204925.80x13.978894.70xMembrane proteinCellular processesDVU030228.3522358.00x37.2121180.00xChemotaxis protein CheXDVU20692.18811.21x3.72658.60xDNA processing protein DprADVU28164.82539.18x8.651502.80xMultidrug resistance proteinDVU28934.599413.30x6.789241.60xFlagellar basal-body rod proteinDVU33263.13782.67x3.682732.60xMultidrug resistance protein, Smr familyDVUA00025.532257.50x12.966384.00xParA family proteinDNA metabolismDVU08783.201967.30x7.1311149.00xdnaK suppressor proteinDVU13532.261516.00x10.31868.27xDNA polymerase III, alpha subunit (dnaE)DVU17031.93906.31x2.69628.53xType I restriction-modification enzyme, R subunitEnergy metabolismDVU00804.422476.60x6.455175.70xFumarate hydratase, class II (fumC)DVU05316.57328.82x11.021157.90xHmc operon protein 6DVU21373.64389.78x4.45338.82xSuccinyl-CoA synthase, beta/alpha subunits (sucCD)Other categoriesDVU01973.46933.13x4.40285.54xPhage portal protein, lambda familyDVU14904.72351.64x5.70372.18xTail tape measure proteinRegulatory functionsDVU00633.743823.20x4.796404.70xTranscriptional regulator, MarR familyDVU13317.546812.50x8.3311619.00xTranscriptional regulator, LysR familyDVU16457.227075.60x14.0323482.00xTranscriptional regulator, ArsR familyDVU30622.341205.40x2.92570.26xSensor histidine kinase/response regulatorDVU309512.3017923.00x14.8623528.00xTranscriptional regulator, Fur familyDVU07214.201085.90x7.19513.28xSensory box histidine kinaseDVU08345.732795.50x12.1713860.00xRibonuclease HII (rnhB)DVU107513.6330772.00x15.4527576.00xRibonuclease P protein component (rnpA)Transport and binding proteinsDVU05932.58992.23x3.281151.30xL-lysine exporterDVU10683.92301.98x5.76270.52xBranched-chain amino acid ABC transporter, permease proteinsDVU23406.0011276.00x7.0011625.00xAmino acid ABC transporter, permease proteinDVU25723.192831.50x5.777937.50xFerrous iron transport protein ADVU25744.625005.50x8.858857.00xFerrous ion transport proteinDVU274417.961201.00x34.423053.50xHigh-affinity branched-chain amino acid ABC transporterDVUA01363.29673.57x3.89346.83xZinc transporter (zupT)Unknown functionDVU01558.3516519.00x12.1518298.00xType I phosphodiesterase/nucleotide pyrophosphatase family proteinsDVU093810.1215409.00x11.7613138.00xIsoamylase N-terminal domain proteinDVU10237.26123.66x8.15144.26xRhomboid family proteinDVU12025.393480.30x5.925395.40xCytidine/deoxycytidylate deaminase family proteinDVU15585.77189.07x6.47234.72xCysteine-rich domain/iron-sulfur cluster-binding domain proteinsDVU17303.161288.70x3.744503.10xDNA-binding proteinDVU17456.55319.33x7.58182.38xDNA-binding proteinDVU17475.08258.33x12.80230.72xATPase, histidine kinase-, DNA gyrase B-, and HSP90-like domainsproteinDVU20364.271086.70x5.07350.61xHelix-turn-helix protein, CopG familyDVUA00393.75269.47x4.27261.42xChain length determinant family protein*Protein abundance is in peptide counts; RNA abundance is in raw fluorescence measurements in microarrays. “x’’ indicated proteins are undetected experimentally.In conclusion, missing values in large-scale proteomic analysis is a frequent problem and in some cases it has created difficulty for accurate interpretation of proteomic data for complex biological systems. However, in spite of its obvious importance, the methodologies to deal with missing proteomic data are still lagging well behind the methods for transcriptomic analysis. In this study, we applied an artificial neural network approach to approximate the relationships between cognate transcriptomic and proteomic datasets ofD. vulgaris, and to predict protein abundance for the proteins not experimentally detected based on relevant predictors, such as mRNA abundance, cellular role, and triple codon counts. The results showed that the coefficients of determination for the trained neural network models ranged from 0.47 to 0.68, providing better description than several previous models [18, 19]. In an attempt to identify possible causes for missing proteomic values, a logistic regression analysis of various experimental measurements and sequence-based features was performed for the experimentally undetected proteins but with high predicted protein abundance values. Although the results are in general consistent with previous knowledge that proteins with low possible expressivity and high instability will have higher chance to be missed, in the study, we for the first time quantified the odds ratio of each of those possible factors in the D. vulgaris datasets. In addition, in a preliminary analysis we demonstrated that biological interpretation could be improved by using the imputed proteomic datasets. Finally, although the model was well validated statistically, caution needs to be exercised when interpreting experimental data based on predicted protein expression values because the predicted abundance values are constrained by the quality of experimental proteomic data used as input, and a dearth of large-scale quantitative predictors which can be included in the model. Nevertheless, the initial success in applying the ANN method is encouraging, and the model could serve as a basis for developing more sophisticated models to further improve the biology interpretation.



---

*Source: 780973-2011-05-04.xml*