NMR Propolis pipeline 3

Propolis peak list data was read and stored in a list, containing 2 elements, the dataset consisting in a list with the samples with their ppm intensities being the elements. The propolis metadata consists on the seasons and regions. Seasons were obtained by the last part of the file name, for example the file name “AC_au” means that the season of the sample is autumn (au). Stations were obtained by the first part of the file name and then assigned to the specific region, for example again the file name “AC_au”, the station is “AC” that was assigned to West region.

setwd("~/Dropbox")
library(metabolomicsUM)
source("Datasets/Propolis/NMR/scripts/propolis_metadata.R")

prop.nmr.metadata.file = "Datasets/Propolis/NMR/metadata/metadata_propolis.csv"
prop.nmr.data.folder = "Datasets/Propolis/NMR/data"

get.metadata(prop.nmr.data.folder, write.file = TRUE, file.name = prop.nmr.metadata.file)
prop.nmr.metadata = read.metadata(prop.nmr.metadata.file)

peaks.lists = read.csvs.folder(prop.nmr.data.folder)

Seasons metadata used, own grouping peaks algorithm used, removed peak groups with less than 25% of values, missing values imputation with low value, and autoscaling.

PREPROCESSING

Own grouping peaks algorithm used with step = 0.03:

# removing resonances in selected regions
peaks.lists = remove.peaks.interval.sample.list(peaks.lists, 0, 0.19)
peaks.lists = remove.peaks.interval.sample.list(peaks.lists, 3.29, 3.31)
peaks.lists = remove.peaks.interval.sample.list(peaks.lists, 4.85, 5)

#group peaks
prop.nmr.ds = group.peaks(peaks.lists, type = "nmr-peaks", metadata = prop.nmr.metadata, description = "NMR propolis", label.x = "ppm", label.values = "intensity")
sum.dataset(prop.nmr.ds)
## Dataset summary:
## Valid dataset
## Description:  NMR propolis 
## Type of data:  nmr-peaks 
## Number of samples:  59 
## Number of data points 293 
## Number of metadata variables:  2 
## Label of x-axis values:  ppm 
## Label of data points:  intensity 
## Number of missing values in data:  5376 
## Mean of data values:  0.09016594 
## Median of data values:  0.0287 
## Standard deviation:  0.1904829 
## Range of values:  0 10 
## Quantiles: 
##      0%     25%     50%     75%    100% 
##  0.0000  0.0081  0.0287  0.0929 10.0000

Peak groups with less than 25% of values were removed:

nsamps = num.samples(prop.nmr.ds)
prop.nmr.ds = remove.variables.by.nas(prop.nmr.ds,  0.75*nsamps)

There are 2659 missing values found in the dataset, which will be replaced with a low value (0.00005).

prop.nmr.na = missingvalues.imputation(prop.nmr.ds, method="value", value = 0.00005)

The data was scaled with autoscaling:

prop.nmr.scal = scaling(prop.nmr.na, "auto")

UNIVARIATE TESTS

An analysis of variance (ANOVA) was conducted over the data with tukey test also, and this is the top 10 results ordered by p-value:

