inst/doc/vignette.R

In FAMetA: Fatty Acid Metabolic Analysis

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----echo=T, eval=F-----------------------------------------------------------
# # Install FAMetA
# install.packages("FAMetA", dependencies = c("Depends", "Imports"))
# 
# # load library
# library(FAMetA)
# 

## ----echo=T, eval=F-----------------------------------------------------------
# library(FAMetA)
# library(tools)
# 
# # Example files can be found at
# # <https://drive.google.com/drive/folders/1-mhqBd6W8VJkIYwVuIOr2Gn4Z9vBUN9-?usp=sharing>.
# # This dataset contains 12 samples of linfocites grown in culture media
# # supplemented with 13Cglc, 13Cglc+iFASN, 13Cglc+iSCD and 13Cglc+iFADS2
# # (3 replicates of each condition), 2 blank samples and 3 injections of a
# # standards mix.
# 
# #==============================================================================#
# # Data pre-processing using LipidMS package
# #==============================================================================#
# 
# #=================#
# # Read metadata
# 
# metadata <- read.csv("samples.csv", header = T, sep=",")
# 
# # check file names (they must include .mzXML)
# if (!all(file_ext(metadata$sample) == "mzXML")){
#   metadata$sample[file_ext(metadata$sample) != "mzXML"] <-
#     paste(metadata$sample[file_ext(metadata$sample) != "mzXML"], ".mzXML", sep="")
# }
# 
# 
# #=================#
# # Set processing parameters
# 
# # Peak-picking
# polarity <- "negative"
# dmzagglom <- 15               # dmz and drt to generate bins/partitions for peak-picking
# drtagglom <- 200              # max rt window for bins
# drtclust <- 100               # drt window for clustering (redefines previous bins)
# minpeak <- 8                  # min number of points to define a peak (MS1, MS2)
# minint <- 100000              # at least minpeak points must have minint intensity
# drtgap <- 5                   # max rt gap to fill missing points in a peak
# drtminpeak <- 8               # min width of a peak when there are more than 1 peak in a EIC
# drtmaxpeak <- 30              # max rt window for a EIC
# recurs <- 10                  # max number of peaks in a EIC
# sb <- 5                       # signal-to-baseline ratio (MS1, MS2)
# sn <- 5                       # signal-to-noise ratio
# weight <- 2                   # weight to assign new peaks
# dmzIso <- 5                   # dmz for isotopes search
# drtIso <- 5                   # drt for isotopes search
# 
# parallel <- TRUE              # parallel processing
# ncores <- 4                   # number of cores
# 
# 
# #=================#
# # Peak-picking
# msbatch <- batchdataProcessing(metadata = metadata,
#                                polarity = polarity,
#                                dmzagglom = dmzagglom,
#                                drtagglom = drtagglom,
#                                drtclust = drtclust,
#                                minpeak = minpeak,
#                                drtgap = drtgap,
#                                drtminpeak = drtminpeak,
#                                drtmaxpeak = drtmaxpeak,
#                                recurs = recurs,
#                                sb = sb,
#                                sn = sn,
#                                minint = minint,
#                                weight = weight,
#                                dmzIso = dmzIso,
#                                drtIso = drtIso,
#                                parallel = parallel,
#                                ncores = ncores)

