*biosigner*: A new method for signature discovery from omics data

In biosigner: Signature discovery from omics data

knitr::opts_chunk$set(fig.width = 6, fig.height = 6, fig.path = 'figures/')

Introduction

High-throughput, non-targeted, technologies such as transcriptomics, proteomics and metabolomics, are widely used to discover molecules which allow to efficiently discriminate between biological or clinical conditions of interest (e.g., disease vs control states). Powerful machine learning approaches such as Partial Least Square Discriminant Analysis (PLS-DA), Random Forest (RF) and Support Vector Machines (SVM) have been shown to achieve high levels of prediction accuracy. Feature selection, i.e., the selection of the few features (i.e., the molecular signature) which are of highest discriminating value, is a critical step in building a robust and relevant classifier [@Guyon2003]: First, dimension reduction is usefull to limit the risk of overfitting and increase the prediction stability of the model; second, intrepretation of the molecular signature is facilitated; third, in case of the development of diagnostic product, a restricted list is required for the subsequent validation steps [@Rifai2006].

Since the comprehensive analysis of all combinations of features is not computationally tractable, several selection techniques have been described, including filter (e.g., p-values thresholding), wrapper (e.g., recursive feature elimination), and embedded (e.g., sparse PLS) approaches [@Saeys2007]. The major challenge for such methods is to be fast and extract restricted and stable molecular signatures which still provide high performance of the classifier [@Gromski2014; @Determan2015].

The biosigner package

The biosigner package implements a new wrapper feature selection algorithm:

1. the dataset is split into training and testing subsets (by bootstraping, controling class proportion),

2. model is trained on the training set and balanced accuracy is evaluated on the test set,

3. the features are ranked according to their importance in the model,

4. the relevant feature subset at level f is found by a binary search: a feature subset is considered relevant if and only if, when randomly permuting the intensities of other features in the test subsets, the proportion of increased or equal prediction accuracies is lower than a defined threshold f,

5. the dataset is restricted to the selected features and steps 1 to 4 are repeated until the selected list of features is stable.

Three binary classifiers have been included in biosigner, namely PLS-DA, RF and SVM, as the performances of each machine learning approach may vary depending on the structure of the dataset [@Determan2015]. The algorithm returns the tier of each feature for the selected classifer(s): tier S corresponds to the final signature, i.e., features which have been found significant in all the selection steps; features with tier A have been found significant in all but the last selection, and so on for tier B to D. Tier E regroup all previous round of selection.

As for a classical classification algorithm, the biosign method takes as input the x samples times features data frame (or matrix) of intensities, and the y factor (or character vector) of class labels (note that only binary classification is currently available). It returns the signature (signatureLs: selected feature names) and the trained model (modelLs) for each of the selected classifier. The plot method for biosign objects enable to visualize the individual boxplots of the selected features. Finally, the predict method allows to apply the trained classifier(s) on new datasets.

The algorithm has been successfully applied to transcriptomics and metabolomics data [@Rinaudo2016; see also the Hands-on section below).

Hands-on

Loading

We first load the biosigner package:

library(biosigner)

We then use the diaplasma metabolomics dataset [@Rinaudo2016] which results from the analysis of plasma samples from 69 diabetic patients were analyzed by reversed-phase liquid chromatography coupled to high-resolution mass spectrometry (LC-HRMS; Orbitrap Exactive) in the negative ionization mode. The raw data were pre-processed with XCMS and CAMERA (5,501 features), corrected for signal drift, log10 transformed, and annotated with an in-house spectral database. The patient's age, body mass index, and diabetic type are recorded [@Rinaudo2016].

data(diaplasma)

We attach diaplasma to the search path and display a summary of the content of the dataMatrix, sampleMetadata and variableMetadata with the view function from the (imported) ropls package:

attach(diaplasma)
library(ropls)
ropls::view(dataMatrix)
ropls::view(sampleMetadata, standardizeL = TRUE)
ropls::view(variableMetadata, standardizeL = TRUE)

We see that the diaplasma list contains three objects:

