BACKGROUND
bioavailability processes were defined in a  <dig> national research council  report  <cit>  as the individual physical, chemical, and biological interactions that determine the exposure of plants and animals to chemicals associated with soils and sediments.
in environmental risk assessment, the amount of chemicals taken up by an animal or plant is termed dose, or, interchangeably tissue residue, body burden or chemical bioavailability.
a distinction exists between dose and exposure as the latter is defined as the amount of chemicals present in the immediate environment where the organism is exposed to.
soils are a major sink for many environmental contaminants including explosives compounds.
soil contamination by military unique compounds is a serious environmental concern that can result in the formation of chemical residues in tissue of exposed organisms  <cit> .
since many site-specific biotic and abiotic factors can modify the form, mobility and availability of these contaminants, the actual exposure risk to ecological receptors may be less than that suggested by their total concentration.
therefore, the extent to which chemicals are bio-available has significant implications for risk management and remedial decision-making at contaminated sites.
a myriad of biological, physical and chemical approaches have been used to evaluate bioavailability of chemicals  <cit> .
for instance, biological approaches include techniques measuring organismal uptake , response, and toxicity.
tissue residues or body burden of chemicals, often determined using analytical chemistry methods, can represent the bio-available or biologically effective concentration at the target site.
however, little consensus exists about optimal approaches for measuring bioavailability.
as summarized in the aforementioned nrc report  <cit>  and a recent review  <cit> , no single tool is free from limitations and none can be applied universally.
an intensive effort to develop mechanistic tools or models based on mechanisms is critical to future development of bioavailability tools  <cit> .
recent advances have prompted the application of genomics-based technologies to ecological risk assessment, including screening, tiered testing, monitoring, remediation, and regulatory decision-making .
toxicogenomics has been extensively applied to assess toxicological effects, especially in biomarker discovery and toxicity mechanistic investigations.
on the other hand, toxicogenomics also possesses a great potential for providing a quantitative measure of chemical exposure.
nevertheless, the application of toxicogenomics to exposure assessment has been relatively under-explored.
previously, we applied microarray technology to profile gene expression in earthworms exposed to explosives  <dig> ,6-trinitrotoluene  or  <dig> ,5-trinitro- <dig> ,5-triazacyclohexane  for  <dig> or 14 days, and built multivariate regression models to quantitatively predict earthworm tissue residue of these two compounds  <cit> .
the models, however, showed only a modest predictive power, explaining close to half of the variance for tnt tissue residue and one-quarter of the variance for rdx tissue residue.
the present study was motivated to improve the quantitative predictive power of regression modeling based on genome-wide gene expression data.
our hypothesis was that small sets of predictor genes could be identified and used to build multivariate regression models for quantitative prediction of tissue residue levels of explosives compounds.
consistent with our previous study, the overall goal of this study was to investigate the feasibility of using gene expression data to assess animal exposure.
to achieve this goal, we re-analyzed the previous datasets by separating the 4-day exposure from the 14-day exposure, and expanding from two regression methods to  <dig> methods.
in addition, we have also generated a new dataset from hmx-exposed earthworms  and applied the same approach to analyzing it.
our results demonstrated that much higher prediction accuracies were attained indicating that microarray gene expression coupled with multivariate regression modeling is a viable approach for assessing chemical bioavailability.
methods
the overall experimental approach is depicted in fig.
<dig>  we first exposed earthworms to explosives compounds and measured tissue residues  and gene expression profiles in both exposed earthworms and unexposed controls.
then we identified small sets of predictor genes  to build regression models , where y’ is the predicted residue and β is a coefficient vector.
these models were applied to predict the bioavailability  in exposed worms using only predictor gene expression data.
finally, the actual measured tissue residues  were compared with the predicted values  to determine prediction accuracy evaluated by correlation coefficients .
our ultimate goal was to apply these prediction models to estimate the unknown bioavailability in worms from their gene expression data .fig.
<dig> the overall experimental approach.
see the methods section for explanation
experimental design
previously, we collected a large 248-sample dataset with measurements of transcriptome-wide gene expression  <cit>  and tissue residue of two explosives chemicals, -labeled tnt and rdx, in the earthworm eisenia fetida .
briefly, three sets of experiments were conducted by exposing mature adult worms bearing a clear clitellum in a pristine sandy loam soil amended with tnt or rdx.
nominal exposure concentrations were as follows:  <dig>   <dig>   <dig>   <dig>   <dig>  96 mg  tnt/kg, or  <dig>   <dig>   <dig>   <dig>  128 mg  rdx/kg for 4 days  or 14 days .
the 4-day exposure was repeated  with the same tnt concentrations but different rdx concentrations .
following the same experimental design, we exposed adult worms to -labeled hmx for  <dig>   <dig> and 28 days at nominal concentrations of  <dig> ,  <dig>   <dig>   <dig>   <dig> and 128 mg/kg in the same silt loam soil as in our previous study  <cit> .
hmx was purchased from dupont nen  with an initial specific activity of  <dig>  mci/mmol.
earthworms were reared in-house in a continuous lab culture as previously described  <cit> .
ten worms were exposed in 250 g  of soil per each treatment.
the nominal soil concentrations were verified to possess less than 10 % variations from target concentrations using a high performance liquid chromatography equipped with a radioactivity flow detector   <cit> .
upon termination, all worms were flash-frozen in liquid nitrogen, stored at−80 °c, lyophilized at −40 °c, and homogenized.
tissue residue analysis for hmx
all ten worms per treatment were analyzed for tissue residue.
for tissue hmx-residue analysis, triplicate subsamples  were digested in 1 ml of  <dig>  n scintigest tissue solubilizer  for 18 h  <cit> .
radioactivity was determined by adding 200 μl of tissue digest to ultima gold scintillation cocktail  using a packard tricarb 2500tr liquid scintillation counter .
prior to the analysis, we validated this method  by comparing it with our previously used methods.
testing results showed that the solubilization method yielded 92 % ± 3 %  of the radioactivity measured previously by using oxidization methods for the same tnt- or rdx-exposed worms  <cit> .
gene expression profiling for hmx-exposed worms
total rna was extracted from five of the ten worms per treatment using an rneasy mini kit .
each rna sample was hybridized to the custom-designed agilent 15 k e. fetida oligo array  previously used for generating the tnt and rdx microarray datasets  <cit> .
this array contained  <dig> unique transcript-targeted 60-mer oligo probes .
details of array hybridization, gene expression data acquisition and pre-processing were described elsewhere  <cit> .
identification of differentially expressed  genes
de genes were identified for all three datasets  among multiple concentrations using a multivariate permutation random-variance t-test  or f-test  implemented in brb-arraytools version  <dig> .
<dig>  <cit> .
a gene was considered statistically significant if it achieved 80 % confidence that the false discovery rate  was less than 10 %.
de genes were derived for individual exposure duration and explosives compound separately.
regression prediction modeling
for residue prediction, we chose  <dig> different multivariate regression models and employed double-looped, 10-fold cross-validation as described in statnikov et al.
<cit>  to assess prediction accuracy.
the inner loop was used to determine the optimal value of parameters  for training in the outer loop.
model performance  was estimated in the outer loop by training on all splits but one, and using the remaining one for testing.
coefficient of determination  was calculated using pearson’s formula to describe prediction accuracy or “goodness of fit”, i.e., how well the tissue residue predicted by a regression model represented the actually measured residue of a worm tissue sample.
the  <dig> models include six linear models ) and  <dig> nonlinear models  .
matlab codes were scripted to implement array data preprocessing, regression, and cross-validation, and are available upon request.
RESULTS
selection of predictor genes from tnt and rdx gene expression datasets
we reanalyzed the tnt and rdx microarray datasets by separating the 4-day sample set from the 14-day sample set.
differentially expressed  genes were identified among different treatments  of 4-day or 14-day earthworm samples using a multivariate permutation test  <cit> .
statistical reanalysis of the 14-day exposure gene expression data resulted in six and three de genes for tnt and rdx, respectively .
for the 4-day tnt exposure,  <dig>   <dig> and  <dig> genes were inferred as de genes from the original exposure , the repeat exposure  and the original vs. repeat controls , respectively .
a group of  <dig> de genes were found to be common between the original and the repeat tnt exposures.
this group was further reduced to  <dig> genes after excluding genes also found to be significantly altered between the original and the repeat controls .
similarly,  <dig> and  <dig> de genes were derived from the original rdx and the repeat rdx exposures , respectively, with  <dig> genes in common .
twenty-six genes out of the 178-gene group remained after excluding the same de genes appearing in the controls comparison .
the final sets of identified de genes are shown in additional file  <dig> and additional file  <dig>  the low numbers of de genes found in 14-day exposures are consistent with our previous report  <cit> .
we obtained more de genes than previously for the 4-day exposures because of the reduced statistics stringency .
however, worms used in the original and the repeat exposures exhibited significant differences that were reflected as nearly 1/ <dig> of all 15 k profiled genes differentially expressed between the two control groups .
therefore, we chose to remove these genes not responding specifically to 4-day tnt or rdx exposure from the final de gene lists .
regression predictive modeling for tnt and rdx tissue residue
using the aforesaid final sets of de genes as predictor genes, the  <dig> regression methods displayed varied power in predicting worm tissue residues of tnt and rdx .
the predictive power was assessed using the coefficient of determination  as a measure of the accuracy of the data model.
no single method was placed as the best performer for all four datasets.
for instance, the multivariate linear regression model was the best performer for the 14-day tnt exposure dataset, but its performance was relatively weak on the 4-day tnt exposure dataset.
quite a few models such as lasso, elastic net, the ridge family models, the svr family models except for svr sigmoid, and nadaraya-watson performed consistently well across all four datasets.
in contrast, some models like svr sigmoid, loglog, and reverse regression performed poorly or even appeared inapplicable to the datasets.table  <dig> performance of  <dig> regression modeling methods on four datasets assessed by coefficient of determination  estimated from ten runs of 10-fold cross-validation with values of the best performing method for each dataset shown in bold
<dig>  ±  <dig>
<dig>  ±  <dig>
<dig>  ±  <dig>
<dig>  ±  <dig>
rdx_d <dig> 4-day rdx exposure, rdx_d <dig> 14-day rdx exposure, tnt_d <dig> 4-day tnt exposure, tnt_d <dig> 14-day tnt exposure, na not available.
see additional file  <dig> for the lists and annotation of predictor genes
the predictive power of the best performers was remarkably improved in comparison to our previously published results  <cit> .
on average, these models explained 75 %  or 66 %  variance of the 4-day samples, and 72 %  or 71 %  variance of the 14-day samples .
in our previous study, the best prediction models explained roughly one-quarter and less than one-half of the variance for rdx and tnt, respectively  <cit> .fig.
<dig> the average predicted versus the measured tissue residues of tnt or rdx in all 4-day or 14-day exposed samples using their respective best performing models
tissue residue in hmx-exposed earthworms
the measured tissue residue of radio-labeled hmx increased with the increasing nominal amendment concentration and also with the duration of exposure .
trace amounts of hmx were detected in some of the two control groups, which are equivalent to the background noise level or the lower detection limit of the analytical method.
the worm tissue residue did not appear to plateau at the highest amendment level, suggesting that more hmx could be taken up from soil by the earthworm, given a higher amendment concentration.fig.
<dig> tissue residue of radio-labeled hmx measured in earthworms exposed for 4-, 14-, and 28-days .
data are represented as mean  + standard deviation  with n =  <dig>  bc = blank control; sc = solvent control
predictor genes selection for hmx gene expression dataset
the hmx microarray gene expression dataset consists of three exposure groups  of  <dig> arrays, i.e.,  <dig> treatments  ×  <dig> replicates.
the two control groups in all three exposures showed no significant difference because no de gene was identifiable at the same settings as above.
so, the two groups were combined as one control group for subsequent statistical analyses.
the three day  <dig> samplings  exhibited little difference with only seven derived de genes, suggesting that worms used in all three exposures were nearly identical and that animal batches had a minimal impact on the test results .
only  <dig>   <dig>  and six genes were inferred as de genes for the 4-, 14-, and 28-day exposures, respectively .
there was no overlap among all four sets of inferred de genes .
using de genes as predictors, regression modeling resulted in a modest predictive power for the 4-day  and 14-day  exposures with the best performer only explaining an average of 55 % variance .
performance of the  <dig> predictor genes on the 28-day dataset was much better with the highest r <dig> of  <dig> .
all these preliminary results are provided in additional file  <dig>
in order to improve the predictive power, a different approach was used to identify additional predictor genes.
we ran a correlation test to determine the degree of dependence between each expressed gene and the tissue residue.
then, all 15 k genes were ranked according to their coefficients of correlation .
with a cutoff of |r| =  <dig> , we obtained  <dig>   <dig> and  <dig> most correlated genes for the 4-d, 14-d and 28-d exposures, respectively, which were considered potential predictor genes .
genes at the top of the most correlated genes lists were most positively correlated with tissue residue, while those at the bottom most negatively correlated with tissue residue.
some of the genes on these lists were also de genes.
we conducted performance analyses of the incremental new sets of predictor genes by adding two or three genes from the top or the bottom of the most correlated genes list at a time to the de genes to form a new predictor gene set.
results  indicate that addition of the top two and the bottom two had an optimal enhancement in prediction accuracy.
therefore, the final revised sets of predictor genes consisted of de genes, the top two and the bottom two genes on the most correlated genes lists.
prediction outcome for hmx gene expression dataset
using the revised sets of predictor genes , the regression models achieved substantial enhancement in prediction accuracy .
particularly, prediction outcomes for the 4- and 14-day exposures improved by nearly 20 % as the best performers for the two datasets both explained 72–73 % of the variance .
the 28-day dataset did not improve as much  largely because five out of the six de genes  were among the most correlated genes list with one  being the most negatively correlated gene  and another  the third most positively correlated gene  .
in contrast, none of the four de genes derived from the 4- and 14-day datasets were on the most correlated genes lists.
models that performed well on the tnt and rdx datasets also did consistently well on the hmx datasets.
three models, svr sigmoid, inverse regression and loglog regression, were not suitable for all datasets , probably because of transformation and normalization operations in data pre-processing.table  <dig> performance of  <dig> regression modeling methods on the three hmx exposure datasets assessed by coefficient of determination  estimated from ten runs of 10-fold cross-validation with values of the best performing method shown in bold
<dig>  ±  <dig>
<dig>  ±  <dig>
<dig>  ±  <dig>
d <dig> 4-day hmx exposure, d <dig> 14-day hmx exposure, d <dig> 28-day hmx exposure, na not available.
see additional file  <dig> for the lists and annotation of predictor genes
fig.
<dig> prediction results of 4-, 14- and 28-day hmx-exposed earthworm tissue residues using the best performing models
discussion
quantitative prediction of a phenotype or trait using high dimensional gene expression data has been long explored in many research fields such as human diseases  <cit> , animal breeding  <cit>  and plant genetics  <cit> .
the phenotypes of prediction interest are broad and have included thrombocytosis etiologies  <cit> , acute myeloid leukemia resistance  <cit>  and breast cancer tumor response  <cit>  to chemotherapy in biomedical research, cattle milk yield in animal breeding  <cit> , and wheat grain yield  <cit>  and plant pathogen infection severity  <cit>  in plant genetics.
the methodology of predictive modeling varies widely from bayesian network-based approaches  <cit>  to bayesian hierarchical regression modeling  <cit>  and ordinary multivariate regression  <cit> .
the reported prediction accuracy also varies from one study to another, implying significant challenges and opportunities co-existing in this field  <cit> .
similar to the aforementioned fields of research, we have also sought to use transcriptomic data coupled with regression modeling to predict a phenotype, i.e., chemical residues in animals  <cit> .
although prediction accuracies were unsatisfactory in our previous study, we have significantly improved them in the present study without adding to the complexity of regression models.
the size of predictor gene sets has been reduced to three to ten genes for all datasets except the 4-day tnt and rdx exposures , which are likely the lowest number of predictors possible  <cit> , given the complicated processes involved in residue formation.
although there is still room for further improvement, the accuracies obtained here are comparable to those reported in other quantitative trait prediction studies .
it has been reported that variable or feature selection is a key step towards accurate prediction  <cit> .
to select features  that have high predictive power as predictors, many different strategies have been developed, including t- or f-test based approaches to identify differentially expressed genes , correlation based methods to rank genes , machine learning based techniques , and some combinations of the three types of strategies .
in the present study, we have used both f-test and correlation based approaches to select for subsets of predictor genes.
a variety of regression models were applied to describe relationships between the response variable  and a set of predictors through a regression function involving some parameter to be estimated from data  <cit> .
given the time-dependent nature of tissue residue  and gene expression levels   <cit> , the datasets were split by exposure duration.
results indicate that this step remarkably improved the prediction outcome for both tnt and rdx exposures.
to be used as predictors, altered levels of a given transcript need not be connected to a specific biologic end point, nor do the specific functions of all the mrna molecules have to be known, although such information would be valuable  <cit> .
bioinformatic annotation of all  <dig> selected earthworm predictor genes using blast2go  <cit>  shows that only  <dig> of them  have meaningful biological functions and that their contribution to residue formation remain largely unknown .
recently, there has been an intense interest in performing variable selection simultaneously with parameter estimation in predictive modeling.
frequently applied approaches include bayesian methods , stochastic search variable selection , adaptive shrinkage with jeffreys’ prior or a laplacian prior, and reversible jump markov chain monte carlo )  <cit>  and lasso, ridge or elastic net regularization methods  <cit> .
given such a large variety, there is no consensus with regard to what method is the universally best performer.
for instance, o’hara and sillanpaa  <cit>  tested several bayesian variable selection methods on both simulated and experimentally collected data and concluded that ssvs, reversible jump mcmc and adaptive shrinkage methods all worked well, but the choice of which method was better depended on the priors that were used, and also on how they were implemented.
similarly, fu et al.
<cit>  concluded that svr, partial least squares regression and multiple linear regression yielded higher prediction accuracies for one dataset but transcriptome-based distances worked better on another dataset.
in the present study, we also observed that lasso, ridge and elastic net regularization methods performed almost equally well on all datasets, whereas inverse and loglog regression methods performed poorly on the datasets in a consistent fashion.
it is worth noting that the approach employed in this study faces several challenges when applied to residue prediction.
first, compared with other existing approaches, it requires a more extensive effort to identify and optimize a set of predictor genes.
second, if the species of interest does not have a transcriptome-wide microarray available, one has to either design the microarray from scratch or use the array of a closely related species , which may limit its applicability.
third, as a novel approach, there is no doubt that it is still in its embryonic phase and that its full potential and limitations are not yet explored thoroughly.
for instance, our approach may be applied to quantitative predictions in drug discovery such as predictive ranking of new drug toxicity and/or potency as they share similarities in high data dimensionality and transcriptomic profiling based on either microarray- or rna-seq/next-generation sequencing-based data.
CONCLUSIONS
chemical residue in exposed animals is an important dose metric in environmental risk assessment.
the formation of tissue residue involves a number of complex biological processes which can be reflected as expression profiles in microarray experiments.
building predictive models of tissue residue based on gene expression would help to accurately assess how much chemical an animal has been exposed to thereby enabling assessment of bio-available toxicant levels in the environment.
it is our belief that the microarray technology coupled with regression modeling provides an innovative and promising tool towards this direction.
the natural next steps are to demonstrate the applicability and prediction power of this new approach in scenarios of contaminants mixtures and also at field contaminated military sites.
ultimately, information collected from such studies will be used to support further development of predictive modeling for toxicogenomic measures of exposure  <cit> .
availability of supporting data
the microarray datasets were deposited in the national center for biotechnology information ’s gene expression omnibus  database as series gse <dig>  and gse <dig> .
the custom-designed earthworm  microarray with 15 k oligo probes is accessible as geo platform gpl <dig>
additional files
additional file 1: description of  <dig> regression models and their references.
additional file 2: identified de genes , their expression data and measured tissue residue results for 4- and 14-day tnt-exposed earthworms.
additional file 3: identified de genes , their expression data and measured tissue residue results for 4- and 14-day rdx-exposed earthworms.
additional file 4: identified de genes and other selected predictor genes, measured tissue residue of hmx, and the performance of  <dig> regression models using only de genes as predictors for the 4-, 14- and 28-day exposures.
additional file 5: functional annotation of  <dig> identified predictor genes for tnt, rdx and hmx residues using blast2go.
abbreviations
dedifferentially expressed
fdrfalse detection rate
geogene expression omnibus
gvsgibbs variable selection
hmxhigh melting explosive
hplc-rfdhigh performance liquid chromatography equipped with a radioactivity flow detector
mcmcmarkov chain monte carlo
ncbinational center for biotechnology information
nrcnational research council
rdxroyal department formula x
ssvsstochastic search variable selection
svrsupport vector regression
tnt <dig> ,6-trinitrotoluene
competing interest
the authors declare that they have no competing interests.
authors’ contributions
pg and ejp conceived the study.
reb, drj and bcs conducted earthworm exposure, chemical analyses and tissue residue data analysis.
ndb and pg performed microarray experiments and statistical analysis of microarray data.
xn, yc and dew carried out regression modeling.
pg coordinated the study and drafted the manuscript.
all authors participated in the revision of the manuscript and read and approved the final manuscript.
this work was funded by the u.s. army environmental quality/installation technology research program as part of projects “toxicogenomics for assessment of munitions and explosives related compounds” and “impact of munitions constituents on biological networks”.
drs.
x. nan, y. chen and d. wilkins were supported in part by the u.s. national science foundation under award number eps  <dig>  we thank dr. guilherme r. lotufo of u.s. army corps of engineers for his suggestion of using scintigest solubilizer for tissue digestion, and a. e. elaine van ommen kloeke of vu university amsterdam for sharing blast2go annotations of earthworm transcriptome.
permission to publish this information was granted by the chief of engineers.

 bioavailability processes were defined in a national research council report as the individual physical, chemical, and biological interactions that determine the exposure of plants and animals to chemicals associated with soils and sediments. in environmental risk assessment, the amount of chemicals taken up by an animal or plant is termed dose, or, interchangeably tissue residue, body burden or chemical bioavailability. a distinction exists between dose and exposure as the latter is defined as the amount of chemicals present in the immediate environment where the organism is exposed to. soils are a major sink for many environmental contaminants including explosives compounds. soil contamination by military unique compounds is a serious environmental concern that can result in the formation of chemical residues in tissue of exposed organisms. since many site-specific biotic and abiotic factors can modify the form, mobility and availability of these contaminants, the actual exposure risk to ecological receptors may be less than that suggested by their total concentration. therefore, the extent to which chemicals are bio-available has significant implications for risk management and remedial decision-making at contaminated sites. a myriad of biological, physical and chemical approaches have been used to evaluate bioavailability of chemicals. for instance, biological approaches include techniques measuring organismal uptake, response, and toxicity. tissue residues or body burden of chemicals, often determined using analytical chemistry methods, can represent the bio-available or biologically effective concentration at the target site. however, little consensus exists about optimal approaches for measuring bioavailability. as summarized in the aforementioned nrc report and a recent review, no single tool is free from limitations and none can be applied universally. an intensive effort to develop mechanistic tools or models based on mechanisms is critical to future development of bioavailability tools. recent advances have prompted the application of genomics-based technologies to ecological risk assessment, including screening, tiered testing, monitoring, remediation, and regulatory decision-making. toxicogenomics has been extensively applied to assess toxicological effects, especially in biomarker discovery and toxicity mechanistic investigations. on the other hand, toxicogenomics also possesses a great potential for providing a quantitative measure of chemical exposure. nevertheless, the application of toxicogenomics to exposure assessment has been relatively under-explored. previously, we applied microarray technology to profile gene expression in earthworms exposed to explosives,6-trinitrotoluene or,5-trinitro-,5-triazacyclohexane for or 14 days, and built multivariate regression models to quantitatively predict earthworm tissue residue of these two compounds. the models, however, showed only a modest predictive power, explaining close to half of the variance for tnt tissue residue and one-quarter of the variance for rdx tissue residue. the present study was motivated to improve the quantitative predictive power of regression modeling based on genome-wide gene expression data. our hypothesis was that small sets of predictor genes could be identified and used to build multivariate regression models for quantitative prediction of tissue residue levels of explosives compounds. consistent with our previous study, the overall goal of this study was to investigate the feasibility of using gene expression data to assess animal exposure. to achieve this goal, we re-analyzed the previous datasets by separating the 4-day exposure from the 14-day exposure, and expanding from two regression methods to methods. in addition, we have also generated a new dataset from hmx-exposed earthworms and applied the same approach to analyzing it. our results demonstrated that much higher prediction accuracies were attained indicating that microarray gene expression coupled with multivariate regression modeling is a viable approach for assessing chemical bioavailability. methods the overall experimental approach is depicted in fig.  we first exposed earthworms to explosives compounds and measured tissue residues and gene expression profiles in both exposed earthworms and unexposed controls. then we identified small sets of predictor genes to build regression models, where y’ is the predicted residue and β is a coefficient vector. these models were applied to predict the bioavailability in exposed worms using only predictor gene expression data. finally, the actual measured tissue residues were compared with the predicted values to determine prediction accuracy evaluated by correlation coefficients. our ultimate goal was to apply these prediction models to estimate the unknown bioavailability in worms from their gene expression data.fig. the overall experimental approach. see the methods section for explanation experimental design previously, we collected a large 248-sample dataset with measurements of transcriptome-wide gene expression and tissue residue of two explosives chemicals, -labeled tnt and rdx, in the earthworm eisenia fetida. briefly, three sets of experiments were conducted by exposing mature adult worms bearing a clear clitellum in a pristine sandy loam soil amended with tnt or rdx. nominal exposure concentrations were as follows: 96 mg tnt/kg, or 128 mg rdx/kg for 4 days or 14 days. the 4-day exposure was repeated with the same tnt concentrations but different rdx concentrations. following the same experimental design, we exposed adult worms to -labeled hmx for and 28 days at nominal concentrations of, and 128 mg/kg in the same silt loam soil as in our previous study. hmx was purchased from dupont nen with an initial specific activity of  mci/mmol. earthworms were reared in-house in a continuous lab culture as previously described. ten worms were exposed in 250 g of soil per each treatment. the nominal soil concentrations were verified to possess less than 10 % variations from target concentrations using a high performance liquid chromatography equipped with a radioactivity flow detector. upon termination, all worms were flash-frozen in liquid nitrogen, stored at−80 °c, lyophilized at −40 °c, and homogenized. tissue residue analysis for hmx all ten worms per treatment were analyzed for tissue residue. for tissue hmx-residue analysis, triplicate subsamples were digested in 1 ml of  n scintigest tissue solubilizer for 18 h. radioactivity was determined by adding 200 μl of tissue digest to ultima gold scintillation cocktail using a packard tricarb 2500tr liquid scintillation counter. prior to the analysis, we validated this method by comparing it with our previously used methods. testing results showed that the solubilization method yielded 92 % ± 3 % of the radioactivity measured previously by using oxidization methods for the same tnt- or rdx-exposed worms. gene expression profiling for hmx-exposed worms total rna was extracted from five of the ten worms per treatment using an rneasy mini kit. each rna sample was hybridized to the custom-designed agilent 15 k e. fetida oligo array previously used for generating the tnt and rdx microarray datasets. this array contained unique transcript-targeted 60-mer oligo probes. details of array hybridization, gene expression data acquisition and pre-processing were described elsewhere. identification of differentially expressed genes de genes were identified for all three datasets among multiple concentrations using a multivariate permutation random-variance t-test or f-test implemented in brb-arraytools version.. a gene was considered statistically significant if it achieved 80 % confidence that the false discovery rate was less than 10 %. de genes were derived for individual exposure duration and explosives compound separately. regression prediction modeling for residue prediction, we chose different multivariate regression models and employed double-looped, 10-fold cross-validation as described in statnikov et al. to assess prediction accuracy. the inner loop was used to determine the optimal value of parameters for training in the outer loop. model performance was estimated in the outer loop by training on all splits but one, and using the remaining one for testing. coefficient of determination was calculated using pearson’s formula to describe prediction accuracy or “goodness of fit”, i.e., how well the tissue residue predicted by a regression model represented the actually measured residue of a worm tissue sample. the models include six linear models ) and nonlinear models. matlab codes were scripted to implement array data preprocessing, regression, and cross-validation, and are available upon request.  selection of predictor genes from tnt and rdx gene expression datasets we reanalyzed the tnt and rdx microarray datasets by separating the 4-day sample set from the 14-day sample set. differentially expressed genes were identified among different treatments of 4-day or 14-day earthworm samples using a multivariate permutation test. statistical reanalysis of the 14-day exposure gene expression data resulted in six and three de genes for tnt and rdx, respectively. for the 4-day tnt exposure, and genes were inferred as de genes from the original exposure, the repeat exposure and the original vs. repeat controls, respectively. a group of de genes were found to be common between the original and the repeat tnt exposures. this group was further reduced to genes after excluding genes also found to be significantly altered between the original and the repeat controls. similarly, and de genes were derived from the original rdx and the repeat rdx exposures, respectively, with genes in common. twenty-six genes out of the 178-gene group remained after excluding the same de genes appearing in the controls comparison. the final sets of identified de genes are shown in additional file and additional file the low numbers of de genes found in 14-day exposures are consistent with our previous report. we obtained more de genes than previously for the 4-day exposures because of the reduced statistics stringency. however, worms used in the original and the repeat exposures exhibited significant differences that were reflected as nearly 1/ of all 15 k profiled genes differentially expressed between the two control groups. therefore, we chose to remove these genes not responding specifically to 4-day tnt or rdx exposure from the final de gene lists. regression predictive modeling for tnt and rdx tissue residue using the aforesaid final sets of de genes as predictor genes, the regression methods displayed varied power in predicting worm tissue residues of tnt and rdx. the predictive power was assessed using the coefficient of determination as a measure of the accuracy of the data model. no single method was placed as the best performer for all four datasets. for instance, the multivariate linear regression model was the best performer for the 14-day tnt exposure dataset, but its performance was relatively weak on the 4-day tnt exposure dataset. quite a few models such as lasso, elastic net, the ridge family models, the svr family models except for svr sigmoid, and nadaraya-watson performed consistently well across all four datasets. in contrast, some models like svr sigmoid, loglog, and reverse regression performed poorly or even appeared inapplicable to the datasets.table performance of regression modeling methods on four datasets assessed by coefficient of determination estimated from ten runs of 10-fold cross-validation with values of the best performing method for each dataset shown in bold  ±   ±   ±   ±  rdx_d 4-day rdx exposure, rdx_d 14-day rdx exposure, tnt_d 4-day tnt exposure, tnt_d 14-day tnt exposure, na not available. see additional file for the lists and annotation of predictor genes the predictive power of the best performers was remarkably improved in comparison to our previously published results. on average, these models explained 75 % or 66 % variance of the 4-day samples, and 72 % or 71 % variance of the 14-day samples. in our previous study, the best prediction models explained roughly one-quarter and less than one-half of the variance for rdx and tnt, respectively.fig. the average predicted versus the measured tissue residues of tnt or rdx in all 4-day or 14-day exposed samples using their respective best performing models tissue residue in hmx-exposed earthworms the measured tissue residue of radio-labeled hmx increased with the increasing nominal amendment concentration and also with the duration of exposure. trace amounts of hmx were detected in some of the two control groups, which are equivalent to the background noise level or the lower detection limit of the analytical method. the worm tissue residue did not appear to plateau at the highest amendment level, suggesting that more hmx could be taken up from soil by the earthworm, given a higher amendment concentration.fig. tissue residue of radio-labeled hmx measured in earthworms exposed for 4-, 14-, and 28-days. data are represented as mean  + standard deviation with n =  bc = blank control; sc = solvent control predictor genes selection for hmx gene expression dataset the hmx microarray gene expression dataset consists of three exposure groups of arrays, i.e., treatments  ×  replicates. the two control groups in all three exposures showed no significant difference because no de gene was identifiable at the same settings as above. so, the two groups were combined as one control group for subsequent statistical analyses. the three day samplings exhibited little difference with only seven derived de genes, suggesting that worms used in all three exposures were nearly identical and that animal batches had a minimal impact on the test results. only and six genes were inferred as de genes for the 4-, 14-, and 28-day exposures, respectively. there was no overlap among all four sets of inferred de genes. using de genes as predictors, regression modeling resulted in a modest predictive power for the 4-day and 14-day exposures with the best performer only explaining an average of 55 % variance. performance of the predictor genes on the 28-day dataset was much better with the highest r of. all these preliminary results are provided in additional file in order to improve the predictive power, a different approach was used to identify additional predictor genes. we ran a correlation test to determine the degree of dependence between each expressed gene and the tissue residue. then, all 15 k genes were ranked according to their coefficients of correlation. with a cutoff of |r| = , we obtained and most correlated genes for the 4-d, 14-d and 28-d exposures, respectively, which were considered potential predictor genes. genes at the top of the most correlated genes lists were most positively correlated with tissue residue, while those at the bottom most negatively correlated with tissue residue. some of the genes on these lists were also de genes. we conducted performance analyses of the incremental new sets of predictor genes by adding two or three genes from the top or the bottom of the most correlated genes list at a time to the de genes to form a new predictor gene set. results indicate that addition of the top two and the bottom two had an optimal enhancement in prediction accuracy. therefore, the final revised sets of predictor genes consisted of de genes, the top two and the bottom two genes on the most correlated genes lists. prediction outcome for hmx gene expression dataset using the revised sets of predictor genes, the regression models achieved substantial enhancement in prediction accuracy. particularly, prediction outcomes for the 4- and 14-day exposures improved by nearly 20 % as the best performers for the two datasets both explained 72–73 % of the variance. the 28-day dataset did not improve as much largely because five out of the six de genes were among the most correlated genes list with one being the most negatively correlated gene and another the third most positively correlated gene. in contrast, none of the four de genes derived from the 4- and 14-day datasets were on the most correlated genes lists. models that performed well on the tnt and rdx datasets also did consistently well on the hmx datasets. three models, svr sigmoid, inverse regression and loglog regression, were not suitable for all datasets, probably because of transformation and normalization operations in data pre-processing.table performance of regression modeling methods on the three hmx exposure datasets assessed by coefficient of determination estimated from ten runs of 10-fold cross-validation with values of the best performing method shown in bold  ±   ±   ±  d 4-day hmx exposure, d 14-day hmx exposure, d 28-day hmx exposure, na not available. see additional file for the lists and annotation of predictor genes fig. prediction results of 4-, 14- and 28-day hmx-exposed earthworm tissue residues using the best performing models discussion quantitative prediction of a phenotype or trait using high dimensional gene expression data has been long explored in many research fields such as human diseases, animal breeding and plant genetics. the phenotypes of prediction interest are broad and have included thrombocytosis etiologies, acute myeloid leukemia resistance and breast cancer tumor response to chemotherapy in biomedical research, cattle milk yield in animal breeding, and wheat grain yield and plant pathogen infection severity in plant genetics. the methodology of predictive modeling varies widely from bayesian network-based approaches to bayesian hierarchical regression modeling and ordinary multivariate regression. the reported prediction accuracy also varies from one study to another, implying significant challenges and opportunities co-existing in this field. similar to the aforementioned fields of research, we have also sought to use transcriptomic data coupled with regression modeling to predict a phenotype, i.e., chemical residues in animals. although prediction accuracies were unsatisfactory in our previous study, we have significantly improved them in the present study without adding to the complexity of regression models. the size of predictor gene sets has been reduced to three to ten genes for all datasets except the 4-day tnt and rdx exposures, which are likely the lowest number of predictors possible, given the complicated processes involved in residue formation. although there is still room for further improvement, the accuracies obtained here are comparable to those reported in other quantitative trait prediction studies. it has been reported that variable or feature selection is a key step towards accurate prediction. to select features that have high predictive power as predictors, many different strategies have been developed, including t- or f-test based approaches to identify differentially expressed genes, correlation based methods to rank genes, machine learning based techniques, and some combinations of the three types of strategies. in the present study, we have used both f-test and correlation based approaches to select for subsets of predictor genes. a variety of regression models were applied to describe relationships between the response variable and a set of predictors through a regression function involving some parameter to be estimated from data. given the time-dependent nature of tissue residue and gene expression levels, the datasets were split by exposure duration. results indicate that this step remarkably improved the prediction outcome for both tnt and rdx exposures. to be used as predictors, altered levels of a given transcript need not be connected to a specific biologic end point, nor do the specific functions of all the mrna molecules have to be known, although such information would be valuable. bioinformatic annotation of all selected earthworm predictor genes using blast2go shows that only of them have meaningful biological functions and that their contribution to residue formation remain largely unknown. recently, there has been an intense interest in performing variable selection simultaneously with parameter estimation in predictive modeling. frequently applied approaches include bayesian methods, stochastic search variable selection, adaptive shrinkage with jeffreys’ prior or a laplacian prior, and reversible jump markov chain monte carlo ) and lasso, ridge or elastic net regularization methods. given such a large variety, there is no consensus with regard to what method is the universally best performer. for instance, o’hara and sillanpaa tested several bayesian variable selection methods on both simulated and experimentally collected data and concluded that ssvs, reversible jump mcmc and adaptive shrinkage methods all worked well, but the choice of which method was better depended on the priors that were used, and also on how they were implemented. similarly, fu et al. concluded that svr, partial least squares regression and multiple linear regression yielded higher prediction accuracies for one dataset but transcriptome-based distances worked better on another dataset. in the present study, we also observed that lasso, ridge and elastic net regularization methods performed almost equally well on all datasets, whereas inverse and loglog regression methods performed poorly on the datasets in a consistent fashion. it is worth noting that the approach employed in this study faces several challenges when applied to residue prediction. first, compared with other existing approaches, it requires a more extensive effort to identify and optimize a set of predictor genes. second, if the species of interest does not have a transcriptome-wide microarray available, one has to either design the microarray from scratch or use the array of a closely related species, which may limit its applicability. third, as a novel approach, there is no doubt that it is still in its embryonic phase and that its full potential and limitations are not yet explored thoroughly. for instance, our approach may be applied to quantitative predictions in drug discovery such as predictive ranking of new drug toxicity and/or potency as they share similarities in high data dimensionality and transcriptomic profiling based on either microarray- or rna-seq/next-generation sequencing-based data.  chemical residue in exposed animals is an important dose metric in environmental risk assessment. the formation of tissue residue involves a number of complex biological processes which can be reflected as expression profiles in microarray experiments. building predictive models of tissue residue based on gene expression would help to accurately assess how much chemical an animal has been exposed to thereby enabling assessment of bio-available toxicant levels in the environment. it is our belief that the microarray technology coupled with regression modeling provides an innovative and promising tool towards this direction. the natural next steps are to demonstrate the applicability and prediction power of this new approach in scenarios of contaminants mixtures and also at field contaminated military sites. ultimately, information collected from such studies will be used to support further development of predictive modeling for toxicogenomic measures of exposure. availability of supporting data the microarray datasets were deposited in the national center for biotechnology information ’s gene expression omnibus database as series gse and gse. the custom-designed earthworm microarray with 15 k oligo probes is accessible as geo platform gpl additional files additional file 1: description of regression models and their references. additional file 2: identified de genes, their expression data and measured tissue residue results for 4- and 14-day tnt-exposed earthworms. additional file 3: identified de genes, their expression data and measured tissue residue results for 4- and 14-day rdx-exposed earthworms. additional file 4: identified de genes and other selected predictor genes, measured tissue residue of hmx, and the performance of regression models using only de genes as predictors for the 4-, 14- and 28-day exposures. additional file 5: functional annotation of identified predictor genes for tnt, rdx and hmx residues using blast2go. abbreviations dedifferentially expressed fdrfalse detection rate geogene expression omnibus gvsgibbs variable selection hmxhigh melting explosive hplc-rfdhigh performance liquid chromatography equipped with a radioactivity flow detector mcmcmarkov chain monte carlo ncbinational center for biotechnology information nrcnational research council rdxroyal department formula x ssvsstochastic search variable selection svrsupport vector regression tnt,6-trinitrotoluene competing interest the authors declare that they have no competing interests. authors’ contributions pg and ejp conceived the study. reb, drj and bcs conducted earthworm exposure, chemical analyses and tissue residue data analysis. ndb and pg performed microarray experiments and statistical analysis of microarray data. xn, yc and dew carried out regression modeling. pg coordinated the study and drafted the manuscript. all authors participated in the revision of the manuscript and read and approved the final manuscript. this work was funded by the u.s. army environmental quality/installation technology research program as part of projects “toxicogenomics for assessment of munitions and explosives related compounds” and “impact of munitions constituents on biological networks”. drs. x. nan, y. chen and d. wilkins were supported in part by the u.s. national science foundation under award number eps we thank dr. guilherme r. lotufo of u.s. army corps of engineers for his suggestion of using scintigest solubilizer for tissue digestion, and a. e. elaine van ommen kloeke of vu university amsterdam for sharing blast2go annotations of earthworm transcriptome. permission to publish this information was granted by the chief of engineers.
