A Critical Assessment of Feature Selection Methods for Biomarker Discovery in Clinical Proteomics

In this paper, we compare the performance of six different feature selection methods for LC-MS-based proteomics and metabolomics biomarker discovery—t test, the Mann–Whitney–Wilcoxon test (mww test), nearest shrunken centroid (NSC), linear support vector machine–recursive features elimination (SVM-RFE), principal component discriminant analysis (PCDA), and partial least squares discriminant analysis (PLSDA)—using human urine and porcine cerebrospinal fluid samples that were spiked with a range of peptides at different concentration levels. The ideal feature selection method should select the complete list of discriminating features that are related to the spiked peptides without selecting unrelated features. Whereas many studies have to rely on classification error to judge the reliability of the selected biomarker candidates, we assessed the accuracy of selection directly from the list of spiked peptides. The feature selection methods were applied to data sets with different sample sizes and extents of sample class separation determined by the concentration level of spiked compounds. For each feature selection method and data set, the performance for selecting a set of features related to spiked compounds was assessed using the harmonic mean of the recall and the precision (f-score) and the geometric mean of the recall and the true negative rate (g-score). We conclude that the univariate t test and the mww test with multiple testing corrections are not applicable to data sets with small sample sizes (n = 6), but their performance improves markedly with increasing sample size up to a point (n > 12) at which they outperform the other methods. PCDA and PLSDA select small feature sets with high precision but miss many true positive features related to the spiked peptides. NSC strikes a reasonable compromise between recall and precision for all data sets independent of spiking level and number of samples. Linear SVM-RFE performs poorly for selecting features related to the spiked compounds, even though the classification error is relatively low.Biomarkers play an important role in advancing medical research through the early diagnosis of disease and prognosis of treatment interventions (1, 2). Biomarkers may be proteins, peptides, or metabolites, as well as mRNAs or other kinds of nucleic acids (e.g. microRNAs) whose levels change in relation to the stage of a given disease and which may be used to accurately assign the disease stage of a patient. The accurate selection of biomarker candidates is crucial, because it determines the outcome of further validation studies and the ultimate success of efforts to develop diagnostic and prognostic assays with high specificity and sensitivity. The success of biomarker discovery depends on several factors: consistent and reproducible phenotyping of the individuals from whom biological samples are obtained; the quality of the analytical methodology, which in turn determines the quality of the collected data; the accuracy of the computational methods used to extract quantitative and molecular identity information to define the biomarker candidates from raw analytical data; and finally the performance of the applied statistical methods in the selection of a limited list of compounds with the potential to discriminate between predefined classes of samples. De novo biomarker research consists of a biomarker discovery part and a biomarker validation part (3). Biomarker discovery uses analytical techniques that try to measure as many compounds as possible in a relatively low number of samples. The goal of subsequent data preprocessing and statistical analysis is to select a limited number of candidates, which are subsequently subjected to targeted analyses in large number of samples for validation.Advanced technology, such as high-performance liquid chromatography–mass spectrometry (LC-MS),¹ is increasingly applied in biomarker discovery research. Such analyses detect tens of thousands of compounds, as well as background-related signals, in a single biological sample, generating enormous amounts of multivariate data. Data preprocessing workflows reduce data complexity considerably by trying to extract only the information related to compounds resulting in a quantitative feature matrix, in which rows and columns correspond to samples and extracted features, respectively, or vice versa. Features may also be related to data preprocessing artifacts, and the ratio of such erroneous features to compound-related features depends on the performance of the data preprocessing workflow (4). Preprocessed LC-MS data sets contain a large number of features relative to the sample size. These features are characterized by their m/z value and retention time, and in the ideal case they can be combined and linked to compound identities such as metabolites, peptides, and proteins. In LC-MS-based proteomics and metabolomics studies, sample analysis is so time consuming that it is practically impossible to increase the number of samples to a level that balances the number of features in a data set. Therefore, the success of biomarker discovery depends on powerful feature selection methods that can deal with a low sample size and a high number of features. Because of the unfavorable statistical situation and the risk of overfitting the data, it is ultimately pivotal to validate the selected biomarker candidates in a larger set of independent samples, preferably in a double-blinded fashion, using targeted analytical methods (1).Biomarker selection is often based on classification methods that are preceded by feature selection methods (filters) or which have built-in feature selection modules (wrappers and embedded methods) that can be used to select a list of compounds/peaks/features that provide the best classification performance for predefined sample groups (e.g. healthy versus diseased) (5). Classification methods are able to classify an unknown sample into a predefined sample class. Univariate feature selection methods such as filters (t test or Wilcoxon–Mann–Whitney tests) cannot be used for sample classification. Other classification methods such as the nearest shrunken centroid method have intrinsic feature selection ability, whereas other classification methods such as principal component discriminant analysis (PCDA) and partial least squares regression coupled with discriminant analysis (PLSDA) should be augmented with a feature selection method. There are classifiers having no feature selection option that perform the classification using all variables, such as support vector machines that use non-linear kernels (6). Classification methods without the ability to select features cannot be used for biomarker discovery, because these methods aim to classify samples into predefined classes but cannot identify the limited number of variables (features or compounds) that form the basis of the classification (6, 7). Different statistical methods with feature selection have been developed according to the complexity of the analyzed data, and these have been extensively reviewed (5, 6, 8, 9). Ways of optimizing such methods to improve sensitivity and specificity are a major topic in current biomarker discovery research and in the many “omics-related” research areas (6, 10, 11). Comparisons of classification methods with respect to their classification and learning performance have been initiated. Van der Walt et al. (12) focused on finding the most accurate classifiers for simulated data sets with sample sizes ranging from 20 to 100. Rubingh et al. (13) compared the influence of sample size in an LC-MS metabolomics data set on the performance of three different statistical validation tools: cross validation, jack-knifing model parameters, and a permutation test. That study concluded that for small sample sets, the outcome of these validation methods is influenced strongly by individual samples and therefore cannot be trusted, and the validation tool cannot be used to indicate problems due to sample size or the representativeness of sampling. This implies that reducing the dimensionality of the feature space is critical when approaching a classification problem in which the number of features exceeds the number of samples by a large margin. Dimensionality reduction retains a smaller set of features to bring the feature space in line with the sample size and thus allow the application of classification methods that perform with acceptable accuracy only when the sample size and the feature size are similar.In this study we compared different classification methods focusing on feature selection in two types of spiked LC-MS data sets that mimic the situation of a biomarker discovery study. Our results provide guidelines for researchers who will engage in biomarker discovery or other differential profiling “omics” studies with respect to sample size and selecting the most appropriate feature selection method for a given data set. We evaluated the following approaches: univariate t test and Mann–Whitney–Wilcoxon test (mww test) with multiple testing correction (14), nearest shrunken centroid (NSC) (15, 16), support vector machine–recursive features elimination (SVM-RFE) (17), PLSDA (18), and PCDA (19). PCDA and PLSDA were combined with the rank-product as a feature selection criterion (20). These methods were evaluated with data sets having three characteristics: different biological background, varying sample size, and varying within- and between-class variability of the added compounds. Data were acquired via LC-MS from human urine and porcine cerebrospinal fluid (CSF) samples that were spiked with a set of known peptides (true positives) at different concentration levels. These samples were then combined in two classes containing peptides spiked at low and high concentration levels. The performance of the classification methods with feature selection was measured based on their ability to select features that were related to the spiked peptides. Because true positives were known in our data set, we compared performance based on the f-score (the harmonic mean of precision and recall) and the g-score (the geometric mean of accuracy).     