## ----echo=T, eval=F-----------------------------------------------------------
# #=================#
# # Batch processing
# dmzalign <- 10                # max dmz and rt to group peaks for alignment
# drtalign <- 60                # max rt window for clustering in alignment
# span <- 0.2                   # span for alignment
# minsamplesfracalign <- 0.50   # min fraction of samples represented in a peak group to be used for alignment
# dmzgroup <- 10                # max dmz and rt to group peaks for grouping
# drtagglomgroup <- 50          # max rt window for clustering in grouping
# drtgroup <- 10                # max rt difference within a peak group
# minsamplesfracgroup <- 0.20   # min fraction of samples represented in a peak group to be kept
# 
# 
# #=================#
# # Alignment
# msbatch <- alignmsbatch(msbatch, dmz = dmzalign, drt = drtalign, span = span,
#                         minsamplesfrac = minsamplesfracalign,
#                         parallel = parallel, ncores = ncores)
# 
# #=================#
# # Grouping
# msbatch <- groupmsbatch(msbatch, dmz = dmzgroup, drtagglom = drtagglomgroup,
#                         drt = drtgroup, minsamplesfrac = minsamplesfracgroup,
#                         parallel = parallel, ncores = ncores)
# 
# 
# #=================#
# # Save msbatch
# save(msbatch, "msbatch.rda.gz", compress = TRUE)
# 
# 
# # If any other external software is used for processing, data can be loaded from
# # a csv file using the following function:
# # fadata <- readfadatafile("externafadata.csv", sep=",", dec=".")
# 
# # In this case, go directly to data correction step.

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # FA annotation
# #==============================================================================#
# 
# #=================#
# # Annotate FA
# msbatch <- annotateFA(msbatch, dmz = 5)
# 
# #=================#
# # plot peaks from identified FAs to check them
# plots <- plotFA(msbatch, dmz = 10)
# 
# pdf("FAs.pdf")
# for (p in 1:length(plots)){
#   print(plots[[p]])
# }
# dev.off()
# 
# #=================#
# # export annotations for curation
# write.csv(msbatch$fas, file="faid.csv", row.names=FALSE)
# 

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # FA curation
# #==============================================================================#
# 
# #=================#
# # read csv file with modified annotations
# faid <- read.csv("faid_curated.csv", sep=",", dec=".")
# 
# #=================#
# # change FA annotations
# msbatch <- curateFAannotations(msbatch, faid)
# 
# #=================#
# # plot FA peaks again to check identities
# plots <- plotFA(msbatch, dmz = 10)
# 
# pdf("FAs_curated.pdf")
# for (p in 1:length(plots)){
#   print(plots[[p]])
# }
# dev.off()

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # Search FA isotopes and get the fadata object
# #==============================================================================#
# fadata <- searchFAisotopes(msbatch, dmz = 10, coelCutoff = 0.6)
# 
# 
# # if you want to save fadata in a csv to subset it for example:
# df <- cbind(rbind(fadata$fattyacids, data.frame(Compound="IS", Label="")),
#             rbind(fadata$intensities, fadata$IS))
# df <- rbind(c("", "sampletype", fadata$metadata$sampletype),
#             c("Compound", "Label", colnames(fadata$intensities)), df)
# write.table(df, file="fadata.csv", sep=",", col.names = FALSE, row.names = FALSE)
# 
# # and then, you could read it again using:
# fadata <- readfadatafile("fadata.csv", sep=",", dec=".")

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # Import FA data
# #==============================================================================#
# inputfile <- "externalfadata.csv"
# fadata <- readfadatafile(inputfile, sep=",", dec=".")

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # Data correction
# #==============================================================================#
# fadata <- dataCorrection(fadata, blankgroup = "Blank")
# 
# # Alternatively, to add external normalization:
# # fadata <- dataCorrection(fadata, blankgroup = "blank",
# # externalnormalization = "protein")

