BACKGROUND
in order to understand how quantitative variation in phenotypic traits is related to the underlying genetic differences between plants, and to differences in gene expression, protein constitution and metabolic variability, an approach is needed in which the combined molecular signature of the plants is shown to be predictive for the phenotypic traits of interest .
we can use high-throughput ~ omics technologies, such as gene expression microarrays  <cit> , mass spectrometry   <cit>  and protein chips  <cit>  to obtain molecular signatures of a population of plants.
in addition we can study phenotypic differences in the same population and hypothesize that differences in the phenotypic trait in the population are related to the variation in these combined molecular profiles across the different data sets  <cit> .
finding the ~ omics variables that are related to the phenotypic traits of interest can then be used in two ways: 1) for prediction of the traits from these molecular profiles.
for example steinfath et al.,  <dig>  <cit>  explored metabolomics data to predict agronomic phenotypes of potato crop plants grown in different environments.
they identified metabolites that can be used as biomarkers for these traits and hence to improve the selection on these traits which can be implemented in breeding programs.
2) to identify functional relationships between traits and molecular networks of the plants .
this means that we are not just interested in interrelating these ~ omics data sets for their own sake to find genetic and metabolic networks, but the networks and their elements should actually be predictive for a phenotypic trait of interest.
on the other hand, in the context of genetics and plant breeding, we also want to be protected against finding relationships between phenotype and ~ omics data that are just caused by environmental or developmental differences.
instead, we are interested to find relationships that have a basis in the genetic differences between plants.
therefore a mapping population is an ideal target for this kind of study: we can study whether an observed relationship between a phenotypic trait and ~ omics variables is also based on genetic differences between the plants in the segregating population, as we can actually map the phenotypic variation as well as the variation in the ~ omics data sets.
a relationship that would be just based on variation in environmental influences or conditions would not result in mapped qtls for the phenotypic traits or the ~ omics variables  <cit> .
we are interested in the association between a number of phenotypic traits related to tuber quality of potato and several ~ omics data sets, some of which were published previously by us .
in addition, we use mapping and genotyping information since the population that we use is a mapping population.
the quality traits considered are 1) potato tuber flesh colour, 2) enzymatic discoloration after peeling, 3) starch gelatinization as measured by differential scanning calorimetry  and 4) tuber shape.
out of these four phenotypes, starch gelatinization as measured by differential scanning calorimetry  and tuber shape phenotypes were not analyzed before.
the other two phenotypes were included to make the link with the previous papers and to permit validation and accuracy testing of our used network analyses.
in this study, we relate transcriptomics, metabolomics and proteomics data to genetic variation in these quality traits in the mapping population.
for this, we use a novel strategy consisting of three steps.
first, using the same approach as in an earlier study  <cit> , we apply random forest regression  to find, per single trait and per individual ~ omics data set, the variables that play a significant role in the prediction of each of the quality traits.
in the second step, we use qtl mapping of the quality traits and of the ~ omics variables to select variables that have a qtl  cosegregating with a quality trait qtl and to remove redundancy in the set of selected variables.
finally, we construct regularized partial correlation networks, for each of the quality traits, of the selected sets of variables over the four ~ omics data sets that have a genetic association with the traits per qtl.
because already much is known about the regulatory genetic and metabolic pathways involved in tuber flesh colour , we used this trait to validate the approach, as we had done for individual data sets in  <cit> .
here we analyze all ~ omics data sets simultaneously to construct networks related to phenotypic traits consisting of features across all these data sets.
we demonstrate that also for the other phenotypic traits we can find networks of small sets of interrelated gene expression profiles, proteins and metabolites that are associated with and predictive for these quality traits.
methods
plant material
we used  <dig> individuals, including the parental clones, of a diploid potato backcross population   <cit> .
this population is derived from an original cross between potato clones c  and e  and is described in detail in  <cit> .
all clones were grown in multi-year repeats in the field, wageningen, the netherlands during the normal potato-growing season in the netherlands   <cit> .
for each genotype, tubers were collected from three plants and representative samples were either used for phenotypic analysis or mechanically peeled and immediately frozen in liquid nitrogen before being ground into a fine powder and stored at −80 °c for metabolomics, transcriptomic and proteomic analyses.
the determination of carotenoids was as described in  <cit> .
this was a targeted metabolic analysis which includes compounds like zeaxanthin, violaxanthin and a compound for which the chemical behaviour is like violaxanthin .
differential scanning calorimetry  is extensively used to study physical properties of starch granules.
dsc-onset measurements  report water temperatures at which starch granules reach their gelatinization state  <cit> .
dsc measurements provided gelatinization onset temperatures of starch granules for  <dig> individuals of the c x e population ranging between  <dig>  °c and  <dig>  °c.
enzymatic discoloration of tubers after peeling and being exposed to air  at different time points such as after 5 min, 30 min, 3 h and difference in discoloration between 3 h and 30 min and 3 h is described in  <cit> .
in this study, we considered only enzymatic discoloration after 5 min since the measurements at later time points are all highly correlated.
potato tuber shape was scored between  <dig>  and  <dig>   <cit> .
omics data sets
we generated different types of ~ omics data to integrate and find novel relationships among them.
parts of the ~ omics data sets were already published individually, for example, the transcriptomics data  <cit> , gc-ms data  <cit> , lc-ms data  <cit> .
only the proteomics data was generated more recently  and these data have not been considered for a combined analysis over all ~ omics technologies simultaneously.
in this section, we briefly summarize the different omics data generation methods and refer to  <cit>  for further details.
the phenotypic data and ~ omics data are provided in a supplementary file .
transcriptomic data
rna was extracted from the  <dig> samples using the hot phenol method described previously  <cit> .
according to the manufacturer’s protocol starting with 2 μg of purified total rna  <cit> ).
for additional data analyses only genes with a pearson correlation coefficient higher than  <dig>  between the cy <dig> and cy <dig> datasets were included resulting in  <dig>  expressed genes.
we took into account only the cy <dig> gene expression signals for further statistical analysis.
for visualization, we used the gene nomenclature in the following way: gene_gene id .
the number refers to the gene id of the supplementary material of  <cit> .
lc-ms, data generation, processing and identification
potato tuber samples were analyzed for variation in semi-polar metabolite composition using an untargeted accurate mass lc-ms approach.
in total  <dig>  mass signals were obtained from mass spectrometry and selected based on the following criteria: signals should be present in at least  <dig> samples and with at least one amplitude higher than  <dig> units .
finally, data redundancy of signals derived from the same compound, i.e.
isotopes, adducts and in-source fragments, was removed by by retention time-dependent clustering.
this clustering of the  <dig> signals resulted in  <dig> reconstructed metabolites .
the untargeted metabolites are represented as centrotype_mass_scan number as in  <cit> .
as an example: 818_795_ <dig> means that the centrotype number is  <dig>  the mass number  <dig> and the scan number  <dig>  for visualization and simplicity we used lc-centrotype number .
for more information on data generation, processing and identification see the methods of  <cit> .
gc-ms data generation, processing and identification
gc-ms data were generated from the same  <dig> genotypes of the cxe population.
detailed materials and methods for gc-ms data generation, processing and identification are as in  <cit> .
in short, raw data were processed by chromatof software  <dig>   and masslynx software , and further analysis was performed using metalign software to extract and align the mass signals.
mass signals that were present in fewer than two samples were discarded.
signal redundancy per metabolite was removed by means of clustering, and mass spectra were reconstructed  <cit> .
this resulted in  <dig> reconstructed polar metabolites .
compounds were subjected to tentative identification by matching to the nist <dig> and wiley spectral libraries.
library hits were manually curated, and a series of commercial standards was used to check annotation.
for visualization and simplicity we used gc-centrotype number.
for example, 7275_20721700_ <dig> means that the centrotype number is  <dig>  the mass number  <dig> and the scan number  <dig>  for visualization only gc_ <dig> is used.
proteomics data generation, processing and identification
protein extraction
total protein was extracted from approximately  <dig>  g of ground tuber material, to which 1 ml of pre-heated  lysis buffer , sds , dtt , pvp-p ) was added.
samples were homogenized and protein amount was measured using the rc/dc assay .
protein labelling
a single lysine per protein molecule was labelled using the fluorescent cydyes from the difference gel electrophoresis  technology  according to the manufacturer’s protocol.
the internal standard was labelled with cy <dig> and consists of an equal mixture of  <dig> randomly chosen samples of the experiment .
every 2d-gel contains one sample labelled with cy <dig>  one labelled with cy <dig> and the internal standard labelled with cy <dig>  this means that every sample on each gel can be compared by using the internal standard sample labelled with cy <dig>
2d-electrophoresis
electrophoresis was performed using first dimension electrophoresis on 24 cm immobilized ph gradient strips  with a linear ph range from  <dig> to  <dig> on an ettan ipgphor isoelectric focusing  system.
cydye labelled samples  were loaded to the strips and focusing was run for 18 h at 20 °c with the following settings: 3 h 150 v, 3 h 300 v, from 300 v to 1000 v in 6 h, from 1000 v to  <dig>  v in 1 h and finally 5 h at  <dig>  v. after equilibration in the dark at room temperature in equilibration buffer , sds 2 % ) containing dtt 1 % the second dimension electrophoresis was run on the ettan dalt  <dig> system on precast  <dig>  % sds polyacrylamide slab gel  and buffers from ge healthcare/amersham biosciences.
the separated cydye-labelled proteins were visualized by scanning with a ettan dige imager , using for cy <dig> an 480 nm laser and an emission filter of 530 nm, for cy <dig> an 540 nm laser and an emission filter of 595 nm and for cy <dig> an 635 nm laser and an emission filter of 680 nm.
gel images were analysed with the decyder software version  <dig> according to decyder 2dv <dig> manual; ge healthcare/amersham biosciences.
detected spots were filtered based on spot volume larger than relative value  <dig>  to exclude spots that could be just background noise or dust particles.
the internal standard in each gel was used to automatically match all images to the reference  value was used.
each of the proteins is presented by pro_x where “x” represents consecutive protein numbers.
protein identification
spots of interest were excised from gel using the ettan spot picker.
peptide mass determinations were carried out using the applied biosystems  <dig> proteomics analyzer.
both pmf and ms/ms in reflectron mode analyses were carried out with the samples.
calibration was carried out with a peptide mass calibration kit.
proteins were identified by searching against the ncbi ‘viridiplantae’ and an est ‘viridiplantae-eudicots’ database using mascot.
random forest regression
random forest  <cit>  was used for regression of the phenotypic trait such as flesh colour, tuber shape, dsc onset, enzymatic discoloration on the transcriptomics, gc-tof-ms signals, lc-ms signals, 2d-dige proteomics.
rf constructs a predictive model for the response using all predictors but quantifies the importance of each, here the metabolites, gene expression and proteins in explaining the variation present in the trait.
rf by itself does not provide significance levels of individual metabolites and does not perform a variable selection to choose a predictive subset of associated metabolites.
therefore, we included a permutation test to indicate significance of the association of a metabolite with a trait.
in each of  <dig> permutations of the trait values we estimated the variance explained by the rf model  and the variable importance of each metabolite in terms of the decrease in node impurities  <cit> .
we ordered node purity values from the permuted data sets and took the  <dig> percentile from the distribution of impurity values as the significance threshold of the individual metabolites.
the same procedure was done for r <dig> values of the model: the  <dig> percentile was taken as a significance threshold for the rf model.
r <dig> in random forest is not just a measure of goodness-of-fit of the data at hand but is determined on left-out samples  so it should be interpreted as a measure for predictive quality  of the random forest on independent samples that have the same properties as the in-bag sample  <cit> .
random forest analysis was done using r statistical software using the randomforest package.
qtl analysis and cis- and trans-eqtls
we mapped expression qtls  from the gene expression data, metabolite qtls  both for lc-ms and gc-ms, and protein qtls  from the proteomics data to find regions on the genome explaining genetic variation in gene expression, metabolite and protein values using the integrated linkage map of the c and e parents for qtl analysis  <cit> .
further, we used the potato genome physical map  <cit>  to investigate eqtl physical positions of genes identified as predictive for phenotypic traits from the random forest analyses.
the potato oligo  microarray  used in the experiments contains  <dig>  features based on a potato unigene set  <cit> .
to allow discrimination between cis- and trans-eqtls all unigenes were blasted against the genome scaffold sequences, predicted coding sequences  and predicted gene regions .
features with a unique and significant hit were assigned to genome scaffolds for which the majority has chromosome information.
identified qtls on the same linkage group as their physical map position are identified as cis-acting while qtls on different linkage groups are defined as trans-acting.
features on the array for which no physical map position could be assigned are classified as unknown  <cit> .
all qtl analyses were done using metanetwork  <cit>  package in r.
network reconstruction
from the rf analyses, for each of the traits, we obtain a list of significant genes, metabolites and proteins.
for visualization in an association network, we reduced the number of genes by checking chromosomal positions of each of the predicted significant genes, metabolites  and proteins.
we took the most significant single gene, lc-ms signal, gc-ms signal and/or protein per chromosome as representative for that ~ omics data set and chromosome and made a full order regularized partial correlation network with these genes, metabolites, proteins and traits as nodes in the network  <cit>  and the strength of the interaction as edges in the network.
partial correlation measures the correlation between two variables after their linear dependence on other variables is removed.
it can distinguish between direct and indirect associations whereas correlation-based network cannot and often yield many spurious edges .
a detailed description of regularized partial correlation is provided in  <cit> .
RESULTS
selection of ~ omics features predictive for quality traits
the prediction of variation in potato tuber flesh colour  <cit>  from the gene expression, lc-ms, and protein data sets was quite high , but much lower for the gc-ms data .
for flesh colour, the microarray data and the lc-ms data were equally good for prediction in terms of the explained variance, but the numbers of significant features were very different: the prediction of flesh colour using  <dig> significant lc-ms features is almost the same as for  <dig> significant gene expression profiles  from the microarray data set.
from the gene expression data, the gene which ranks first and third with respect to variable importance for predicting flesh colour was a beta-carotene hydroxylase .
two oligos were present on the array targeting the same beta-carotene hydroxylase gene, hence the two high ranks for the same gene.
another gene from the carotenoid pathway, zeaxanthin epoxidase  ranked forty-fourth.
based on our current knowledge of potato tuber flesh colour and carotenoid content  <cit> , these two genes were expected to be associated with flesh colour.
from the gc-ms data, out of six significant metabolites, four had an annotation: malic acid, 2-4-5-trihydroxypentanoic acid, glucopyranose, 2-butenedioic acid- and bis ester.table  <dig> percentage variance explained  in out-of-bag  prediction by random forest  models using all genes, lc-peaks, gc-peaks or proteins separately
non-significant models are indicated as ns
gene + lc + gc + protein;
gene + gc + protein;
gene + lc + protein;
the numbers of significant  genes, lc-peaks, gc-peaks or proteins are between brackets.
non-significant models are indicated as ns
for the three other traits the gene expression data was more predictive than the metabolomics and proteomics data sets, and for tuber shape actually no significant features were found for gc-ms, lc-ms and the protein data, while significant gene expression differences explain 55 % of the variation in tuber shape.
for starch gelatinization  no lc-ms features were significant, for enzymatic discoloration no gc-ms features were significant.
associated genomic regions
using random forest regression, we generated lists of candidate genes, metabolite features and proteins that are predictive for quality traits of potato tubers .
however, from this prediction we cannot conclude that these associations necessarily have a basis in genetic differences between the plants since they could also be caused by environmental variation in both the trait and the levels of the features, or in developmental differences.
therefore we also investigated qtl positions for both the phenotypic traits and the ~ omics features.
when co-segregation of an eqtl, mqtl or pqtl with a trait qtl is observed, this could imply a functional relationship or identify causal genes/proteins/metabolites.
however, this is not necessarily the case since any linked but functionally unrelated qtls will also show this co-segregation.
still, these sets of predictive genes highlight genomic regions of interest, comparable to standard qtl analysis, and can provide additional information or help to narrow down the region of the causative polymorphism.
for tuber flesh colour the  <dig> eqtls of significantly predictive gene expression profiles from rf mapped to eight chromosomes:  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>  a large number of these  eqtls were mapped to chromosomes  <dig> and  <dig>  significant gc-ms peaks mapped to two chromosomes,  <dig> and  <dig>  qtls for significant lc-ms peaks mapped to chromosomes  <dig> and  <dig> and qtls for significant protein spots mapped to chromosome  <dig>  using  <dig> variables  as predictor set for flesh colour the rf out-of-bag prediction explains 73 % of the phenotypic variation.
for tuber shape the expression qtls of the significant genes mapped to chromosome  <dig> to  <dig> but not to chromosome  <dig>  however, the largest number  expression qtls of the significant genes map to chromosome  <dig>  using one representative gene from each of those  <dig> chromosomes as a predictor set for tuber shape, the rf prediction explains 53 % of the variation in the trait.
no proteins or gc- or lc-variables were included since none were significantly associated to tuber shape.
only the chromosome  <dig> representative gene explains already 38 % of the phenotypic variation.
for starch gelatinization , the eqtls of the significant genes mapped to all the  <dig> chromosomes but the largest number  mapped to chromosome  <dig>  qtls for significant gc-ms peaks mapped to five chromosomes:  <dig>   <dig>   <dig>   <dig> and  <dig>  qtls for significant proteins were mapped to chromosomes  <dig> and  <dig>  these  <dig> selected variables explain 45 % of the variation in starch gelatinization.
for enzymatic discoloration the expression qtls of the significant genes mapped to four chromosomes:  <dig> , <dig> and  <dig>  qtls for significant lc-ms peaks mapped to chromosomes  <dig> and  <dig> and for significant protein spots to chromosomes  <dig> and  <dig>  these  <dig> selected variables  explain 43 % of the variation in enzymatic discoloration.
integration of ~ omics data and network visualization
with random forest we selected genes, metabolites and proteins that have a significant association with quality traits.
many of these have qtls on the same genomic regions and can be considered as redundant.
therefore we selected a single representative feature per qtl region per data set for network visualization .
the networks show partial correlations with a lasso penalty.
the nodes in the network show the phenotypic trait of interest and selected genes, metabolites and proteins .
positive partial correlations are shown in solid lines, negative partial correlations in dotted lines.fig.
<dig>  a partial correlation network of the phenotypic trait tuber flesh colour  with gene expression features , metabolites from lc-ms , metabolites from gc-ms  and proteins .
the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients.
bch = beta-carotene hydroxylase, lc_x represents metabolites derived from lc-ms with centrotype number x, gc_x represents metabolites derived from gc-ms with centrotype x, gene_x = gene with gene id x. pro_x represents a protein with protein id x.  shows the existing published part of the carotenoid pathway  <cit> , and some of the genes: bch and zep are identified by our data
fig.
<dig> a partial correlation network of tuber shape  with gene expression features .
the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients.
gene_x = gene with gene id x
fig.
<dig> a partial correlation network of dsc onset  with gene expression features , metabolites from gc-ms  and proteins .
the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients.
gc_x represents metabolites derived from gc-ms with centrotype x, gene_x = gene with gene id x. pro_x represents proteins with protein id x
fig.
<dig> a partial correlation network of enzymatic discoloration  with gene expression features , metabolites from lc-ms  and proteins .
the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients.
lc_x represents metabolites derived from lc-ms with centrotype number x, gene_x = gene with gene id x. pro_x represents proteins with protein id x
networks reconstructed from the integrated ~ omics technology platforms can provide a formal framework for investigating plant metabolism.
for example, the network shown in fig.
1a is completely based on the ~ omics data, the one in fig.
1b is based on the known carotenoid pathway as in  <cit> .
in fig.
1b we can see that bch interacts with zeaxanthin and that violoxanthin is produced.
by using only the ~ omics data  we could reproduce part of the carotenoid pathway.
this finding is extended to other traits as well in order to hypothesize how genes/proteins/metabolites are regulated.
discussion
we used random forest regression for integrating transcriptomics, metabolomics and proteomics data for prediction of four quality traits of potato: tuber flesh colour, dsc onset, tuber shape and enzymatic discoloration.
for each of these traits, we selected sets of genes, metabolites and proteins that were significant in explaining variation in the trait.
we then quantified the amount of variance explained in the prediction using these selected sets of ~ omics features, and we constructed partial correlation networks for subsets of genes, metabolites and proteins using qtl mapping information.
flesh colour
for validation of our approach, we took advantage of what was known about the relationship between the trait tuber flesh colour and previously published genes associated to flesh colour.
in our analyses, we found that beta-carotene hydroxylase  <cit>  and zeaxanthin epoxidase   <cit>  were ranked first and forty-fourth respectively for prediction of flesh colour and both compounds have previously been associated with flesh colour in potato tubers  <cit> .
in addition, for the gc-ms data we found that six metabolites were significant and out of these six, four were annotated.
it was observed that one of the compounds was identified as malic acid.
although there is no direct link between tuber flesh colour and malic acid,  <cit>  reported a correlation of malic acid to skin colour in apricot.
for the other three metabolites we did not find any connection with carotenoids or with flesh colour in the literature.
using all the significant genes across the technology platforms, the combined lc-peaks, gc-peaks and proteins, the oob variation explained  was 75 %, only slightly higher than what gene expression or lc-ms data explain by themselves  which indicates that there are correlations among these variables across the data sets.
this is indeed the case: for example, from the proteomics data, pro_ <dig> and 1684_644_ <dig> from the lc-ms data are strongly correlated with beta-carotene hydroxylase gene expression, with pearson correlation coefficients of − <dig>  and  <dig>   and, partial correlations of − <dig>  and  <dig> , respectively.
based on the correlations among the selected ~ omics data, two main clusters of ~ omics features related to flesh colour are visible .
tuber shape
for tuber shape regressed on the gene expression, lc-ms, gc-ms and proteomics data sets separately, only gene expression data was found to explain significant variation.
more than half of the eqtls of significant genes are mapped to chromosome  <dig>  tuber shape is thought to be regulated by a single locus ro on chromosome  <dig>  where round  is dominant over long   <cit> .
at the ro-locus a series of multiple alleles can explain all intermediate shapes between round  and long  <cit> .
the large number of eqtls of genes predictive for tuber shape can be due to linkage to this ro locus.
no metabolites and proteins were found as significantly predictive of tuber shape, which indicates that for this particular trait, gene expression data is more informative than metabolomics or proteomics data.
dsc onset
for dsc onset, we found  <dig> significant gene expression,  <dig> metabolite levels  and  <dig> proteins that are associated with the trait.
using those  <dig> significant variables, the variation explained was 45 %.
in order to find what genomic regions the selected genes are regulating, eqtl analyses were performed.
the analyses showed many associations with genomic regions in chromosome  <dig> with also the highest explained variation compared to other chromosomes.
for gene contig micro.
<dig> c <dig> the eqtl analysis revealed a large qtl region at the bottom of chromosome  <dig> between 73 cm and 86 cm.
eqtl analysis for est bf_lbchxxxx_0013b11_t3m.scf produced a single qtl on the same chromosomal location.
for the gc-ms data, out of five metabolites we could identify three: proline, glucopyranose and 2-piperidone that were associated to dsc onset.
the two proteins that were associated with dsc onset could not be identified.
enzymatic discoloration
transcriptomics and metabolomics analysis was previously done for enzymatic discoloration after 5 min and resulted in  <dig> significant genes and  <dig> significant lc metabolites, among which two were putatively identified as caffeoylquinic acid methyl ester and tyrosine  <cit> .
in this paper we used genetic information through qtl analysis on the one hand and prediction of the traits using rf analysis from transcriptomics, metabolomics and proteomics analysis on the other hand.
from the qtl analyses, we can identify the map position of the qtls for gene expression or metabolite signals but we expect that functional genes and genes that are only linked and that influence other pathways will show similar correlations.
for example: we obtain a large number of significant genes for prediction of flesh colour mapped to chromosome  <dig> but from previous published results  <cit> , we know that only bch on chromosome  <dig> is responsible for flesh colour, so it is likely that the remaining genes are all associated because of linkage to bch.
in rf regression, the prediction of phenotype from metabolomics, transcriptomics or proteomics data is possible in a way that genes, metabolites and/or proteins might be linked with a phenotype but independent of the genetic information  <cit> .
in this study where we combine the prediction of the phenotype with qtl information, the selected genes are not only significant for predicting the phenotype of interest but also co-localized with the phenotypic trait qtls.
further, for the predictive genes we checked their location in the genome sequence to find genomic positions of those genes irrespective of qtl analysis.
the key advantage of eqtl, mqtl or pqtl mapping in addition to the traditional mapping of phenotypic qtls, is that it connects variation at the level of rna expression, metabolite or protein abundance to variation at the level of dna.
the latter provides versatile tools for breeding whereas the first can reveal information on the biology of a trait and can direct to new candidate genes.
mapping of eqtl to the gene itself indicates that cis-regulation is responsible for the different expression levels, whereas map positions of eqtl different from the position of the corresponding genes indicate trans-regulation which allows deriving regulatory networks of genes  <cit> .
in this study, rf regression was used as a tool for data integration of metabolites, gene expression and protein profiles relating to a phenotypic trait of interest where it was used to identify leads for further exploration.
for example: flesh colour was associated with metabolites and after putative identification of the metabolic peaks  <dig> -megastigmadiene- <dig> -diol-glucoside and  <dig> -dihydroxy-4-megastigmen-9-one-glucoside were identified as carotenoid derived compounds.
such approaches give us leads for further research on the metabolites and help to hypothesize which components  are in a specific pathway of interest and the genetic basis of the genes, metabolites or proteins involved in the pathway.
metabolite peaks that are not identified but that also show an association to the trait, could in many cases be breakdown products of the carotenoid pathway  <cit> .
CONCLUSIONS
in this study, we devised and evaluated a strategy to integrate multiple ~ omics data with phenotypic traits of interest and to select sets of co-expressed genes, metabolites and proteins.
prediction of significant features was obtained through rf regression, then through a genetical genomics study  <cit>  we mapped those qtls from gene transcripts , metabolites  and proteins .
we selected a single gene transcript, metabolite and/or protein per chromosome if multiple features mapped to the same position to use for an integrated network analysis and visualization using lasso regularized partial correlations; this network can subsequently be used as a predictive network for the trait of interest.
by doing so, we selected genes, proteins and metabolites which are not only predictive for a trait of interest but also explain variation of the trait of interest.
further, we explored previously published results on the same mapping population and extended our integrative approach, previously applied to transcriptomics and metabolomics data sets  <cit>  on the combined data over all the ~ omics platforms .
for the selected combined data sets as stated before the results are interesting because there is an improvement of the r <dig> prediction value as compared to the unselected data.
this improvement is most likely due to filtering out noise variables from the data set, as is explained in  <cit>  in more detail.
for dsc onset, none of the lc-ms peaks were significant which might indicate that primary metabolism is more important for this trait compared to secondary metabolism, whereas for enzymatic discoloration the situation is the other way around.
now questions arise to which type of ~ omics data is important to study a trait.
our study indicates that, depending on the type of trait, different data sets are important.
for example, if the trait is a quality trait then metabolomics data will be useful to investigate further the biochemical pathways and to arrive at potential candidate genes   <cit> .
we selected significant genes, metabolites and proteins based on permutation tests and finally selected representative individual peaks for networks based on eqtl, mqtl or pqtl information.
network analysis was done to interpret how a particular trait is associated with gene expression, metabolite and protein data.
for these networks, we used regularized partial correlation coefficients because these quantify the correlation between two variables  while conditioning on one or several other variables  <cit> .
those genes, metabolites and proteins might be considered as leads with connections to the phenotype.
although the identified chromosomal regions do not lead automatically to genes or metabolites directly involved in the trait, and it might be necessary to know more about the metabolic pathway and indirect acting genes and/or metabolites the locations surely will help us to zoom in on potential candidate genes.
in other words, even though we are not finding the genes directly, this procedure is still helpful in guiding us towards metabolic pathways or genetic regulation of these pathways.
however, further study would be needed to analyze the combined effect of multiple qtl regions over different chromosomes.
an important limitation of this approach is validation in the absence of prior knowledge regarding genes, metabolites or proteins of the trait under investigation; a second limitation is the difficulty in identification of metabolites and especially proteins.
further, validation of such predictive networks is necessary to show that statistical association is also pointing to functional relationships in biology.
additional files
additional file 1: list of the data sets used in this article.
additional file 2: hierarchical clustering with the pearson correlation.
abbreviations
bchbeta-carotene hydroxylase
cy3cyanine 3
cy5cyanine 5
dscdifferential scanning calorimetry
eqtlexpression quantitative trait loci
gc-msgas chromatography–mass spectrometry
lc-msliquid chromatography–mass spectrometry
mqtlmetabolic quantitative trait loci
oobout-of-bag error from random forest
pqtlprotein quantitative trait loci
qtlquantitative trait loci
rfrandom forest
zepzeaxanthin epoxidase
competing interests
the authors declare that they have no competing interests.
authors’ contributions
aa and cm conceived the statistical analysis plan and arranged data for analysis.
aa performed the necessary data analyses to integrate omics data.
bk performed experiments and helped interpret the results in terms of plant biology.
aa, rgfv and cm wrote the manuscript.
all authors read and approved the final manuscript.

 in order to understand how quantitative variation in phenotypic traits is related to the underlying genetic differences between plants, and to differences in gene expression, protein constitution and metabolic variability, an approach is needed in which the combined molecular signature of the plants is shown to be predictive for the phenotypic traits of interest. we can use high-throughput ~ omics technologies, such as gene expression microarrays, mass spectrometry and protein chips to obtain molecular signatures of a population of plants. in addition we can study phenotypic differences in the same population and hypothesize that differences in the phenotypic trait in the population are related to the variation in these combined molecular profiles across the different data sets. finding the ~ omics variables that are related to the phenotypic traits of interest can then be used in two ways: 1) for prediction of the traits from these molecular profiles. for example steinfath et al., explored metabolomics data to predict agronomic phenotypes of potato crop plants grown in different environments. they identified metabolites that can be used as biomarkers for these traits and hence to improve the selection on these traits which can be implemented in breeding programs. 2) to identify functional relationships between traits and molecular networks of the plants. this means that we are not just interested in interrelating these ~ omics data sets for their own sake to find genetic and metabolic networks, but the networks and their elements should actually be predictive for a phenotypic trait of interest. on the other hand, in the context of genetics and plant breeding, we also want to be protected against finding relationships between phenotype and ~ omics data that are just caused by environmental or developmental differences. instead, we are interested to find relationships that have a basis in the genetic differences between plants. therefore a mapping population is an ideal target for this kind of study: we can study whether an observed relationship between a phenotypic trait and ~ omics variables is also based on genetic differences between the plants in the segregating population, as we can actually map the phenotypic variation as well as the variation in the ~ omics data sets. a relationship that would be just based on variation in environmental influences or conditions would not result in mapped qtls for the phenotypic traits or the ~ omics variables. we are interested in the association between a number of phenotypic traits related to tuber quality of potato and several ~ omics data sets, some of which were published previously by us. in addition, we use mapping and genotyping information since the population that we use is a mapping population. the quality traits considered are 1) potato tuber flesh colour, 2) enzymatic discoloration after peeling, 3) starch gelatinization as measured by differential scanning calorimetry and 4) tuber shape. out of these four phenotypes, starch gelatinization as measured by differential scanning calorimetry and tuber shape phenotypes were not analyzed before. the other two phenotypes were included to make the link with the previous papers and to permit validation and accuracy testing of our used network analyses. in this study, we relate transcriptomics, metabolomics and proteomics data to genetic variation in these quality traits in the mapping population. for this, we use a novel strategy consisting of three steps. first, using the same approach as in an earlier study, we apply random forest regression to find, per single trait and per individual ~ omics data set, the variables that play a significant role in the prediction of each of the quality traits. in the second step, we use qtl mapping of the quality traits and of the ~ omics variables to select variables that have a qtl cosegregating with a quality trait qtl and to remove redundancy in the set of selected variables. finally, we construct regularized partial correlation networks, for each of the quality traits, of the selected sets of variables over the four ~ omics data sets that have a genetic association with the traits per qtl. because already much is known about the regulatory genetic and metabolic pathways involved in tuber flesh colour, we used this trait to validate the approach, as we had done for individual data sets in. here we analyze all ~ omics data sets simultaneously to construct networks related to phenotypic traits consisting of features across all these data sets. we demonstrate that also for the other phenotypic traits we can find networks of small sets of interrelated gene expression profiles, proteins and metabolites that are associated with and predictive for these quality traits. methods plant material we used individuals, including the parental clones, of a diploid potato backcross population. this population is derived from an original cross between potato clones c and e and is described in detail in. all clones were grown in multi-year repeats in the field, wageningen, the netherlands during the normal potato-growing season in the netherlands. for each genotype, tubers were collected from three plants and representative samples were either used for phenotypic analysis or mechanically peeled and immediately frozen in liquid nitrogen before being ground into a fine powder and stored at −80 °c for metabolomics, transcriptomic and proteomic analyses. the determination of carotenoids was as described in. this was a targeted metabolic analysis which includes compounds like zeaxanthin, violaxanthin and a compound for which the chemical behaviour is like violaxanthin. differential scanning calorimetry is extensively used to study physical properties of starch granules. dsc-onset measurements report water temperatures at which starch granules reach their gelatinization state. dsc measurements provided gelatinization onset temperatures of starch granules for individuals of the c x e population ranging between  °c and  °c. enzymatic discoloration of tubers after peeling and being exposed to air at different time points such as after 5 min, 30 min, 3 h and difference in discoloration between 3 h and 30 min and 3 h is described in. in this study, we considered only enzymatic discoloration after 5 min since the measurements at later time points are all highly correlated. potato tuber shape was scored between and. omics data sets we generated different types of ~ omics data to integrate and find novel relationships among them. parts of the ~ omics data sets were already published individually, for example, the transcriptomics data, gc-ms data, lc-ms data. only the proteomics data was generated more recently and these data have not been considered for a combined analysis over all ~ omics technologies simultaneously. in this section, we briefly summarize the different omics data generation methods and refer to for further details. the phenotypic data and ~ omics data are provided in a supplementary file. transcriptomic data rna was extracted from the samples using the hot phenol method described previously. according to the manufacturer’s protocol starting with 2 μg of purified total rna ). for additional data analyses only genes with a pearson correlation coefficient higher than between the cy and cy datasets were included resulting in expressed genes. we took into account only the cy gene expression signals for further statistical analysis. for visualization, we used the gene nomenclature in the following way: gene_gene id. the number refers to the gene id of the supplementary material of. lc-ms, data generation, processing and identification potato tuber samples were analyzed for variation in semi-polar metabolite composition using an untargeted accurate mass lc-ms approach. in total mass signals were obtained from mass spectrometry and selected based on the following criteria: signals should be present in at least samples and with at least one amplitude higher than units. finally, data redundancy of signals derived from the same compound, i.e. isotopes, adducts and in-source fragments, was removed by by retention time-dependent clustering. this clustering of the signals resulted in reconstructed metabolites. the untargeted metabolites are represented as centrotype_mass_scan number as in. as an example: 818_795_ means that the centrotype number is the mass number and the scan number for visualization and simplicity we used lc-centrotype number. for more information on data generation, processing and identification see the methods of. gc-ms data generation, processing and identification gc-ms data were generated from the same genotypes of the cxe population. detailed materials and methods for gc-ms data generation, processing and identification are as in. in short, raw data were processed by chromatof software and masslynx software, and further analysis was performed using metalign software to extract and align the mass signals. mass signals that were present in fewer than two samples were discarded. signal redundancy per metabolite was removed by means of clustering, and mass spectra were reconstructed. this resulted in reconstructed polar metabolites. compounds were subjected to tentative identification by matching to the nist and wiley spectral libraries. library hits were manually curated, and a series of commercial standards was used to check annotation. for visualization and simplicity we used gc-centrotype number. for example, 7275_20721700_ means that the centrotype number is the mass number and the scan number for visualization only gc_ is used. proteomics data generation, processing and identification protein extraction total protein was extracted from approximately g of ground tuber material, to which 1 ml of pre-heated lysis buffer, sds, dtt, pvp-p ) was added. samples were homogenized and protein amount was measured using the rc/dc assay. protein labelling a single lysine per protein molecule was labelled using the fluorescent cydyes from the difference gel electrophoresis technology according to the manufacturer’s protocol. the internal standard was labelled with cy and consists of an equal mixture of randomly chosen samples of the experiment. every 2d-gel contains one sample labelled with cy one labelled with cy and the internal standard labelled with cy this means that every sample on each gel can be compared by using the internal standard sample labelled with cy 2d-electrophoresis electrophoresis was performed using first dimension electrophoresis on 24 cm immobilized ph gradient strips with a linear ph range from to on an ettan ipgphor isoelectric focusing system. cydye labelled samples were loaded to the strips and focusing was run for 18 h at 20 °c with the following settings: 3 h 150 v, 3 h 300 v, from 300 v to 1000 v in 6 h, from 1000 v to  v in 1 h and finally 5 h at  v. after equilibration in the dark at room temperature in equilibration buffer, sds 2 % ) containing dtt 1 % the second dimension electrophoresis was run on the ettan dalt system on precast  % sds polyacrylamide slab gel and buffers from ge healthcare/amersham biosciences. the separated cydye-labelled proteins were visualized by scanning with a ettan dige imager, using for cy an 480 nm laser and an emission filter of 530 nm, for cy an 540 nm laser and an emission filter of 595 nm and for cy an 635 nm laser and an emission filter of 680 nm. gel images were analysed with the decyder software version according to decyder 2dv manual; ge healthcare/amersham biosciences. detected spots were filtered based on spot volume larger than relative value to exclude spots that could be just background noise or dust particles. the internal standard in each gel was used to automatically match all images to the reference value was used. each of the proteins is presented by pro_x where “x” represents consecutive protein numbers. protein identification spots of interest were excised from gel using the ettan spot picker. peptide mass determinations were carried out using the applied biosystems proteomics analyzer. both pmf and ms/ms in reflectron mode analyses were carried out with the samples. calibration was carried out with a peptide mass calibration kit. proteins were identified by searching against the ncbi ‘viridiplantae’ and an est ‘viridiplantae-eudicots’ database using mascot. random forest regression random forest was used for regression of the phenotypic trait such as flesh colour, tuber shape, dsc onset, enzymatic discoloration on the transcriptomics, gc-tof-ms signals, lc-ms signals, 2d-dige proteomics. rf constructs a predictive model for the response using all predictors but quantifies the importance of each, here the metabolites, gene expression and proteins in explaining the variation present in the trait. rf by itself does not provide significance levels of individual metabolites and does not perform a variable selection to choose a predictive subset of associated metabolites. therefore, we included a permutation test to indicate significance of the association of a metabolite with a trait. in each of permutations of the trait values we estimated the variance explained by the rf model and the variable importance of each metabolite in terms of the decrease in node impurities. we ordered node purity values from the permuted data sets and took the percentile from the distribution of impurity values as the significance threshold of the individual metabolites. the same procedure was done for r values of the model: the percentile was taken as a significance threshold for the rf model. r in random forest is not just a measure of goodness-of-fit of the data at hand but is determined on left-out samples so it should be interpreted as a measure for predictive quality of the random forest on independent samples that have the same properties as the in-bag sample. random forest analysis was done using r statistical software using the randomforest package. qtl analysis and cis- and trans-eqtls we mapped expression qtls from the gene expression data, metabolite qtls both for lc-ms and gc-ms, and protein qtls from the proteomics data to find regions on the genome explaining genetic variation in gene expression, metabolite and protein values using the integrated linkage map of the c and e parents for qtl analysis. further, we used the potato genome physical map to investigate eqtl physical positions of genes identified as predictive for phenotypic traits from the random forest analyses. the potato oligo microarray used in the experiments contains features based on a potato unigene set. to allow discrimination between cis- and trans-eqtls all unigenes were blasted against the genome scaffold sequences, predicted coding sequences and predicted gene regions. features with a unique and significant hit were assigned to genome scaffolds for which the majority has chromosome information. identified qtls on the same linkage group as their physical map position are identified as cis-acting while qtls on different linkage groups are defined as trans-acting. features on the array for which no physical map position could be assigned are classified as unknown. all qtl analyses were done using metanetwork package in r. network reconstruction from the rf analyses, for each of the traits, we obtain a list of significant genes, metabolites and proteins. for visualization in an association network, we reduced the number of genes by checking chromosomal positions of each of the predicted significant genes, metabolites and proteins. we took the most significant single gene, lc-ms signal, gc-ms signal and/or protein per chromosome as representative for that ~ omics data set and chromosome and made a full order regularized partial correlation network with these genes, metabolites, proteins and traits as nodes in the network and the strength of the interaction as edges in the network. partial correlation measures the correlation between two variables after their linear dependence on other variables is removed. it can distinguish between direct and indirect associations whereas correlation-based network cannot and often yield many spurious edges. a detailed description of regularized partial correlation is provided in.  selection of ~ omics features predictive for quality traits the prediction of variation in potato tuber flesh colour from the gene expression, lc-ms, and protein data sets was quite high, but much lower for the gc-ms data. for flesh colour, the microarray data and the lc-ms data were equally good for prediction in terms of the explained variance, but the numbers of significant features were very different: the prediction of flesh colour using significant lc-ms features is almost the same as for significant gene expression profiles from the microarray data set. from the gene expression data, the gene which ranks first and third with respect to variable importance for predicting flesh colour was a beta-carotene hydroxylase. two oligos were present on the array targeting the same beta-carotene hydroxylase gene, hence the two high ranks for the same gene. another gene from the carotenoid pathway, zeaxanthin epoxidase ranked forty-fourth. based on our current knowledge of potato tuber flesh colour and carotenoid content, these two genes were expected to be associated with flesh colour. from the gc-ms data, out of six significant metabolites, four had an annotation: malic acid, 2-4-5-trihydroxypentanoic acid, glucopyranose, 2-butenedioic acid- and bis ester.table percentage variance explained in out-of-bag prediction by random forest models using all genes, lc-peaks, gc-peaks or proteins separately non-significant models are indicated as ns gene + lc + gc + protein; gene + gc + protein; gene + lc + protein; the numbers of significant genes, lc-peaks, gc-peaks or proteins are between brackets. non-significant models are indicated as ns for the three other traits the gene expression data was more predictive than the metabolomics and proteomics data sets, and for tuber shape actually no significant features were found for gc-ms, lc-ms and the protein data, while significant gene expression differences explain 55 % of the variation in tuber shape. for starch gelatinization no lc-ms features were significant, for enzymatic discoloration no gc-ms features were significant. associated genomic regions using random forest regression, we generated lists of candidate genes, metabolite features and proteins that are predictive for quality traits of potato tubers. however, from this prediction we cannot conclude that these associations necessarily have a basis in genetic differences between the plants since they could also be caused by environmental variation in both the trait and the levels of the features, or in developmental differences. therefore we also investigated qtl positions for both the phenotypic traits and the ~ omics features. when co-segregation of an eqtl, mqtl or pqtl with a trait qtl is observed, this could imply a functional relationship or identify causal genes/proteins/metabolites. however, this is not necessarily the case since any linked but functionally unrelated qtls will also show this co-segregation. still, these sets of predictive genes highlight genomic regions of interest, comparable to standard qtl analysis, and can provide additional information or help to narrow down the region of the causative polymorphism. for tuber flesh colour the eqtls of significantly predictive gene expression profiles from rf mapped to eight chromosomes: and a large number of these eqtls were mapped to chromosomes and significant gc-ms peaks mapped to two chromosomes, and qtls for significant lc-ms peaks mapped to chromosomes and and qtls for significant protein spots mapped to chromosome using variables as predictor set for flesh colour the rf out-of-bag prediction explains 73 % of the phenotypic variation. for tuber shape the expression qtls of the significant genes mapped to chromosome to but not to chromosome however, the largest number expression qtls of the significant genes map to chromosome using one representative gene from each of those chromosomes as a predictor set for tuber shape, the rf prediction explains 53 % of the variation in the trait. no proteins or gc- or lc-variables were included since none were significantly associated to tuber shape. only the chromosome representative gene explains already 38 % of the phenotypic variation. for starch gelatinization, the eqtls of the significant genes mapped to all the chromosomes but the largest number mapped to chromosome qtls for significant gc-ms peaks mapped to five chromosomes: and qtls for significant proteins were mapped to chromosomes and these selected variables explain 45 % of the variation in starch gelatinization. for enzymatic discoloration the expression qtls of the significant genes mapped to four chromosomes:, and qtls for significant lc-ms peaks mapped to chromosomes and and for significant protein spots to chromosomes and these selected variables explain 43 % of the variation in enzymatic discoloration. integration of ~ omics data and network visualization with random forest we selected genes, metabolites and proteins that have a significant association with quality traits. many of these have qtls on the same genomic regions and can be considered as redundant. therefore we selected a single representative feature per qtl region per data set for network visualization. the networks show partial correlations with a lasso penalty. the nodes in the network show the phenotypic trait of interest and selected genes, metabolites and proteins. positive partial correlations are shown in solid lines, negative partial correlations in dotted lines.fig. a partial correlation network of the phenotypic trait tuber flesh colour with gene expression features, metabolites from lc-ms, metabolites from gc-ms and proteins. the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients. bch = beta-carotene hydroxylase, lc_x represents metabolites derived from lc-ms with centrotype number x, gc_x represents metabolites derived from gc-ms with centrotype x, gene_x = gene with gene id x. pro_x represents a protein with protein id x. shows the existing published part of the carotenoid pathway, and some of the genes: bch and zep are identified by our data fig. a partial correlation network of tuber shape with gene expression features. the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients. gene_x = gene with gene id x fig. a partial correlation network of dsc onset with gene expression features, metabolites from gc-ms and proteins. the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients. gc_x represents metabolites derived from gc-ms with centrotype x, gene_x = gene with gene id x. pro_x represents proteins with protein id x fig. a partial correlation network of enzymatic discoloration with gene expression features, metabolites from lc-ms and proteins. the dotted lines represent negative partial correlation coefficients, solid lines represent positive partial correlation coefficients. lc_x represents metabolites derived from lc-ms with centrotype number x, gene_x = gene with gene id x. pro_x represents proteins with protein id x networks reconstructed from the integrated ~ omics technology platforms can provide a formal framework for investigating plant metabolism. for example, the network shown in fig. 1a is completely based on the ~ omics data, the one in fig. 1b is based on the known carotenoid pathway as in. in fig. 1b we can see that bch interacts with zeaxanthin and that violoxanthin is produced. by using only the ~ omics data we could reproduce part of the carotenoid pathway. this finding is extended to other traits as well in order to hypothesize how genes/proteins/metabolites are regulated. discussion we used random forest regression for integrating transcriptomics, metabolomics and proteomics data for prediction of four quality traits of potato: tuber flesh colour, dsc onset, tuber shape and enzymatic discoloration. for each of these traits, we selected sets of genes, metabolites and proteins that were significant in explaining variation in the trait. we then quantified the amount of variance explained in the prediction using these selected sets of ~ omics features, and we constructed partial correlation networks for subsets of genes, metabolites and proteins using qtl mapping information. flesh colour for validation of our approach, we took advantage of what was known about the relationship between the trait tuber flesh colour and previously published genes associated to flesh colour. in our analyses, we found that beta-carotene hydroxylase and zeaxanthin epoxidase were ranked first and forty-fourth respectively for prediction of flesh colour and both compounds have previously been associated with flesh colour in potato tubers. in addition, for the gc-ms data we found that six metabolites were significant and out of these six, four were annotated. it was observed that one of the compounds was identified as malic acid. although there is no direct link between tuber flesh colour and malic acid, reported a correlation of malic acid to skin colour in apricot. for the other three metabolites we did not find any connection with carotenoids or with flesh colour in the literature. using all the significant genes across the technology platforms, the combined lc-peaks, gc-peaks and proteins, the oob variation explained was 75 %, only slightly higher than what gene expression or lc-ms data explain by themselves which indicates that there are correlations among these variables across the data sets. this is indeed the case: for example, from the proteomics data, pro_ and 1684_644_ from the lc-ms data are strongly correlated with beta-carotene hydroxylase gene expression, with pearson correlation coefficients of − and and, partial correlations of − and, respectively. based on the correlations among the selected ~ omics data, two main clusters of ~ omics features related to flesh colour are visible. tuber shape for tuber shape regressed on the gene expression, lc-ms, gc-ms and proteomics data sets separately, only gene expression data was found to explain significant variation. more than half of the eqtls of significant genes are mapped to chromosome tuber shape is thought to be regulated by a single locus ro on chromosome where round is dominant over long. at the ro-locus a series of multiple alleles can explain all intermediate shapes between round and long. the large number of eqtls of genes predictive for tuber shape can be due to linkage to this ro locus. no metabolites and proteins were found as significantly predictive of tuber shape, which indicates that for this particular trait, gene expression data is more informative than metabolomics or proteomics data. dsc onset for dsc onset, we found significant gene expression, metabolite levels and proteins that are associated with the trait. using those significant variables, the variation explained was 45 %. in order to find what genomic regions the selected genes are regulating, eqtl analyses were performed. the analyses showed many associations with genomic regions in chromosome with also the highest explained variation compared to other chromosomes. for gene contig micro. c the eqtl analysis revealed a large qtl region at the bottom of chromosome between 73 cm and 86 cm. eqtl analysis for est bf_lbchxxxx_0013b11_t3m.scf produced a single qtl on the same chromosomal location. for the gc-ms data, out of five metabolites we could identify three: proline, glucopyranose and 2-piperidone that were associated to dsc onset. the two proteins that were associated with dsc onset could not be identified. enzymatic discoloration transcriptomics and metabolomics analysis was previously done for enzymatic discoloration after 5 min and resulted in significant genes and significant lc metabolites, among which two were putatively identified as caffeoylquinic acid methyl ester and tyrosine. in this paper we used genetic information through qtl analysis on the one hand and prediction of the traits using rf analysis from transcriptomics, metabolomics and proteomics analysis on the other hand. from the qtl analyses, we can identify the map position of the qtls for gene expression or metabolite signals but we expect that functional genes and genes that are only linked and that influence other pathways will show similar correlations. for example: we obtain a large number of significant genes for prediction of flesh colour mapped to chromosome but from previous published results, we know that only bch on chromosome is responsible for flesh colour, so it is likely that the remaining genes are all associated because of linkage to bch. in rf regression, the prediction of phenotype from metabolomics, transcriptomics or proteomics data is possible in a way that genes, metabolites and/or proteins might be linked with a phenotype but independent of the genetic information. in this study where we combine the prediction of the phenotype with qtl information, the selected genes are not only significant for predicting the phenotype of interest but also co-localized with the phenotypic trait qtls. further, for the predictive genes we checked their location in the genome sequence to find genomic positions of those genes irrespective of qtl analysis. the key advantage of eqtl, mqtl or pqtl mapping in addition to the traditional mapping of phenotypic qtls, is that it connects variation at the level of rna expression, metabolite or protein abundance to variation at the level of dna. the latter provides versatile tools for breeding whereas the first can reveal information on the biology of a trait and can direct to new candidate genes. mapping of eqtl to the gene itself indicates that cis-regulation is responsible for the different expression levels, whereas map positions of eqtl different from the position of the corresponding genes indicate trans-regulation which allows deriving regulatory networks of genes. in this study, rf regression was used as a tool for data integration of metabolites, gene expression and protein profiles relating to a phenotypic trait of interest where it was used to identify leads for further exploration. for example: flesh colour was associated with metabolites and after putative identification of the metabolic peaks -megastigmadiene- -diol-glucoside and -dihydroxy-4-megastigmen-9-one-glucoside were identified as carotenoid derived compounds. such approaches give us leads for further research on the metabolites and help to hypothesize which components are in a specific pathway of interest and the genetic basis of the genes, metabolites or proteins involved in the pathway. metabolite peaks that are not identified but that also show an association to the trait, could in many cases be breakdown products of the carotenoid pathway.  in this study, we devised and evaluated a strategy to integrate multiple ~ omics data with phenotypic traits of interest and to select sets of co-expressed genes, metabolites and proteins. prediction of significant features was obtained through rf regression, then through a genetical genomics study we mapped those qtls from gene transcripts, metabolites and proteins. we selected a single gene transcript, metabolite and/or protein per chromosome if multiple features mapped to the same position to use for an integrated network analysis and visualization using lasso regularized partial correlations; this network can subsequently be used as a predictive network for the trait of interest. by doing so, we selected genes, proteins and metabolites which are not only predictive for a trait of interest but also explain variation of the trait of interest. further, we explored previously published results on the same mapping population and extended our integrative approach, previously applied to transcriptomics and metabolomics data sets on the combined data over all the ~ omics platforms. for the selected combined data sets as stated before the results are interesting because there is an improvement of the r prediction value as compared to the unselected data. this improvement is most likely due to filtering out noise variables from the data set, as is explained in in more detail. for dsc onset, none of the lc-ms peaks were significant which might indicate that primary metabolism is more important for this trait compared to secondary metabolism, whereas for enzymatic discoloration the situation is the other way around. now questions arise to which type of ~ omics data is important to study a trait. our study indicates that, depending on the type of trait, different data sets are important. for example, if the trait is a quality trait then metabolomics data will be useful to investigate further the biochemical pathways and to arrive at potential candidate genes. we selected significant genes, metabolites and proteins based on permutation tests and finally selected representative individual peaks for networks based on eqtl, mqtl or pqtl information. network analysis was done to interpret how a particular trait is associated with gene expression, metabolite and protein data. for these networks, we used regularized partial correlation coefficients because these quantify the correlation between two variables while conditioning on one or several other variables. those genes, metabolites and proteins might be considered as leads with connections to the phenotype. although the identified chromosomal regions do not lead automatically to genes or metabolites directly involved in the trait, and it might be necessary to know more about the metabolic pathway and indirect acting genes and/or metabolites the locations surely will help us to zoom in on potential candidate genes. in other words, even though we are not finding the genes directly, this procedure is still helpful in guiding us towards metabolic pathways or genetic regulation of these pathways. however, further study would be needed to analyze the combined effect of multiple qtl regions over different chromosomes. an important limitation of this approach is validation in the absence of prior knowledge regarding genes, metabolites or proteins of the trait under investigation; a second limitation is the difficulty in identification of metabolites and especially proteins. further, validation of such predictive networks is necessary to show that statistical association is also pointing to functional relationships in biology. additional files additional file 1: list of the data sets used in this article. additional file 2: hierarchical clustering with the pearson correlation. abbreviations bchbeta-carotene hydroxylase cy3cyanine 3 cy5cyanine 5 dscdifferential scanning calorimetry eqtlexpression quantitative trait loci gc-msgas chromatography–mass spectrometry lc-msliquid chromatography–mass spectrometry mqtlmetabolic quantitative trait loci oobout-of-bag error from random forest pqtlprotein quantitative trait loci qtlquantitative trait loci rfrandom forest zepzeaxanthin epoxidase competing interests the authors declare that they have no competing interests. authors’ contributions aa and cm conceived the statistical analysis plan and arranged data for analysis. aa performed the necessary data analyses to integrate omics data. bk performed experiments and helped interpret the results in terms of plant biology. aa, rgfv and cm wrote the manuscript. all authors read and approved the final manuscript. 
