NMR Propolis pipeline 2
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 agroregions.
setwd("~/Dropbox")
library(metabolomicsUM)
source("Datasets/Propolis/NMR/scripts/propolis_metadata.R")
prop.nmr.metadata.file = "Datasets/Propolis/NMR/metadata/metadata_propolis_agro.csv"
prop.nmr.data.folder = "Datasets/Propolis/NMR/data"
get.metadata.agro(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)Agroregions metadata used, own grouping peaks algorithm used, removed peak groups with less than 25% of values, missing values imputation with low value, and log transformation.
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.0000Peak 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 transformed by log transformation:
prop.nmr.log = transform.data(prop.nmr.na, method="log")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.log = aov.all.vars(prop.nmr.log, "agroregions")
anova.prop.nmr.log[1:20,]## pvalues logs fdr
## 2.11 3.778168e-06 5.422719 0.0009143166
## 5.17 9.140972e-06 5.039008 0.0011060577
## 2.2 3.759972e-05 4.424815 0.0023624243
## 5.42 4.314793e-05 4.365040 0.0023624243
## 2.86 4.881042e-05 4.311487 0.0023624243
## 2.98 2.108492e-04 3.676028 0.0060283467
## 1.49 2.134870e-04 3.670629 0.0060283467
## 2.5 2.232885e-04 3.651134 0.0060283467
## 1.73 2.241947e-04 3.649375 0.0060283467
## 1.87 3.294498e-04 3.482211 0.0079726846
## 6.13 3.854040e-04 3.414084 0.0084788889
## 6.17 5.532303e-04 3.257094 0.0106587346
## 5.62 5.725766e-04 3.242166 0.0106587346
## 0.56 6.991081e-04 3.155456 0.0120845832
## 9.51 8.282134e-04 3.081858 0.0130755406
## 7.09 8.644986e-04 3.063236 0.0130755406
## 0.66 9.664493e-04 3.014821 0.0137576897
## 2.02 1.550503e-03 2.809528 0.0208456451
## 1.18 1.811065e-03 2.742066 0.0215165417
## 5.88 1.860076e-03 2.730469 0.0215165417
## tukey
## 2.11 Plain-Highlands; Plateau-Plain
## 5.17 Plain-Highlands; Plateau-Plain
## 2.2 Plain-Highlands; Plateau-Plain
## 5.42 Plain-Highlands; Plateau-Plain
## 2.86 Plain-Highlands; Plateau-Plain
## 2.98 Plain-Highlands; Plateau-Plain
## 1.49 Plain-Highlands; Plateau-Plain
## 2.5 Plain-Highlands; Plateau-Highlands; Plateau-Plain
## 1.73 Plain-Highlands; Plateau-Plain
## 1.87 Plain-Highlands; Plateau-Plain
## 6.13 Plain-Highlands; Plateau-Highlands; Plateau-Plain
## 6.17 Plain-Highlands; Plateau-Plain
## 5.62 Plain-Highlands; Plateau-Plain
## 0.56 Plain-Highlands; Plateau-Plain
## 9.51 Plain-Highlands; Plateau-Highlands
## 7.09 Plain-Highlands; Plateau-Plain
## 0.66 Plain-Highlands; Plateau-Plain
## 2.02 Plain-Highlands; Plateau-Plain
## 1.18 Plain-Highlands; Plateau-Plain
## 5.88 Plain-Highlands; Plateau-PlainA heatmap with the correlations between all the variables is shown below:
correl.prop.nmr.log = correlations.dataset(prop.nmr.log, method = "pearson")
heatmap(correl.prop.nmr.log, 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.log = clustering(prop.nmr.log, method = "hc", distance = "euclidean")
dendrogram.plot.col(prop.nmr.log, hc.prop.nmr.log, "agroregions")
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.log = clustering(prop.nmr.log, method = "kmeans", num.clusters = 4)
kmeans.plot(prop.nmr.log, kmeans.prop.nmr.log)
kmeans.df = kmeans.result.df(kmeans.prop.nmr.log, 4)
kmeans.df## cluster
## 1 1
## 2 2
## 3 3
## 4 4
## samples
## 1 CE_au CN_au IT_au AC_sm AN_sm BR_sm CE_sm CN_sm FP_sm IT_sm JB_sm PU_sm SJC_sm SJ_sm UR_sm VR_sm
## 2 AC_au BR_au XX_sm BR_sp CE_sp CN_sp DC_sp FP_sp IT_sp JB_sp PU_sp SA_sp SJC_sp SJ_sp UR_sp VR_sp AN_wi BR_wi CE_wi CN_wi PU_wi SA_wi UR_wi XX_wi
## 3 AN_au AC_sp AN_sp AC_wi DC_wi FP_wi JB_wi
## 4 DC_au JB_au PU_au SA_au SJ_au SJC_au UR_au VR_au XX_au DC_sm SA_sm SJ_wiPCA
Principal components analysis was performed on the data and some plots are shown below:
pca.analysis.result = pca.analysis.dataset(prop.nmr.log)
pca.pairs.plot(prop.nmr.log, pca.analysis.result, "agroregions")
pca.screeplot(pca.analysis.result)
pca.scoresplot2D(prop.nmr.log, pca.analysis.result, "agroregions", ellipses = T)
pca.kmeans.plot2D(prop.nmr.log, pca.analysis.result, kmeans.result = kmeans.prop.nmr.log, 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.log, c("pls", "J48", "JRip", "svmLinear", "rf"), "agroregions", "repeatedcv", num.folds = 10, num.repeats = 10, tunelength = 20, metric = "ROC")
ml.prop.nmr$performance## Accuracy Kappa Sensitivity Specificity ROC AccuracySD
## pls 0.7436190 0.4725440 0.6461111 0.8148333 0.8127407 0.1344380
## J48 0.6128333 0.2563641 0.5113889 0.7528333 0.6452523 0.1782898
## JRip 0.5615357 0.1089182 0.4252778 0.6990556 0.5829167 0.1635796
## svmLinear 0.6711429 0.2698650 0.5002778 0.7483333 0.8065370 0.1377969
## rf 0.7487143 0.4555116 0.6250000 0.8074444 0.7892685 0.1540536
## KappaSD SensitivitySD SpecificitySD ROCSD
## pls 0.2903929 0.1957675 0.09697075 0.1610395
## J48 0.3299237 0.2251870 0.11383720 0.1963986
## JRip 0.3146361 0.1983414 0.10703248 0.1561798
## svmLinear 0.3135632 0.1956032 0.10066341 0.1478531
## rf 0.3495580 0.2296561 0.11576995 0.1633538Also the confusion matrices and a plot using the first 3 PCs, showing the separation of the classes (agroregions) 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 Highlands Plain Plateau
## Highlands 8.3 0.9 1.7
## Plain 0.1 11.5 4.9
## Plateau 11.7 6.2 54.5
##
##
## $J48
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction Highlands Plain Plateau
## Highlands 9.9 2.9 8.8
## Plain 1.0 5.7 6.7
## Plateau 9.2 10.1 45.7
##
##
## $JRip
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction Highlands Plain Plateau
## Highlands 3.7 0.9 4.8
## Plain 1.4 6.3 10.2
## Plateau 15.0 11.5 46.1
##
##
## $svmLinear
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction Highlands Plain Plateau
## Highlands 4.5 0.0 2.4
## Plain 1.2 7.1 3.1
## Plateau 14.6 11.6 55.6
##
##
## $rf
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction Highlands Plain Plateau
## Highlands 10.0 0.0 1.6
## Plain 0.0 9.2 3.9
## Plateau 10.1 9.6 55.7pls.model = ml.prop.nmr$final.models$pls
pca.plot.3d(prop.nmr.log, pls.model, "agroregions")
And the variable importance in the four season classes for all models:
summary.var.importance(ml.prop.nmr, 10)## $pls
## Highlands Plain Plateau Mean
## 7.5 64.64380 31.751225 100.00000 65.46501
## 3.01 59.44834 36.683049 79.35365 58.49501
## 1.15 54.38420 33.393921 84.46511 57.41441
## 2.86 49.30220 50.688959 72.22487 57.40535
## 1.68 71.07022 13.170967 87.12019 57.12046
## 9.51 58.69869 7.202725 89.18363 51.69502
## 0.56 31.33796 60.296445 61.18120 50.93853
## 3.07 48.28552 41.381811 60.97409 50.21381
## 1.64 43.75482 36.874803 65.75909 48.79624
## 2.65 58.92064 9.622731 74.65386 47.73241
##
## $J48
## Highlands Plain Plateau Mean
## 2.11 100.00000 100.00000 88.92950 96.30983
## 5.17 96.24021 96.24021 85.16971 92.55004
## 1.87 100.00000 100.00000 77.44125 92.48042
## 1.9 97.49347 97.49347 72.56745 89.18480
## 5.28 94.98695 94.98695 77.44125 89.13838
## 5.42 96.24021 96.24021 74.72585 89.06876
## 5.62 96.24021 96.24021 73.05483 88.51175
## 2.2 93.73368 93.73368 78.06789 88.51175
## 2.86 75.56136 91.64491 91.64491 86.28372
## 0.66 86.21410 86.21410 86.21410 86.21410
##
## $JRip
## Overall Mean
## 2.11 100 100
## 3.9 100 100
## 3.99 100 100
## 0.27 0 0
## 0.31 0 0
## 0.34 0 0
## 0.41 0 0
## 0.44 0 0
## 0.48 0 0
## 0.5 0 0
##
## $svmLinear
## Highlands Plain Plateau Mean
## 2.11 100.00000 100.00000 88.92950 96.30983
## 5.17 96.24021 96.24021 85.16971 92.55004
## 1.87 100.00000 100.00000 77.44125 92.48042
## 1.9 97.49347 97.49347 72.56745 89.18480
## 5.28 94.98695 94.98695 77.44125 89.13838
## 5.42 96.24021 96.24021 74.72585 89.06876
## 5.62 96.24021 96.24021 73.05483 88.51175
## 2.2 93.73368 93.73368 78.06789 88.51175
## 2.86 75.56136 91.64491 91.64491 86.28372
## 0.66 86.21410 86.21410 86.21410 86.21410
##
## $rf
## Overall Mean
## 2.86 100.00000 100.00000
## 5.31 78.53738 78.53738
## 5.17 61.84556 61.84556
## 3.25 56.91384 56.91384
## 2.65 54.08224 54.08224
## 2.11 53.00751 53.00751
## 3.19 44.71398 44.71398
## 6.17 39.66502 39.66502
## 6.13 39.19544 39.19544
## 6.75 37.34783 37.34783FEATURE 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.log, "agroregions", 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.5701 0.1757 0.1747 0.3265
## 4 0.6241 0.2606 0.1710 0.3287
## 8 0.6746 0.3080 0.1609 0.3577
## 16 0.7069 0.3981 0.1668 0.3576
## 32 0.7560 0.4687 0.1599 0.3626 *
## 64 0.7400 0.4291 0.1422 0.3352
## 242 0.7464 0.4407 0.1396 0.3355
##
## The top 5 variables (out of 32):
## X6.17, X2.11, X5.21, X6.13, X2.86plot(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.log, "agroregions", 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.7432 0.4371 0.1302 0.3274
##
## Using the training set, 70 variables were selected:
## X0.5, X0.56, X0.61, X0.66, X0.99...
##
## During resampling, the top 5 selected variables (out of a possible 117):
## X0.56 (100%), X0.61 (100%), X0.66 (100%), X1.02 (100%), X1.18 (100%)
##
## On average, 65.4 variables were selected (min = 50, max = 91)