1. dataMatrix: 69 samples x 5,501 matrix of numeric type containing the intensity profiles (log10 transformed),

2. sampleMetadata: a 69 x 3 data frame, with the patients'

   • type: diabetic type, factor

   • age: numeric

   • bmi: body mass index, numeric

3. variableMetadata: a 5,501 x 8 data frame, with the median m/z ('mzmed', numeric) and the median retention time in seconds ('rtmed', numeric) from XCMS, the 'isotopes' (character), 'adduct' (character) and 'pcgroups' (numeric) annotations from CAMERA, the names of the m/z and RT matching compounds from an in-house spectra of commercial metabolites ('name_hmdb', character), and the p-values resulting from the non-parametric hypothesis testing of difference in medians between types ('type_wilcox_fdr', numeric), and correlation with age ('age_spearman_fdr', numeric) and body mass index ('bmi_spearman_fdr', numeric), all corrected for multiple testing (False Discovery Rate).

We can observe that the 3 clinical covariates (diabetic type, age, and bmi) are stronlgy associated:

with(sampleMetadata,
plot(age, bmi, cex = 1.5, col = ifelse(type == "T1", "blue", "red"), pch = 16))
legend("topleft", cex = 1.5, legend = paste0("T", 1:2),
text.col = c("blue", "red"))

Figure 1: age, body mass index (bmi), and diabetic type of the patients from the diaplasma cohort.

Molecular signatures

Let us look for signatures of type in the diaplasma dataset by using the biosign method. To speed up computations in this demo vignette, we restrict the number of features (from 5,501 to about 500) and the number of bootstraps (5 instead of 50 [default]); the selection on the whole dataset, 50 bootstraps, and the 3 classifiers, takes around 10 min.

featureSelVl <- variableMetadata[, "mzmed"] >= 450 &
variableMetadata[, "mzmed"] < 500
sum(featureSelVl)
dataMatrix <- dataMatrix[, featureSelVl]
variableMetadata <- variableMetadata[featureSelVl, ]

diaSign <- biosigner::biosign(dataMatrix, sampleMetadata[, "type"], bootI = 5)

Figure 2: Relevant signatures for the PLS-DA, Random Forest, and SVM classifiers extracted from the diaplasma dataset. The S tier corresponds to the final metabolite signature, i.e., metabolites which passed through all the selection steps.

The arguments are:

• x: the numerical matrix (or data frame) of intensities (samples as rows, variables as columns),

• y: the factor (or character) specifying the sample labels from the 2 classes,

• methodVc: the classifier(s) to be used; here, the default all value means that all classifiers available (plsda, randomforest, and svm) are selected,

• bootI: the number of bootstraps is set to 5 to speed up computations when generating this vignette; we however recommend to keep the default 50 value for your analyzes (otherwise signatures may be less stable).

• The set.seed argument ensures that the results from this vignette can be reproduced exactly; by choosing alternative seeds (and the default bootI = 50), similar signatures are obtained, showing the stability of the selection.

Note:

• If some features from the x matrix/data frame contain missing values (NA), these features will be removed prior to modeling with Random Forest and SVM (in contrast, the NIPALS algorithm from PLS-DA can handle missing values),

The resulting signatures for the 3 selected classifiers are both printed and plotted as tiers from S, A, up to E by decreasing relevance. The (S) tier corresponds to the final signature, i.e. features which passed through all the backward selection steps. In contrast, features from the other tiers were discarded during the last (A) or previous (B to E) selection rounds.

Note that tierMaxC = 'A' argument in the print and plot methods can be used to view the features from the larger S+A signatures (especially when no S features have been found, or when the performance of the S model is much lower than the S+A model).

The performance of the model built with the input dataset (balanced accuracy: mean of the sensitivity and specificity), or the subset restricted to the S or S+A signatures are shown. We see that with 1 to 5 S feature signatures (i.e., less than 1% of the input), the 3 classifiers achieve good performances (even higher than the full Random Forest and SVM models). Furthermore, reducing the number of features decreases the risk of building non-significant models (i.e., models which do not perform significantly better than those built after randomly permuting the labels). The signatures from the 3 classifiers have some distinct features, which highlights the interest of comparing various machine learning approaches.