anova.prop.nmr.scal = aov.all.vars(prop.nmr.scal, "seasons")
anova.prop.nmr.scal[1:20,]
##           pvalues     logs          fdr                      tukey
## 5.99 2.314264e-07 6.635587 5.600519e-05        sp-au; sp-sm; wi-sp
## 4.66 1.746928e-06 5.757725 2.113783e-04        sm-au; sp-sm; wi-sm
## 4.5  1.843384e-05 4.734384 1.486996e-03 sp-au; wi-au; sp-sm; wi-sm
## 4.45 2.617632e-05 4.582091 1.507751e-03        sp-au; sp-sm; wi-sm
## 4.55 3.115188e-05 4.506516 1.507751e-03        sm-au; sp-sm; wi-sm
## 6.28 5.288837e-05 4.276640 2.133164e-03        sm-au; sp-au; wi-au
## 4.58 6.524368e-05 4.185462 2.170869e-03        sm-au; sp-sm; wi-sm
## 4.63 7.256888e-05 4.139250 2.170869e-03        sm-au; sp-sm; wi-sm
## 4.08 8.073480e-05 4.092939 2.170869e-03               sp-au; wi-au
## 0.79 3.135256e-04 3.503727 7.587319e-03               sm-au; sp-au
## 4.05 4.321872e-04 3.364328 9.508119e-03               sp-sm; wi-sm
## 4.71 4.751568e-04 3.323163 9.582328e-03        sm-au; sp-sm; wi-sm
## 3.9  5.205975e-04 3.283498 9.691123e-03        wi-au; wi-sm; wi-sp
## 6.25 5.828928e-04 3.234411 1.007572e-02        sm-au; sp-au; wi-au
## 6.2  8.174775e-04 3.087524 1.318864e-02               sm-au; wi-au
## 4.53 9.637365e-04 3.016042 1.435688e-02               sp-sm; wi-sm
## 1.71 1.008541e-03 2.996306 1.435688e-02               sm-au; sp-au
## 6.03 1.484145e-03 2.828524 1.995351e-02        sp-au; sp-sm; wi-sp
## 7.12 1.941417e-03 2.711881 2.472594e-02               sm-au; sp-au
## 4.02 2.049582e-03 2.688335 2.472594e-02               sp-au; sp-sm

A heatmap with the correlations between all the variables is shown below:

correl.prop.nmr.scal = correlations.dataset(prop.nmr.scal, method = "pearson")
heatmap(correl.prop.nmr.scal, col =  topo.colors(256))

CLUSTERING

Hierarchical clustering with euclidean distance and complete method was performed on the data and the resulting dendrogram is shown below:

hc.prop.nmr.scal = clustering(prop.nmr.scal, method = "hc", distance = "euclidean")
dendrogram.plot.col(prop.nmr.scal, hc.prop.nmr.scal, "seasons")

K-Means was performed on the data also with 4 centers and the results and the plot giving for each cluster the median of the samples in blue, and in grey the values of all samples in that cluster are shown below:

kmeans.prop.nmr.scal = clustering(prop.nmr.scal, method = "kmeans", num.clusters = 4)
kmeans.plot(prop.nmr.scal, kmeans.prop.nmr.scal)

kmeans.df = kmeans.result.df(kmeans.prop.nmr.scal, 4)
kmeans.df
##   cluster
## 1       1
## 2       2
## 3       3
## 4       4
##                                                                                                                           samples
## 1        AN_au JB_au AC_sm AN_sm BR_sm CE_sm CN_sm FP_sm JB_sm SJC_sm SJ_sm UR_sm VR_sm AC_sp AN_sp DC_sp JB_sp DC_wi FP_wi JB_wi
## 2 BR_au CE_au CN_au IT_au SJC_au XX_au IT_sm PU_sm XX_sm BR_sp CE_sp CN_sp IT_sp PU_sp SJC_sp BR_wi CE_wi CN_wi PU_wi SA_wi XX_wi
## 3                                                         AC_au DC_au PU_au SA_au VR_au DC_sm SA_sm FP_sp SA_sp VR_sp AC_wi AN_wi
## 4                                                                                             SJ_au UR_au SJ_sp UR_sp SJ_wi UR_wi

PCA

Principal components analysis was performed on the data and some plots are shown below:

pca.analysis.result = pca.analysis.dataset(prop.nmr.scal)

pca.pairs.plot(prop.nmr.scal, pca.analysis.result, "seasons")

pca.screeplot(pca.analysis.result)

pca.scoresplot2D(prop.nmr.scal, pca.analysis.result, "seasons", ellipses = T)

pca.kmeans.plot2D(prop.nmr.scal, pca.analysis.result, kmeans.result = kmeans.prop.nmr.scal, ellipses = T)

MACHINE LEARNING

For classification models and prediction the following parameters were used: - models: PLS, J48, JRip, SVM and Random Forests - validation method: repeated cross-validation - number of folds: 5 - number of repeats: 10

Below are some results with the best tune for each model:

ml.prop.nmr = train.models.performance(prop.nmr.scal, c("pls", "J48", "JRip", "svmLinear", "rf"), "seasons", "repeatedcv", num.folds = 10, num.repeats = 10, tunelength = 20, metric = "ROC")
ml.prop.nmr$performance
##            Accuracy     Kappa Sensitivity Specificity       ROC AccuracySD
## pls       0.7977381 0.7241769     0.79250   0.9325833 0.9341667  0.1554731
## J48       0.7396071 0.6471210     0.72875   0.9145833 0.8527813  0.1435735
## JRip      0.5464762 0.3943389     0.53500   0.8514583 0.7340208  0.1920258
## svmLinear 0.7352619 0.6393040     0.73125   0.9113333 0.8940417  0.1663753
## rf        0.8316905 0.7714573     0.83125   0.9439167 0.9398125  0.1690439
##             KappaSD SensitivitySD SpecificitySD      ROCSD
## pls       0.2111746     0.1696885    0.05206651 0.09613388
## J48       0.1933682     0.1549428    0.04816534 0.11855880
## JRip      0.2467743     0.1889297    0.06200668 0.13945184
## svmLinear 0.2257624     0.1797647    0.05573229 0.10712427
## rf        0.2279423     0.1716812    0.05675603 0.08480285

Also the confusion matrices and a plot using the first 3 PCs, showing the separation of the four classes (seasons) are shown below:

ml.prop.nmr$confusion.matrices
## $pls
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix 
## 
## (entries are percentages of table totals)
##  
##           Reference
## Prediction   au   sm   sp   wi
##         au 18.8  0.0  0.8  1.7
##         sm  0.8 25.5  2.7  1.6
##         sp  4.3  1.7 19.4  2.8
##         wi  1.4  0.0  2.3 16.0
## 
## 
## $J48
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix 
## 
## (entries are percentages of table totals)
##  
##           Reference
## Prediction   au   sm   sp   wi
##         au 17.0  0.0  1.2  8.1
##         sm  0.0 25.5  0.0  0.0
##         sp  2.5  1.5 22.0  4.6
##         wi  5.8  0.2  2.1  9.4
## 
## 
## $JRip
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix 
## 
## (entries are percentages of table totals)
##  
##           Reference
## Prediction   au   sm   sp   wi
##         au  9.5  3.4  2.2  4.1
##         sm  7.2 20.7  3.4  5.8
##         sp  4.6  2.4 18.1  5.6
##         wi  4.2  0.8  1.7  6.4
## 
## 
## $svmLinear
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix 
## 
## (entries are percentages of table totals)
##  
##           Reference
## Prediction   au   sm   sp   wi
##         au 19.3  0.0  4.5  3.3
##         sm  1.1 24.7  2.3  0.2
##         sp  2.7  2.4 15.2  4.3
##         wi  2.4  0.0  3.3 14.2
## 
## 
## $rf
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix 
## 
## (entries are percentages of table totals)
##  
##           Reference
## Prediction   au   sm   sp   wi
##         au 19.3  0.0  1.5  1.9
##         sm  0.4 25.5  0.0  0.0
##         sp  1.9  1.7 22.0  3.9
##         wi  3.7  0.0  1.7 16.3
pls.model = ml.prop.nmr$final.models$pls
pca.plot.3d(prop.nmr.scal, pls.model, "seasons")

And the variable importance in the four season classes for all models:

summary.var.importance(ml.prop.nmr, 10)
## $pls
##            au       sm        sp       wi     Mean
## 5.99 45.55981 32.95526 100.00000 87.35786 66.46823
## 3.9  40.66413 13.55412  57.40562 83.78794 48.85295
## 6.03 37.55165 20.99708  64.15965 71.13740 48.46145
## 9.9  21.66274 33.19419  68.71583 59.48347 45.76406
## 4.66 51.25612 59.81745  35.07289 27.65622 43.45067
## 4.84 48.17582 30.76596  26.25941 67.32524 43.13161
## 6.28 42.91731 30.03494  50.37414 47.96650 42.82323
## 4.63 49.35203 55.26174  36.51380 27.04217 42.04244
## 1.27 34.92026 33.07798  52.13195 46.07721 41.55185
## 5.78 29.06201 18.36214  56.28591 61.33341 41.26086
## 
## $J48
##             au  sm  sp        wi      Mean
## 4.66 100.00000 100 100 100.00000 100.00000
## 4.58  96.00000 100 100 100.00000  99.00000
## 4.55  83.00000 100 100 100.00000  95.75000
## 4.5   89.00000 100 100  92.69231  95.42308
## 4.38  89.00000 100 100  88.84615  94.46154
## 4.08  89.00000 100 100  88.07692  94.26923
## 4.17  89.00000 100 100  88.07692  94.26923
## 4.31  83.66667 100 100  91.92308  93.89744
## 4.45  83.66667 100 100  91.15385  93.70513
## 4.63  88.33333  95  95  95.00000  93.33333
## 
## $JRip
##      Overall Mean
## 0.76     100  100
## 4.02     100  100
## 5.6      100  100
## 5.62     100  100
## 5.99     100  100
## 0.27       0    0
## 0.31       0    0
## 0.34       0    0
## 0.41       0    0
## 0.44       0    0
## 
## $svmLinear
##             au  sm  sp        wi      Mean
## 4.66 100.00000 100 100 100.00000 100.00000
## 4.58  96.00000 100 100 100.00000  99.00000
## 4.55  83.00000 100 100 100.00000  95.75000
## 4.5   89.00000 100 100  92.69231  95.42308
## 4.38  89.00000 100 100  88.84615  94.46154
## 4.08  89.00000 100 100  88.07692  94.26923
## 4.17  89.00000 100 100  88.07692  94.26923
## 4.31  83.66667 100 100  91.92308  93.89744
## 4.45  83.66667 100 100  91.15385  93.70513
## 4.63  88.33333  95  95  95.00000  93.33333
## 
## $rf
##         Overall       Mean
## 4.66 100.000000 100.000000
## 5.99  39.408282  39.408282
## 4.58  17.862694  17.862694
## 6.03  13.799246  13.799246
## 2.32  13.280028  13.280028
## 6.28  12.929364  12.929364
## 4.71  11.137461  11.137461
## 4.08  10.251426  10.251426
## 4.38   8.331068   8.331068
## 4.2    7.757915   7.757915

FEATURE SELECTION

Using recursive feature selection, various subsets were used with random forests classifier. The results are shown below:

feature.selection.result = feature.selection(prop.nmr.scal, "seasons", method="rfe", functions = rfFuncs, validation = "repeatedcv", repeats = 5, subsets = 2^(1:6))
feature.selection.result
## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (10 fold, repeated 5 times) 
## 
## Resampling performance over subset size:
## 
##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
##          2   0.6542 0.5341     0.1342  0.1729         
##          4   0.7002 0.5939     0.1235  0.1583         
##          8   0.7227 0.6246     0.1216  0.1604         
##         16   0.7798 0.7042     0.1593  0.2113         
##         32   0.8261 0.7646     0.1156  0.1547         
##         64   0.8444 0.7896     0.1057  0.1419        *
##        242   0.8219 0.7581     0.1295  0.1759         
## 
## The top 5 variables (out of 64):
##    X4.66, X4.58, X4.08, X4.17, X4.5
plot(feature.selection.result, type=c("g","o"))

Also selection by filter was used with the results shown below:

feature.selection.result2 = feature.selection(prop.nmr.scal, "seasons", method="filter", functions = rfSBF, validation = "repeatedcv", repeats = 5, subsets = 2^(1:6))
feature.selection.result2
## 
## Selection By Filter
## 
## Outer resampling method: Cross-Validated (10 fold, repeated 5 times) 
## 
## Resampling performance:
## 
##  Accuracy  Kappa AccuracySD KappaSD
##    0.8256 0.7629     0.1474  0.1989
## 
## Using the training set, 84 variables were selected:
##    X0.41, X0.54, X0.66, X0.79, X1.12...
## 
## During resampling, the top 5 selected variables (out of a possible 157):
##    X0.79 (100%), X1.15 (100%), X1.71 (100%), X1.96 (100%), X3.87 (100%)
## 
## On average, 76.2 variables were selected (min = 54, max = 100)