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)