The individual boxplots of the features from the complete signature can be visualized with:

biosigner::plot(diaSign, typeC = "boxplot")

Figure 3: Individual boxplots of the features selected for at least one of the classification methods. Features selected for a single classifier are colored (red for PLS-DA, green for Random Forest and blue for SVM).

Let us see the metadata of the complete signature:

variableMetadata[getSignatureLs(diaSign)[["complete"]], ]

Predictions

Let us split the dataset into a training (the first 4/5th of the 183 samples) and a testing subsets, and extract the relevant features from the training subset:

trainVi <- 1:floor(0.8 * nrow(dataMatrix))
testVi <- setdiff(1:nrow(dataMatrix), trainVi)

diaTrain <- biosigner::biosign(dataMatrix[trainVi, ], sampleMetadata[trainVi, "type"],
bootI = 5)

Figure 4: Signatures from the training data set.

We extract the fitted types on the training dataset restricted to the S signatures:

diaFitDF <- biosigner::predict(diaTrain)

We then print the confusion tables for each classifier:

lapply(diaFitDF, function(predFc) table(actual = sampleMetadata[trainVi,
"type"], predicted = predFc))

and the corresponding balanced accuracies:

sapply(diaFitDF, function(predFc) {
conf <- table(sampleMetadata[trainVi, "type"], predFc)
conf <- sweep(conf, 1, rowSums(conf), "/")
round(mean(diag(conf)), 3)
})

Note that these values are slightly different from the accuracies returned by biosign because the latter are computed by using the resampling scheme selected by the bootI (or crossvalI) arguments:

round(biosigner::getAccuracyMN(diaTrain)["S", ], 3)

Finally, we can compute the performances on the test subset:

diaTestDF <- biosigner::predict(diaTrain, newdata = dataMatrix[testVi, ])
sapply(diaTestDF, function(predFc) {
conf <- table(sampleMetadata[testVi, "type"], predFc)
conf <- sweep(conf, 1, rowSums(conf), "/")
round(mean(diag(conf)), 3)
})

Working on ExpressionSet omics objects from bioconductor

The ExpressionSet class from the Biobase bioconductor package has been developed to conveniently handle preprocessed omics objects, including the variables x samples matrix of intensities, and data frames containing the sample and variable metadata [@Huber2015]. The matrix and the two data frames can be accessed by the exprs, pData and fData respectively (note that the data matrix is stored in the object with samples in columns).

The biosign method can be applied to an ExpressionSet object, by using the object as the x argument, and by indicating as the y argument the name of the phenoData column to be used as the response.

In the example below, we will first build an ExpressionSet object from the diaplasma data set, we display its content with the view method from the ropls package, and we then perform the signature discovery:

library(Biobase)

diaSet <- Biobase::ExpressionSet(assayData = t(dataMatrix),
                                 phenoData = new("AnnotatedDataFrame",
                                                 data = sampleMetadata),
                                 featureData = new("AnnotatedDataFrame",
                                                 data = variableMetadata))

ropls::view(diaSet)

biosigner::biosign(diaSet, "type", bootI = 5)

Note: Because of identical names in the two packages, the combine method from package Biobase is replaced by the combine method from package randomForest.

The ExpressionSet can be displayed

ropls::view(diaSet)

Before moving to new data sets, we detach diaplasma from the search path:

detach(diaplasma)

Extraction of biomarker signatures from other omics datasets

In this section, biosign is applied to two metabolomics and one transcriptomics data sets. Please refer to @Rinaudo2016 for a full discussion of the methods and results.

Physiological variations of the human urine metabolome (metabolomics)

The sacurine LC-HRMS dataset from the dependent ropls package can also be used [@Thevenot2015]: Urine samples from a cohort of 183 adults were analyzed by using an LTQ Orbitrap in the negative ionization mode. A total of 109 metabolites were identified or annotated at the MSI level 1 or 2. Signal drift and batch effect were corrected, and each urine profile was normalized to the osmolality of the sample. Finally, the data were log10 transformed (see the ropls vignette for further details and examples).