## ----echo=T, eval=F-----------------------------------------------------------
# #==============================================================================#
# # Metabolic analysis
# #==============================================================================#
# 
# #=================#
# # Synthesis analysis
# fadata <- synthesisAnalysis(fadata=fadata, R2Thr = 0.95, maxiter = 1e3,
#                             maxconvergence = 100, startpoints = 5)
# 
# # If inhibitors have been used, make sure D2 has not been underestimated. If so,
# # D2 could be set as the one calculated for 13-Glc Control samples to improve
# # the results:
# 
# # D2 <- fadata$synthesis$results$D2[fadata$synthesis$results$FA == "FA(16:0)"]
# # fadata$synthesis$results$Group[fadata$synthesis$results$FA == "FA(16:0)"]
# 
# # D2[4:12] <- rep(mean(D2[1:3]))
# 
# # relaunch synthesis analysis fixing D2
# # fadata <- synthesisAnalysis(fadata=fadata, R2Thr = 0.95, maxiter = 1e3,
# #                             maxconvergence = 100, startpoints = 5, D2 = D2)
# 
# 
# # Explore results
# View(fadata$synthesis$results)
# View(fadata$synthesis$predictedValues)
# pdf("plotsSynthesis.pdf")
# for (f in 1:length(fadata$synthesis$plots)){
#   for (s in 1:length(fadata$synthesis$plots[[f]])){
#     print(fadata$synthesis$plots[[f]][[s]])
#   }
# }
# dev.off()
# 
# # to use multinomial distribution without over dispersion, set P parameter to 0
# fadata <- synthesisAnalysis(fadata=fadata, R2Thr = 0.95, maxiter = 1e3,
#                              maxconvergence = 100, startpoints = 5, P = 0)
# 

## ----echo=T, eval=F-----------------------------------------------------------
# #=================#
# # Elongation analysis
# fadata <- elongationAnalysis(fadata, R2Thr = 0.95, maxiter = 1e4,
#                              maxconvergence=100, startpoints = 5, DThr = 0.1)
# 
# 
# # Explore results
# View(fadata$elongation$results)
# View(fadata$elongation$predictedValues)
# pdf("plotsElongation.pdf")
# for (f in 1:length(fadata$elongation$plots)){
#   for (s in 1:length(fadata$elongation$plots[[f]])){
#     print(fadata$elongation$plots[[f]][[s]])
#   }
# }
# dev.off()

## ----echo=T, eval=F-----------------------------------------------------------
# #=================#
# # Desaturation analysis
# fadata <- desaturationAnalysis(fadata)
# 
# # Explore results
# View(fadata$desaturations$results)