in this study, we made a strategy for selecting and integrating multiple ~ omics data using random forest method and selected representative individual peaks for networks based on eqtl, mqtl or pqtl information. we used random forest regression for integrating multiple ~ omics data for prediction of four quality traits of potato: tuber flesh colour, dsc onset, tuber shape and enzymatic discoloration. for tuber shape regressed on the gene expression, lc-ms, gc-ms and proteomics data sets separately, only gene expression data was found to explain significant variation. this gives us networks of relatively small sets of interrelated ~ omics variables that can predict, with higher accuracy, a quality trait of interest. transcriptomics and metabolomics analysis on enzymatic discoloration after 5 min resulted in significant genes and significant lc metabolites, among which two were putatively identified as caffeoylquinic acid methyl ester and tyrosine. combining all the significant genes, lc-peaks, gc-peaks and proteins, the variation explained was 75 %, only slightly more than what gene expression or lc-ms data explain by themselves which indicates that there are correlations among the variables across data sets. in order to find genetic and metabolic pathways related to phenotypic traits of interest, we analyzed gene expression data, metabolite data obtained with gc-ms and lc-ms, proteomics data and a selected set of tuber quality phenotypic data from a diploid segregating mapping population of potato. in this study we present an approach to integrate these ~ omics data sets for the purpose of predicting phenotypic traits. for tuber flesh colour beta-carotene hydroxylase and zeaxanthin epoxidase were ranked first and forty-fourth respectively both of which have previously been associated with flesh colour in potato tubers. expression qtl analyses showed many associations with genomic regions in chromosome with also the highest explained variation compared to other chromosomes. network analysis was done to interpret how a particular trait is associated with gene expression, metabolite and protein data. using those significant variables, the variation explained was 45 %. for dsc onset, we found significant gene expression, metabolite levels and proteins that are associated with the trait.
 in order to find genetic and metabolic pathways related to phenotypic traits of interest, we analyzed gene expression data, metabolite data obtained with gc-ms and lc-ms, proteomics data and a selected set of tuber quality phenotypic data from a diploid segregating mapping population of potato. in this study we present an approach to integrate these ~ omics data sets for the purpose of predicting phenotypic traits. this gives us networks of relatively small sets of interrelated ~ omics variables that can predict, with higher accuracy, a quality trait of interest.  we used random forest regression for integrating multiple ~ omics data for prediction of four quality traits of potato: tuber flesh colour, dsc onset, tuber shape and enzymatic discoloration. for tuber flesh colour beta-carotene hydroxylase and zeaxanthin epoxidase were ranked first and forty-fourth respectively both of which have previously been associated with flesh colour in potato tubers. combining all the significant genes, lc-peaks, gc-peaks and proteins, the variation explained was 75 %, only slightly more than what gene expression or lc-ms data explain by themselves which indicates that there are correlations among the variables across data sets. for tuber shape regressed on the gene expression, lc-ms, gc-ms and proteomics data sets separately, only gene expression data was found to explain significant variation. for dsc onset, we found significant gene expression, metabolite levels and proteins that are associated with the trait. using those significant variables, the variation explained was 45 %. expression qtl analyses showed many associations with genomic regions in chromosome with also the highest explained variation compared to other chromosomes. transcriptomics and metabolomics analysis on enzymatic discoloration after 5 min resulted in significant genes and significant lc metabolites, among which two were putatively identified as caffeoylquinic acid methyl ester and tyrosine.  in this study, we made a strategy for selecting and integrating multiple ~ omics data using random forest method and selected representative individual peaks for networks based on eqtl, mqtl or pqtl information. network analysis was done to interpret how a particular trait is associated with gene expression, metabolite and protein data. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords data integrationgenetical genomicsnetworksrandom foreststatistical methods for omics data integration and analysis smodia 2014heraklion, crete, greece 10- november http://smodia com/issue-copyright-statement© the author 2016 