Genomic prediction based on data from three layer lines: a comparison between linear methods

Background

The prediction accuracy of several linear genomic prediction models, which have previously been used for within-line genomic prediction, was evaluated for multi-line genomic prediction.

Methods

Compared to a conventional BLUP (best linear unbiased prediction) model using pedigree data, we evaluated the following genomic prediction models: genome-enabled BLUP (GBLUP), ridge regression BLUP (RRBLUP), principal component analysis followed by ridge regression (RRPCA), BayesC and Bayesian stochastic search variable selection. Prediction accuracy was measured as the correlation between predicted breeding values and observed phenotypes divided by the square root of the heritability. The data used concerned laying hens with phenotypes for number of eggs in the first production period and known genotypes. The hens were from two closely-related brown layer lines (B1 and B2), and a third distantly-related white layer line (W1). Lines had 1004 to 1023 training animals and 238 to 240 validation animals. Training datasets consisted of animals of either single lines, or a combination of two or all three lines, and had 30 508 to 45 974 segregating single nucleotide polymorphisms.

Results

Genomic prediction models yielded 0.13 to 0.16 higher accuracies than pedigree-based BLUP. When excluding the line itself from the training dataset, genomic predictions were generally inaccurate. Use of multiple lines marginally improved prediction accuracy for B2 but did not affect or slightly decreased prediction accuracy for B1 and W1. Differences between models were generally small except for RRPCA which gave considerably higher accuracies for B2. Correlations between genomic predictions from different methods were higher than 0.96 for W1 and higher than 0.88 for B1 and B2. The greater differences between methods for B1 and B2 were probably due to the lower accuracy of predictions for B1 (~0.45) and B2 (~0.40) compared to W1 (~0.76).

Conclusions