## ----echo=T, eval=F-----------------------------------------------------------
# #=================#
# # Summarize results
# fadata <- summarizeResults(fadata, controlgroup = "Control13Cglc")
# 
# #=================#
# # Save fadata
# save(fadata, file="fadata.rda")
# 
# 
# 
# #=================#
# # Export results
# write.csv(fadata$results$results, file = "results.csv", row.names=FALSE)
# write.csv(fadata$results$summary, file = "summary.csv")
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$mid),
#                   fadata$mid),
#             file = "mid.csv", sep=",", col.names = FALSE)
# write.table(rbind(fadata$metadata$sampletype,
#                     colnames(fadata$synthesis$predictedValues),
#                     fadata$synthesis$predictedValues),
#               file = "predictedmid.csv", sep=",", col.names = FALSE)
# 
# 
# pdf("relativepoolsizeRaw.pdf")
# print(fadata$results$heatmaps$relativepoolsize$raw$plot)
# dev.off()
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$heatmaps$relativepoolsize$raw$values),
#                   fadata$results$heatmaps$relativepoolsize$raw$values),
#             file = "relativepoolsizeRaw.csv", sep=",", col.names = FALSE)
# 
# 
# pdf("relativepoolsizeZscore.pdf")
# print(fadata$results$heatmaps$relativepoolsize$zscore$plot)
# dev.off()
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$heatmaps$relativepoolsize$raw$values),
#                   fadata$results$heatmaps$relativepoolsize$zscore$values),
#             file = "relativepoolsizeZscore.csv", sep=",", col.names = FALSE)
# 
# if ("log2FC" %in% names(fadata$results$heatmaps$relativepoolsize)){
#   pdf("relativepoolsizeLog2FC.pdf")
#   print(fadata$results$heatmaps$relativepoolsize$log2FC$plot)
#   dev.off()
# 
#   write.table(rbind(fadata$metadata$sampletype,
#                     colnames(fadata$results$heatmaps$relativepoolsize$log2FC$values),
#                     fadata$results$heatmaps$relativepoolsize$log2FC$values),
#               file = "relativepoolsizeLog2FC.csv", sep=",", col.names = FALSE)
# }
# 
# pdf("resultsRaw_endogenouslySynthesized.pdf")
# print(fadata$results$heatmaps$synthesized$raw$plot)
# dev.off()
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$heatmaps$synthesized$raw$values),
#                   fadata$results$heatmaps$synthesized$raw$values),
#             file = "resultsRaw_endogenouslySynthesized.csv", sep=",",
#             col.names = FALSE)
# 
# 
# if ("log2FC" %in% names(fadata$results$heatmaps$synthesized)){
#   pdf("resultsLog2FC_endogenouslySynthesized.pdf")
#   print(fadata$results$heatmaps$synthesized$log2FC$plot)
#   dev.off()
# 
#   write.table(rbind(fadata$metadata$sampletype,
#                     colnames(fadata$results$heatmaps$synthesized$log2FC$values),
#                     fadata$results$heatmaps$synthesized$log2FC$values),
#               file = "resultsLog2FC_endogenouslySynthesized.csv", sep=",",
#               col.names = FALSE)
# }
# 
# 
# #=================#
# # Isotopologue distributions: observed vs predicted
# 
# pdf("isotopologueDistributions.pdf")
# for (f in 1:length(fadata$synthesis$plots)){
#   for (s in 1:length(fadata$synthesis$plots[[f]])){
#     print(fadata$synthesis$plots[[f]][[s]])
#   }
# }
# for (f in 1:length(fadata$elongation$plots)){
#     for (s in 1:length(fadata$elongation$plots[[f]])){
#       print(fadata$elongation$plots[[f]][[s]])
#     }
#   }
# dev.off()
# 
# 
# 
# #=================#
# # Reorganized tables for synthesis and elongation parameters (S16, E1, E2, E3,
# # E4 and E5)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$S16),
#                   fadata$results$allparameters$S16),
#             file = "S16.csv", sep=",", col.names = FALSE)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$E1),
#                   fadata$results$allparameters$E1),
#             file = "E1.csv", sep=",", col.names = FALSE)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$E2),
#                   fadata$results$allparameters$E2),
#             file = "E2.csv", sep=",", col.names = FALSE)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$E3),
#                   fadata$results$allparameters$E3),
#             file = "E3.csv", sep=",", col.names = FALSE)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$E4),
#                   fadata$results$allparameters$E4),
#             file = "E4.csv", sep=",", col.names = FALSE)
# 
# write.table(rbind(fadata$metadata$sampletype,
#                   colnames(fadata$results$allparameters$E5),
#                   fadata$results$allparameters$E5),
#             file = "E5.csv", sep=",", col.names = FALSE)

## ----echo=T, eval=F-----------------------------------------------------------
# #=================#
# # Customize parameters database
# parameters <- FAMetA::parameters
# 
# # Add a new unknown FA(18:1)
# newrow <- data.frame(FattyAcid = "FA(18:1)nv",
#                      M = 18,
#                      S16 = 1,
#                      E1 = 1,
#                      E2 = 0,
#                      E3 = 0,
#                      E4 = 0,
#                      E5 = 0)
# 
# 
# parameters <- data.frame(rbind(parameters, newrow))
# parameters <- parameters[order(parameters$FattyAcid),]
# View(parameters)
# 
# # Change fatty acid settings: add E1 step for FA(18:3)n6
# parameters$E1[parameters$FattyAcid == "FA(18:3)n6"] <- 1
# 
# 
# # Then add the parameters argument to elongationAnalysis and summarizeResults
# # functions
# fadata <- elongationAnalysis(fadata, R2Thr = 0.95, maxiter = 1e4,
#                              maxconvergence=100, startpoints = 5, D2Thr = 0.1,
#                              parameters = parameters)
# 
# fadata <- FAMetA:::summarizeResults(fadata, controlgroup = "H460_13Cglc",
#                                      parameters = parameters)