We can for instance look for signatures of the gender:

data(sacurine)
sacSign <- biosigner::biosign(sacurine[["dataMatrix"]],
sacurine[["sampleMetadata"]][, "gender"],
methodVc = "plsda")

Figure 5: PLS-DA signature from the 'sacurine' data set.

Apples spikes with known compounds (metabolomics)

The spikedApples dataset was obtained by LC-HRMS analysis (SYNAPT Q-TOF, Waters) of one control and three spiked groups of 10 apples each. The spiked mixtures consists in 2 compounds which were not naturally present in the matrix and 7 compounds aimed at achieving a final increase of 20%, 40% or 100% of the endogeneous concentrations. The authors identified 22 features (out of the 1,632 detected in the positive ionization mode; i.e. 1.3%) which came from the spiked compounds. The dataset is included in the BioMark R package [@Franceschi2012]. Let us use the control and group1 samples (20 in total) in this study.

library(BioMark)
data(SpikePos)
group1Vi <- which(SpikePos[["classes"]] %in% c("control", "group1"))
appleMN <- SpikePos[["data"]][group1Vi, ]
spikeFc <- factor(SpikePos[["classes"]][group1Vi])
annotDF <- SpikePos[["annotation"]]
rownames(annotDF) <- colnames(appleMN)

We can check, by using the opls method from the ropls package for multivariate analysis, that:

1. no clear separation can be observed by PCA:

biomark.pca <- ropls::opls(appleMN, fig.pdfC = "none")
ropls::plot(biomark.pca, parAsColFcVn = spikeFc)

1. PLS-DA modeling with the full dataset is not significant (as seen on the top left plot: 7 out of 20 models trained after random permutations of the labels have Q2 values higher than the model trained with the true labels):

biomark.pls <- ropls::opls(appleMN, spikeFc)

Let us now extract the molecular signatures:

appleSign <- biosigner::biosign(appleMN, spikeFc)

The 449.1/327 corresponds to the Cyanidin-3-galactoside (absent in the control; @Franceschi2012).

annotDF <- SpikePos[["annotation"]]
rownames(annotDF) <- colnames(appleMN)
annotDF[getSignatureLs(appleSign)[["complete"]], c("adduct", "found.in.standards")]

Bone marrow from acute leukemia patients (transcriptomics)

Samples from 47 patients with acute lymphoblastic leukemia (ALL) and 25 patients with acute myeloid leukemia (AML) have been analyzed using Affymetrix Hgu6800 chips resulting in expression data of 7,129 gene probes [@Golub1999]. The golub dataset is available in the golubEsets package from Bioconductor. Let us compute for example the SVM signature (to speed up this demo example, the number of features is restricted to 500):

library(golubEsets)
data(Golub_Merge)
golubMN <- t(exprs(Golub_Merge))
leukemiaFc <- pData(Golub_Merge)[["ALL.AML"]]
table(leukemiaFc)
varSubVi <- 1501:2000
golubSign <- biosigner::biosign(golubMN[, varSubVi], leukemiaFc, methodVc = "svm")

Figure 6: SVM signature from the golub data set.

The computation results in a signature of 4 features only and a sparse SVM model performing even better (95.9% accuracy) than the model trained on the dataset of 500 variables (95.5% accuracy).

The hu6800.db bioconductor package can be used to get the annotation of the selected probes

library(hu6800.db)
sapply(biosigner::getSignatureLs(golubSign)[["complete"]],
       function(probeC)
       get(probeC, env = hu6800GENENAME))

Cystatin C is part of the 50 gene signature selected by Golub and colleagues on the basis of a metric derived from the Student's statistic of mean differences between the AML and ALL groups [@Golub1999]. Interestingly, the third probe, myeloperoxidase, is a cytochemical marker for the diagnosis (and also potentially the prognosis) of acute myeloid leukemia (AML).

rm(list = ls())

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

sessionInfo()

References

biosigner documentation built on Nov. 24, 2020, 2 a.m.