Multi-line genomic prediction did not affect or slightly improved prediction accuracy for closely-related lines. For distantly-related lines, multi-line genomic prediction yielded similar or slightly lower accuracies than single-line genomic prediction. Bayesian variable selection and GBLUP generally gave similar accuracies. Overall, RRPCA yielded the greatest accuracies for two lines, suggesting that using PCA helps to alleviate the “n ≪ p” problem in genomic prediction.

Electronic supplementary material

The online version of this article (doi:10.1186/s12711-014-0057-5) contains supplementary material, which is available to authorized users.     

A high-throughput processing service for retention time alignment of complex proteomics and metabolomics LC-MS data

Warp2D is a novel time alignment approach, which uses the overlapping peak volume of the reference and sample peak lists to correct misleading peak shifts. Here, we present an easy-to-use web interface for high-throughput Warp2D batch processing time alignment service using the Dutch Life Science Grid, reducing processing time from days to hours. This service provides the warping function, the sample chromatogram peak list with adjusted retention times and normalized quality scores based on the sum of overlapping peak volume of all peaks. Heat maps before and after time alignment are created from the arithmetic mean of the sum of overlapping peak area rearranged with hierarchical clustering, allowing the quality control of the time alignment procedure. Taverna workflow and command line tool are provided for remote processing of local user data. AVAILABILITY: online data processing service is available at http://www.nbpp.nl/warp2d.html. Taverna workflow is available at myExperiment with title '2D Time Alignment-Webservice and Workflow' at http://www.myexperiment.org/workflows/1283.html. Command line tool is available at http://www.nbpp.nl/Warp2D_commandline.zip. CONTACT: p.l.horvatovich@rug.nl SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.     

MicroRNA transcriptome profiles during swine skeletal muscle development

Background

Dietary polyunsaturated fatty acids (PUFA), in particular the long chain marine fatty acids docosahexaenoic (DHA) and eicosapentaenoic (EPA), are linked to many health benefits in humans and in animal models. Little is known of the molecular response to DHA and EPA of the small intestine, and the potential contribution of this organ to the beneficial effects of these fatty acids. Here, we assessed gene expression changes induced by DHA and EPA in the wildtype C57BL/6J murine small intestine using whole genome microarrays and functionally characterized the most prominent biological process.

Results

The main biological process affected based on gene expression analysis was lipid metabolism. Fatty acid uptake, peroxisomal and mitochondrial beta-oxidation, and omega-oxidation of fatty acids were all increased. Quantitative real time PCR, and -in a second animal experiment- intestinal fatty acid oxidation measurements confirmed significant gene expression differences and showed in a dose-dependent manner significant changes at biological functional level. Furthermore, no major changes in the expression of lipid metabolism genes were observed in the colon.

Conclusion

We show that marine n-3 fatty acids regulate small intestinal gene expression and increase fatty acid oxidation. Since this organ contributes significantly to whole organism energy use, this effect on the small intestine may well contribute to the beneficial physiological effects of marine PUFAs under conditions that will normally lead to development of obesity, insulin resistance and diabetes.     

Ultrastructure of the sodium pump: Comparison of thin sectioning, negative staining, and freeze-fracture of purified, membrane-bound (Na+, K+)-ATPase

Purified (Na+, K+)-ATPase was studied by electron microscopy after thin sectioning, negative staining, and freeze-fracturing, particular emphasis being paid to the dimensions and frequencies of substructures in the membranes. Ultrathin sections show exclusively flat or cup-shaped membrane fragments which are triple-layered along much of their length and have diameters of 0.1-0.6 μm. Negative staining revealed a distinct substructure of particles with diameters between 30 and 50 A and with a frequency of 12,500 +/- 2,400 (SD) per μm(2). Comparisons with sizes of the protein components suggest that each surface particle contains as its major component one large catalytic chain with mol wt close to 100,000 and that two surface particles unite to form the unit of (Na+,K+)-ATPase which binds one molecule of ATP or ouabain. The further observations that the surface particles protrude from the membrane surface and are observed on both membrane surfaces in different patterns and degrees of clustering suggest that protein units span the membrane and are capable of lateral mobility. Freeze-fracturing shows intramembranous particles with diameters of 90-110 A and distributed on both concave and convex fracture faces with a frequency of 3,410 +/- 370 per μm(2) and 390 +/- 170 per μm(2), respectively. The larger diameters and three to fourfold smaller frequency of the intramembranous particles as compared to the surface particles seen after negative staining may reflect technical differences between methods, but it is more likely that the intramembranous particle is an oliogomer composed of two or even more of the protein units which form the surface particles.     

Solvent-induced free energy landscape and solute-solvent dynamic coupling in a multielement solute