chemical bioavailability is an important dose metric in environmental risk assessment. our results demonstrate that a prediction accuracy of r2 =  – was achievable at a relatively low model complexity with as few as 3– predictor genes per model. in the present study, we further compared different regression models and performed variable selection simultaneously with parameter estimation. these results are much more encouraging than our previous ones. previously, we developed a new genomics-based approach that integrated microarray technology and regression modeling for predicting bioavailability of explosives compounds in exposed earthworms.
 chemical bioavailability is an important dose metric in environmental risk assessment. although many approaches have been used to evaluate bioavailability, not a single approach is free from limitations. previously, we developed a new genomics-based approach that integrated microarray technology and regression modeling for predicting bioavailability of explosives compounds in exposed earthworms. in the present study, we further compared different regression models and performed variable selection simultaneously with parameter estimation.  this refined approach was applied to both previously collected and newly acquired earthworm microarray gene expression datasets for three explosive compounds. our results demonstrate that a prediction accuracy of r2 =  – was achievable at a relatively low model complexity with as few as 3– predictor genes per model. these results are much more encouraging than our previous ones.  this study has demonstrated that our approach is promising for bioavailability measurement, which warrants further studies of mixed contamination scenarios in field settings electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords tissue residueglobal gene expression profilingpredictor genespredictive regression modelingtnt rdx hmx earthworm http://dx.doi.org/ /100000001national science foundation eps0903787u.s. army environmental quality/installation technology research programhttp://dx.doi.org/ /100006435center for selective c-h functionalization, national science foundation eps0903787chen yixin u.s. army enviornmental quality/installation technology research programissue-copyright-statement© the author 2016 