NMR Propolis pipeline 1
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 no normalization.
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)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.na = aov.all.vars(prop.nmr.na, "seasons")
anova.prop.nmr.na[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-smA heatmap with the correlations between all the variables is shown below:
correl.prop.nmr.na = correlations.dataset(prop.nmr.na, method = "pearson")
heatmap(correl.prop.nmr.na, 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.na = clustering(prop.nmr.na, method = "hc", distance = "euclidean")
dendrogram.plot.col(prop.nmr.na, hc.prop.nmr.na, "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.na = clustering(prop.nmr.na, method = "kmeans", num.clusters = 4)
kmeans.plot(prop.nmr.na, kmeans.prop.nmr.na)
kmeans.df = kmeans.result.df(kmeans.prop.nmr.na, 4)
kmeans.df## cluster
## 1 1
## 2 2
## 3 3
## 4 4
## samples
## 1 XX_sm
## 2 BR_au CE_au CN_au IT_au SJC_au XX_au IT_sm PU_sm BR_sp CE_sp CN_sp IT_sp PU_sp SJC_sp BR_wi CE_wi CN_wi PU_wi SA_wi SJ_wi XX_wi
## 3 DC_au SA_au SJ_au UR_au SJ_sp UR_sp UR_wi
## 4 AC_au AN_au JB_au PU_au VR_au AC_sm AN_sm BR_sm CE_sm CN_sm DC_sm FP_sm JB_sm SA_sm SJC_sm SJ_sm UR_sm VR_sm AC_sp AN_sp DC_sp FP_sp JB_sp SA_sp VR_sp AC_wi AN_wi DC_wi FP_wi JB_wiPCA
Principal components analysis was performed on the data and some plots are shown below:
pca.analysis.result = pca.analysis.dataset(prop.nmr.na)
pca.pairs.plot(prop.nmr.na, pca.analysis.result, "seasons")
pca.screeplot(pca.analysis.result)
pca.scoresplot2D(prop.nmr.na, pca.analysis.result, "seasons", ellipses = T)
pca.kmeans.plot2D(prop.nmr.na, pca.analysis.result, kmeans.result = kmeans.prop.nmr.na, 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.na, 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.7347738 0.6395083 0.72750 0.9107917 0.9059792 0.1806787
## J48 0.7335476 0.6363272 0.71625 0.9108750 0.8588854 0.1414981
## JRip 0.5489762 0.3951371 0.54625 0.8504167 0.7212812 0.1944324
## svmLinear 0.7388333 0.6443269 0.74125 0.9125000 0.8907500 0.1491218
## rf 0.8304405 0.7691382 0.82000 0.9436250 0.9428437 0.1369526
## KappaSD SensitivitySD SpecificitySD ROCSD
## pls 0.2452844 0.1883441 0.06150896 0.10455008
## J48 0.1920518 0.1526853 0.04795150 0.12101891
## JRip 0.2541152 0.1902409 0.06456466 0.13385105
## svmLinear 0.2051869 0.1630817 0.05160513 0.11458479
## rf 0.1850074 0.1541923 0.04526335 0.08480958Also 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 17.3 0.0 0.6 2.3
## sm 3.4 25.0 2.4 2.3
## sp 3.1 2.2 18.4 4.8
## wi 1.5 0.0 3.8 12.7
##
##
## $J48
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction au sm sp wi
## au 18.7 0.0 2.0 9.8
## sm 0.0 25.3 0.0 0.0
## sp 1.1 1.7 21.1 4.1
## wi 5.6 0.0 2.4 8.2
##
##
## $JRip
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction au sm sp wi
## au 11.0 3.3 2.8 3.1
## sm 5.5 19.9 1.8 3.6
## sp 4.6 1.7 17.7 9.2
## wi 4.1 2.3 3.1 6.3
##
##
## $svmLinear
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction au sm sp wi
## au 19.7 0.0 5.0 3.1
## sm 1.1 25.1 2.2 0.6
## sp 2.6 2.1 14.3 3.5
## wi 2.1 0.0 3.9 14.8
##
##
## $rf
## Cross-Validated (10 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentages of table totals)
##
## Reference
## Prediction au sm sp wi
## au 19.5 0.0 1.6 2.2
## sm 0.2 25.4 0.0 0.0
## sp 2.7 1.7 21.9 3.6
## wi 3.1 0.0 1.7 16.2pls.model = ml.prop.nmr$final.models$pls
pca.plot.3d(prop.nmr.na, 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
## 1.15 75.01267 27.27721 100.00000 64.81564 66.77638
## 3.6 47.65108 31.44005 90.73953 62.62301 58.11342
## 1.27 70.45369 68.69512 43.31358 41.52516 55.99689
## 3.84 67.68811 57.47973 42.94971 35.75304 50.96765
## 1.71 87.56977 42.85494 39.03791 29.87543 49.83451
## 0.99 63.07606 34.95108 55.59388 45.34850 49.74238
## 4.84 26.95760 94.15108 24.07930 46.27577 47.86593
## 3.78 56.29958 38.79599 51.28186 34.37013 45.18689
## 3.87 29.86127 35.38285 58.79089 43.12568 41.79017
## 1.73 59.08808 54.17174 25.06019 27.59580 41.47895
##
## $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.41 100 100
## 1.3 100 100
## 1.36 100 100
## 1.4 100 100
## 1.46 100 100
## 2.14 100 100
## 2.8 100 100
## 3.9 100 100
## 4.02 100 100
## 0.27 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 44.844114 44.844114
## 4.58 24.614271 24.614271
## 6.03 15.948973 15.948973
## 2.32 14.428795 14.428795
## 4.71 12.131944 12.131944
## 6.28 11.629699 11.629699
## 4.17 10.812097 10.812097
## 4.2 10.432548 10.432548
## 4.08 8.678565 8.678565FEATURE 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.na, "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.6328 0.4963 0.1497 0.1972
## 4 0.6933 0.5823 0.1557 0.2010
## 8 0.7248 0.6244 0.1623 0.2123
## 16 0.7447 0.6561 0.1763 0.2316
## 32 0.8037 0.7326 0.1697 0.2305
## 64 0.8181 0.7527 0.1614 0.2174 *
## 242 0.8152 0.7473 0.1736 0.2349
##
## The top 5 variables (out of 64):
## X4.66, X4.5, X4.58, X4.17, X4.38plot(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.na, "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.8227 0.76 0.1505 0.2001
##
## 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 156):
## X0.79 (100%), X1.15 (100%), X1.71 (100%), X3.87 (100%), X3.9 (100%)
##
## On average, 76.4 variables were selected (min = 58, max = 111)