Molecular dynamics simulations using a simple multielement model solute with internal degrees of freedom and accounting for solvent-induced interactions to all orders in explicit water are reported. The potential energy landscape of the solute is flat in vacuo. However, the sole untruncated solvent-induced interactions between apolar (hydrophobic) and charged elements generate a rich landscape of potential of mean force exhibiting typical features of protein landscapes. Despite the simplicity of our solute, the depth of minima in this landscape is not far in size from free energies that stabilize protein conformations. Dynamical coupling between configurational switching of the system and hydration reconfiguration is also elicited. Switching is seen to occur on a time scale two orders of magnitude longer than that of the reconfiguration time of the solute taken alone, or that of the unperturbed solvent. Qualitatively, these results are unaffected by a different choice of the water-water interaction potential. They show that already at an elementary level, solvent-induced interactions alone, when fully accounted for, can be responsible for configurational and dynamical features essential to protein folding and function.     

Empirical and deterministic accuracies of across-population genomic prediction

Background

Differences in linkage disequilibrium and in allele substitution effects of QTL (quantitative trait loci) may hinder genomic prediction across populations. Our objective was to develop a deterministic formula to estimate the accuracy of across-population genomic prediction, for which reference individuals and selection candidates are from different populations, and to investigate the impact of differences in allele substitution effects across populations and of the number of QTL underlying a trait on the accuracy.

Methods

A deterministic formula to estimate the accuracy of across-population genomic prediction was derived based on selection index theory. Moreover, accuracies were deterministically predicted using a formula based on population parameters and empirically calculated using simulated phenotypes and a GBLUP (genomic best linear unbiased prediction) model. Phenotypes of 1033 Holstein-Friesian, 105 Groninger White Headed and 147 Meuse-Rhine-Yssel cows were simulated by sampling 3000, 300, 30 or 3 QTL from the available high-density SNP (single nucleotide polymorphism) information of three chromosomes, assuming a correlation of 1.0, 0.8, 0.6, 0.4, or 0.2 between allele substitution effects across breeds. The simulated heritability was set to 0.95 to resemble the heritability of deregressed proofs of bulls.

Results

Accuracies estimated with the deterministic formula based on selection index theory were similar to empirical accuracies for all scenarios, while accuracies predicted with the formula based on population parameters overestimated empirical accuracies by ~25 to 30%. When the between-breed genetic correlation differed from 1, i.e. allele substitution effects differed across breeds, empirical and deterministic accuracies decreased in proportion to the genetic correlation. Using a multi-trait model, it was possible to accurately estimate the genetic correlation between the breeds based on phenotypes and high-density genotypes. The number of QTL underlying the simulated trait did not affect the accuracy.

Conclusions

The deterministic formula based on selection index theory estimated the accuracy of across-population genomic predictions well. The deterministic formula using population parameters overestimated the across-population genomic accuracy, but may still be useful because of its simplicity. Both formulas could accommodate for genetic correlations between populations lower than 1. The number of QTL underlying a trait did not affect the accuracy of across-population genomic prediction using a GBLUP method.     

Phagocytosis of bacteria by polymorphonuclear leukocytes: a freeze-fracture, scanning electron microscope, and thin-section investigation of membrane structure

The changes in membrane structure of rabbit polymorphonuclear (PMN) leukocytes during bacterial phagocytosis was investigated with scanning electron microscope (SEM), thin-section, and freeze-fracture techniques. SEM observations of bacterial attachment sites showed the involvement of limited areas of PMN membrane surface (0.01-0.25μm(2)). Frequently, these areas of attachment were located on membrane extensions. The membrane extensions were present before, during, and after the engulfment of bacteria, but were diminished in size after bacterial engulfment. In general, the results obtained with SEM and thin-section techniques aided in the interpretation of the three-dimensional freeze-fracture replicas. Freeze-fracture results revealed the PMN leukocytes had two fracture faces as determined by the relative density of intramembranous particles (IMP). Membranous extensions of the plasma membrane, lysosomes, and phagocytic vacuoles contained IMP's with a distribution and density similar to those of the plasma membrane. During phagocytosis, IMPs within the plasma membrane did not undergo a massive aggregation. In fact, structural changes within the membranes were infrequent and localized to regions such as the attachment sites of bacteria, the fusion sites on the plasma membrane, and small scale changes in the phagocytic vacuole membrane during membrane fusion. During the formation of the phagocytic vacuole, the IMPs of the plasma membrane appeared to move in with the lipid bilayer while maintaining a distribution and density of IMPs similar to those of the plasma membranes. Occasionally, IMPs were aligned to linear arrays within phagocytic vacuole membranes. This alignment might be due to an interaction with linearly arranged motile structures on the side of the phagocytic vacuole membranes. IMP-free regions were observed after fusion of lysosomes with the phagocytic vacuoles or plasma membrane. These IMP-free areas probably represent sites where membrane fusion occurred between lysosomal membrane and phagocytic vacuole membrane or plasma membrane. Highly symmetrical patterns of IMPs were not observed during lysosomal membrane fusion.