Untargeted Metabolomics

#Library

library(timeSeries)

## Cargando paquete requerido: timeDate

## 
## Adjuntando el paquete: 'timeSeries'

## The following objects are masked from 'package:graphics':
## 
##     lines, points

library(reshape2)
library(cluster)
library(gplots)

## 
## Adjuntando el paquete: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library(ggfortify)

## Cargando paquete requerido: ggplot2

library(RColorBrewer)
library(ggplot2)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks timeSeries::filter(), stats::filter()
## ✖ dplyr::lag()    masks timeSeries::lag(), stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(base)
library(graphics)
library(stats)
library(fBasics)
library(pvclust)
#library(pca3d)
library(randomForest)

## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
## 
## Adjuntando el paquete: 'randomForest'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
## 
## The following object is masked from 'package:timeSeries':
## 
##     outlier

library(varSelRF)

## Cargando paquete requerido: parallel

library(e1071)

## 
## Adjuntando el paquete: 'e1071'
## 
## The following objects are masked from 'package:fBasics':
## 
##     kurtosis, skewness
## 
## The following objects are masked from 'package:timeDate':
## 
##     kurtosis, skewness

library(VennDiagram)

## Cargando paquete requerido: grid
## Cargando paquete requerido: futile.logger

library(googlesheets4)

Call the data

normalized <- read.csv(file = "data_normalized.csv")
  
pca.data <- t(normalized) %>% 
  as.data.frame()


colnames(pca.data) <- as.character(unlist(pca.data[1, ]))


pca.data <- pca.data[-1, ]

df.autoscale <- pca.data[, -1] 

df.autoscale[] <- lapply(df.autoscale, function(x) as.numeric(as.character(x)))
 

pca.data <- cbind(row_names = rownames(pca.data), pca.data)

Non-Supervised

autoplot(prcomp(df.autoscale), data = pca.data, colour = 'Label', frame = TRUE, frame.type = 'norm')

#Hierarchical Cluster Analysis

metaData.1 <- pca.data %>%
unite("Sample_Label", row_names, Label, sep = "_")
HCA.data <- df.autoscale
rownames(HCA.data) <- metaData.1[,1]
distEucl<-dist(HCA.data)
hclust(distEucl)

## 
## Call:
## hclust(d = distEucl)
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 82

distcorr<-as.dist(1-cor(t(HCA.data)))
hclust(distcorr)

## 
## Call:
## hclust(d = distcorr)
## 
## Cluster method   : complete 
## Number of objects: 82

plot(hclust(distEucl)) 

plot(hclust(distcorr))

#Run just one time
if (!requireNamespace("BiocManager", quietly = TRUE))
 install.packages("BiocManager")
BiocManager::install("ropls")

## Bioconductor version 3.19 (BiocManager 1.30.25), R 4.4.1 (2024-06-14 ucrt)

## Warning: package(s) not installed when version(s) same as or greater than current; use
##   `force = TRUE` to re-install: 'ropls'

## Installation paths not writeable, unable to update packages
##   path: C:/Program Files/R/R-4.4.1/library
##   packages:
##     boot, foreign, MASS, nlme, survival

## Old packages: 'bitops', 'digest', 'glue', 'timeDate', 'timeSeries', 'xfun'

library(ropls)

## 
## Adjuntando el paquete: 'ropls'

## The following object is masked from 'package:tibble':
## 
##     view

opls.data <- cbind(pca.data[,1:2], df.autoscale)

oplsda.T5_T6T8 <- opls(opls.data[,-c(1:2)], opls.data[, "Label"], algoC = "default", scaleC = "none",predI = 1)

## PLS-DA
## 82 samples x 1188 variables and 1 response
## none scaling of predictors and standard scaling of response(s)
##       R2X(cum) R2Y(cum) Q2(cum) RMSEE pre ort pR2Y  pQ2
## Total     0.45    0.541   0.528 0.325   1   0 0.05 0.05

## Warning: Single component model: only 'overview' and 'permutation' (in case of
## single response (O)PLS(-DA)) plots available

getSummaryDF(oplsda.T5_T6T8)

##       R2X(cum) R2Y(cum) Q2(cum) RMSEE pre ort pR2Y  pQ2
## Total     0.45    0.541   0.528 0.325   1   0 0.05 0.05

VIP.T5_T6T8 <- getVipVn(oplsda.T5_T6T8)
VIPT5_T6T8 <-head(sort(VIP.T5_T6T8, decreasing = TRUE), 20)
write.csv(VIPT5_T6T8, file = "VIPT5_T6T8.csv")