Modeling the growth of tomato fruits based on enzyme activity profiles

To assess the influence of the environment on fruit metabolism, tomato (Solanum lycopersicum ‘Moneymaker’) plants were grown under contrasting conditions (optimal for commercial, water limited, or shaded production) and locations. Samples were harvested at nine stages of development, and 36 enzyme activities of central metabolism were measured as well as protein, starch, and major metabolites, such as hexoses, sucrose, organic acids, and amino acids.

Online documentation

Frim1 dataset

ODAM

Metabolism Team



Initialisation

  • Import Functions
source("Functions.R")
  • Init ODAM object : retrieve the dataset structure along with their metadata
# FRIM1 dataset managed by ODAM
dh <- new('odamws', 'https://pmb-bordeaux.fr/getdata/', 'frim1')

# Show dataset structure
options(width=128)
show(dh)
##                        levelName SetID Identifier WSEntry                             Description Count
## 1  plants                            1    PlantID   plant                          Plant features   552
## 2   °--samples                       2   SampleID  sample                         Sample features  1288
## 3       ¦--aliquots                  3  AliquotID aliquot                       Aliquots features   530
## 4       ¦   ¦--cellwall_metabo       4  AliquotID aliquot      Cell wall Compound quantifications    75
## 5       ¦   ¦--cellwall_metaboFW     5  AliquotID aliquot Cell Wall Compound quantifications (FW)    75
## 6       ¦   ¦--activome              6  AliquotID aliquot                       Activome Features   266
## 7       ¦   ¦--plato_hexosesP       10  AliquotID aliquot                       Hexoses Phosphate   266
## 8       ¦   ¦--lipids_AG            11  AliquotID aliquot                               Lipids AG    57
## 9       ¦   °--AminoAcid            12  AliquotID aliquot                             Amino Acids    69
## 10      °--pools                     7     PoolID    pool                Pools of remaining pools   195
## 11          ¦--qMS_metabo            8     PoolID    pool             MS Compounds quantification    25
## 12          °--qNMR_metabo           9     PoolID    pool            NMR Compounds quantification    64


  • Show ‘samples’ factors
dh$getWS('(samples)/factor')
  • Show ‘samples’ quantitative variables
dh$getWS('(samples)/quantitative')



Modelisation of the plant growth

Analysing and modelling plant growth is an important interdisciplinary field of plant science. The use of relative growth rates, involving the analysis of plant growth relative to plant size, has more or less independently emerged in different research groups and at different times and has provided powerful tools for assessing the growth performance and growth efficiency of plants and plant populations.

To measure the bioproductivity of a plant, the component of immediate interest is the net primary production or total yield. The plant weight — usually the dry weight - is a needed measurements for growth analysis. First we have to model the plant growth. The optimization parameters are based on a Levenberg-Marquardt algorithm.

NLS problems with the Levenberg-Marquardt algorithm


Models

  • Models : Weight = f(time-dependent .i.e Day After Anthesis)

Model 1 - Single Sigmoidal Model

  • \(\Large f(t) = d + \frac{a}{1+e^{-b(t-c)}}\)

Model 2 - Sum of two Sigmoidal Model

  • \(\Large f(t) = d + \frac{a}{1+e^{-b(t-c)}} + \frac{f}{1+e^{-g(t-h)}}\)


  • Performation of the modelisation
set.seed(1674)

# Sample subset name
setName <- 'samples'

# Time variable for applying the modelisation
Tname <- 'FruitAge'

# Y variable for applying the modelisation
Yname <- 'FruitDW'

# Additionnal condition
condition <- 'treatment/Control'

# Performation of the modelisation
system.time( fitObj <- fitSigmoid(dh, setName, Tname, Yname, condition, model=1, info=TRUE) )
## Getting ... Fitting ... R2 = 0.9496223  OK
##    user  system elapsed 
##    0.26    0.02    1.76
  • Print the coefficients of the model
print_fittedParams(fitObj)
## Single Sigmoid 
##          a         b        c         d
## 1 6.584323 0.1076615 29.77704 0.1051375
  • Plot both ‘RGR’ & ‘Weigth’ curves
par(mfrow=c(1,2))
plot_fittedCurve(fitObj)
plot_RGRCurve(fitObj)



Explanation of the Relative Growth Rate (RGR)

Relative growth rate is a standardised measure of growth with the benefit of avoiding, as far as possible, the inherent differences in scale between contrasting organisms so that their performances can be compared on an equitable basis
* The Relative Growth Rate (RGR) of a plant at an instant in time (t) is defined as the increase of plant material per unit of material present per unit of time. * The equation is written as: \(\Large RGR = \frac{1}{W}\frac{dW}{dt}\) - W stands for Weight * To calculate the RGR curve, we will leverage on the previous modelisation of the Weight curves (i.e. fitObj object obtained by the fitSigmoid function) * Then we will apply linear modeling using the ‘cv.glmnet’ in order to determine which variable(s) could explain the fruit growth.

  • References
    • Arne Pommerening & Anders Muszta (2015) Methods of modelling relative growth rate, Forest Ecosystems - doi:org/10.1186/s40663-015-0029-4
    • Beauvoit et al (2014) Model-Assisted Analysis of Sugar Metabolism throughout Tomato Fruit Development Reveals Enzyme and Carrier Properties in Relation to Vacuole Expansion, Plant Cell, American Society of Plant Biologists, epub ahead of print. 10.1105/tpc.114.127761 hal-01058814

Linear Modeling with cv.glmnet

run_cvglmnet :Application of linear modeling using the ‘cv.glmnet’ function of the R glmnet package

  • for the Relative Growth Rate (RGR)
  • using the models previously fitted by the fitSigmoid function (fitObj)
  • by selecting some one or more data subsets (setNameList - the data subsets will be merged if severals)
  • alpha=1 => lasso penalty, alpha=0 => ridge penalty, 0 < alpha < 1 => elastic-net penalty

GLMNET


Linear Modeling with cv.glmnet

run_cvglmnet :Application of linear modeling using the ‘cv.glmnet’ function of the R glmnet package

  • for the Relative Growth Rate (RGR)
  • using the models previously fitted by the fitSigmoid function (fitObj)
  • by selecting some one or more data subsets (setNameList - the data subsets will be merged if severals)
  • alpha=1 => lasso penalty, alpha=0 => ridge penalty, 0 < alpha < 1 => elastic-net penalty

GLMNET

# List of data subset names for explaining the RGR
setNameList <- c('activome') # qNMR_metabo

system.time( fitRGR <- run_cvglmnet(fitObj, dh, Tname, setNameList, alpha=1, info=TRUE) )
##  Getting & Merging ... 
## Fitting ... R2 = 0.9994468  OK
##    user  system elapsed 
##    0.47    0.00    2.39
  • CV Glmnet : RGR = F(activome)
print.table( as.matrix(fitRGR$cvfitList), digits=4, zero.print=".", na.print='-' )
##                          1
## (Intercept)      5.730e-01
## PGM             -1.040e-01
## cFBPase          2.321e-03
## PyrK                     .
## CitS                     .
## PFP                      .
## Aconitase                .
## PFK             -6.313e-03
## FruK                     .
## pFBPase                  .
## GluK             1.687e-01
## NAD_ISODH                .
## Enolase          2.971e-01
## NADP_ISODH      -9.027e-05
## PEPC                     .
## Aldolase        -3.589e-02
## Succ_CoA_ligase  2.045e-01
## NAD_MalDH        1.836e-01
## AlaAT            1.827e-03
## Fumarase                 .
## AspAT                    .
## NADP_GluDH               .
## NAD_GAPDH                .
## NADP_GAPDH      -5.692e-03
## NAD_GluDH       -1.480e-01
## TPI                      .
## PGK                      .
## Neutral_Inv     -1.623e-02
## Acid_Inv         3.795e-02
## G6PDH           -2.069e-01
## UGPase                   .
## SuSy                     .
## NAD_ME           1.495e-02
## ShiDH                    .
## NADP_ME                  .
## PGI             -9.186e-02
## StarchS          6.404e-03
## AGPase                   .
## SPS             -8.475e-02
  • ‘RGR’ fitted curve - activome
plot_fittedRGRCurve(fitRGR$cvfitOut)



Plot variables

  • activome
plot_vars(dh, Tname, setNameList[1], smoothtype='lowess', Gmax=20, ncol=4, cex.axis=2, cex.lab=2, cex.main=4)






Annexe: R source code

General

options(width=80)
options(warn=-1)
options(stringsAsFactors=FALSE)

repos='http://cran.rstudio.com'

packages <- c("impute", "pcaMethods","limma")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
   if (paste(R.Version()$major,R.Version()$minor, sep=".") > "3.5") {
      if (!requireNamespace('BiocManager', quietly = TRUE))
         install.packages('BiocManager', repos = repos);
         BiocManager::install(setdiff(packages, rownames(installed.packages())));
   } else {
      source('http://bioconductor.org/biocLite.R');
      biocLite(setdiff(packages, rownames(installed.packages())));
   }
}

packages <- c('Rodam', 'minpack.lm', 'glmnet', 'gplots')
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(packages, rownames(installed.packages())), repos=repos)
}

#-----

library(pcaMethods)
library(gplots)
library(limma)
library(minpack.lm)
library(glmnet)
library(Rodam)
library(UpSetR)

#-----
# Plot functions
#-----

plot.with.errorbars <- function(x, y, err, ylim=NULL, xlab=NULL, ylab=NULL, ...)
{
  if (is.null(ylim))
     ylim <- c(min(y-err), max(y+err))
  plot(x, y, ylim=ylim, pch=19, col='blue', xlab=xlab, ylab=ylab, ...)
  arrows(x, y-err, x, y+err, length=0.05, angle=90, code=3, col='red')
}

plot_vars <- function(dh, Tname, setName, smoothtype="lowess", Gmax=NULL, ncol=4, margin=c(5,6,4,2), ...)
{
   ds1 <- dh$getSubsetByName(setName)
   Vars <- ds1$varnames
   Facs <- ds1$facnames
   data <- ds1$data[, c(Facs, Vars) ]

   X <- data[, Vars ]
   resNIPALS <- pca(as.matrix(X), method = "nipals", center = FALSE)
   X <- resNIPALS@completeObs
   data[ , Vars ] <- X

   t <- as.numeric(gsub("[^0-9\\.]+", "",data[, Tname]))
   M <- as.data.frame(cbind(t, X))
   avg <- aggregate(. ~ t, M, mean)
   sdev <- aggregate(. ~ t, M, sd)

   V <- 1:length(Vars)
   N <- length(Vars)
   i <- 0
   if (is.null(Gmax)) Gmax <- N
   repeat {
      v <- V[(i+1):(i+min(N-i,Gmax))]
      n <- length(v)
      l2 <- ncol
      l1 <- round(n/l2 + 0.4)

    # Plot the variables curves with error bars
      par(mfrow=c(1,1))
      layout(matrix(c(1:(l1*l2)), l1, l2, byrow = TRUE))
      par(mar=margin)
      varDesc <- ds1$LABELS[ds1$LABELS$Attribute %in% Vars,]$Description
      Units <- gsub( ")", "", simplify2array(strsplit(varDesc, "\\(" ))[2,])
      for( k in v ) {
         plot.with.errorbars(avg[,1], avg[,k+1], sdev[,k+1], main=Vars[k], xlab="DPA", ylab=Units[k], ...)
         if (smoothtype=='spline') {
            m <- spline(avg[,1], avg[,k+1])
            lines(m$x, m$y, col='green', lty=1, lwd=1)
         } else if (smoothtype=='lowess') {
            df=data.frame(x=avg[,1], y=avg[,k+1])
            mod=loess(formula = y~x, data=df)
            yfit=predict(mod, newdata=df$x)
            lines(df$x,yfit, col="green",lty=1, lwd=2)
         } else if (smoothtype!='none') {
            lines(avg[,1], avg[,k+1], col='green', lty=1, lwd=1)
         } else {
            lines(avg[,1], smooth(avg[,k+1], smoothtype), col='green', lty=1, lwd=1)
         }
      }
      i <- i + Gmax
      if (i>=N) break
   }
}


Growth Modelisation

# Model 1 - Single Sigmoidal Model
# f(x) = d + a/(1+exp(-b*(x-c)))
sigmoid1 <- function(x,par)
{
   par[4]+par[1]/(1+exp(-par[2]*(x-par[3])))
}
deriv_sigmoid1 <- function(x,par)
{
   (par[1]*par[2]*exp(-par[2]*(x-par[3])))/(1+exp(-par[2]*(x-par[3])))^2
}

# Model 2 - Sum of two Sigmoidal Model
# f(x) = d + a/(1+exp(-b*(x-c))) + f/(1+exp(-g*(x-h)))
sigmoid2 <- function(x,par)
{
   sigmoid1(x,par[1:4]) + sigmoid1(x,c(par[5:7],0))
}
deriv_sigmoid2 <- function(x,par)
{
   deriv_sigmoid1(x,par[1:3]) + deriv_sigmoid1(x,par[5:7])
}

# Normalisation of Data before modeling
normalizeData <-function(dataInput, dataInputName = NA)
{
   timeData <- dataInput$time
   timeRange <- max(timeData,na.rm = T)
   timeData <- timeData / timeRange

   intensityMin <- min(dataInput$intensity,na.rm = T)
   intensityMax <- max(dataInput$intensity,na.rm = T)
   intensityData <- dataInput$intensity - intensityMin
   intensityRange <- max(intensityData,na.rm = T)
   intensityData <- intensityData / intensityRange

   dataOutput <- data.frame(time = timeData, intensity = intensityData)
   return(list(timeIntensityData = dataOutput,
               dataScalingParameters = c(timeRange = timeRange,
                                          intensityMin = intensityMin,
                                          intensityMax = intensityMax,
                                          intensityRange = intensityRange),
               dataInputName = dataInputName))
}

# Fit the sigmoidal model on dataInput (Model 1 & 2)
SigmoidalFitModel <- function(dataInput, model=1, noptim=100, vmax=80, fsign=1)
{
    # The one-sigmoidal model (Model 1)
    oneSigmoid <- function(time,a,b,c) { a/(1+exp(-b*(time-c))) }

    # The double-sigmoidal model (Model 2)
    sumSigmoid <- function(time,a,b,c,f,g,h) {
                    oneSigmoid(time,a,b,c)+oneSigmoid(time,f,g,h)
    }

    normalizedInput <- normalizeData(dataInput = dataInput, dataInputName = "")
    dataFrameInput <- normalizedInput$timeIntensityData

    if (model==1) {
       lowerBounds <- c(a=0.3, b=0.01, c=0); upperBounds <- c(a=1, b=vmax, c=0.95)
       startList <- function() {
            list(a=rnorm(1,1,0.1), b=rnorm(1,1,0.1), c=rnorm(1,0.5,0.1))}
       formula <- intensity ~ oneSigmoid(time,a,b,c)
    } else {
       lowerBounds <- c(a=0.1, b=0.001, c=0.1, f=fsign*0.1, g=0.001, h=0.2)
       upperBounds <- c(a=10, b=vmax, c=1, f=fsign*1, g=vmax, h=1)
       startList <- function() {
            list(a=rnorm(1,0.25,0.05), b=rnorm(1,0.1,0.01), c=rnorm(1,0.25,0.05),
            f=fsign*rnorm(1,0.25,0.05), g=rnorm(1,0.1,0.05), h=rnorm(1,0.5,0.05)) }
       formula <- intensity ~ sumSigmoid(time,a,b,c,f,g,h)
    }
    fitCurve <- list()
    for( idx in 1:noptim ) {
       repeat {
          theFitResult <- try(minpack.lm::nlsLM(formula, dataFrameInput, start = startList(),
                              control = list(maxiter = 1000, minFactor = 1/2^20),
                              lower = lowerBounds, upper = upperBounds, trace=F),
                           silent = TRUE)
          if (class(theFitResult)=="nls") break
       }
       m<-as.data.frame(t(theFitResult$m$getPars()))
       o <- normalizedInput$dataScalingParameters
       if (model==1) {
          p <- c( m$a*o[4], m$b/o[1], m$c*o[1], o[2] )
          R2 <- cor(dataFrameInput$intensity, oneSigmoid(dataFrameInput$time, m$a, m$b, m$c))
       } else {
          p <- c( m$a*o[4], m$b/o[1], m$c*o[1], o[2], m$f*o[4], m$g/o[1], m$h*o[1] )
          R2 <- cor(dataFrameInput$intensity, sumSigmoid(dataFrameInput$time,
                                                         m$a, m$b, m$c, m$f, m$g, m$h))
       }
       fitCurve[[idx]] <- list(params=p, R2=R2)
    }
    R2.best <- 0
    idx.best <- 1
    for( idx in 1:noptim ) {
        R2.idx <- fitCurve[[idx]]$R2
        if (R2.idx>R2.best) { R2.best=R2.idx; idx.best=idx; }
    }
    fitCurve[[idx.best]]
}

# Fit the Growth curve
fitSigmoid <- function(dh, setName, Tname, Yname, model=1, noptim=100, info=FALSE)
{
  if (info) cat("Getting ... ")
     options(stringsAsFactors=FALSE)
     ds1 <- dh$getSubsetByName(setName)
     M <- ds1$data[ , c(Tname, Yname) ]
     M <- M[ ! is.na(M[Yname]), ]
     M[,Yname] <- as.numeric(M[,Yname])
     formula <- as.formula(paste(". ~",Tname))
     AVG <- aggregate(formula, M, mean)
     SDEV <- aggregate(formula, M, sd)
     AVG[,1] <- as.numeric(gsub("[^0-9\\.]+", "",AVG[,1]))
     dat <- data.frame(x=AVG[,1], y=AVG[,2], sdev=SDEV[,2])

  #-----

  if (info) cat("Fitting ... ")
     dataInput <- data.frame(time=dat$x, intensity=dat$y)
     if ( model==1 ) {
        fitObj <- SigmoidalFitModel(dataInput, model=model, noptim=noptim, vmax=50)
        fsig <- sigmoid1
        dfsig <- deriv_sigmoid1
     }
     if ( model==2 ) {
        fsign <- 1
        vmax <-  50
        fitObj <- SigmoidalFitModel(dataInput, model=model, noptim=noptim, fsign=fsign, vmax=vmax)
        fsig <- sigmoid2
        dfsig <- deriv_sigmoid2
     }

     par1 <- as.numeric(fitObj$params)
     R2 <- fitObj$R2
     yest <- fsig(dat$x,par1)
     dyest <- dfsig(dat$x,par1)
     RGR <- dyest/yest
     dat <- data.frame(x=dat$x, y=dat$y, sdev=dat$sdev,
                       ymodel=yest, dymodel=dyest, RGR=RGR)
  if (info) cat("R2 =",R2," OK\n")
     list(data=dat, params=par1, model=model, R2=R2, fsig=fsig, dfsig=dfsig)
}

plot_fittedCurve <- function(fitObj, title="fittedCurve")
{
   dat    <- fitObj$data
   params <- fitObj$params
   fsig   <- fitObj$fsig
   m <- spline(dat$x, dat$y)
   plot.with.errorbars(dat$x, dat$y, dat$sdev, main=title)
   lines(m$x, fsig(m$x,params), col='magenta')
}

plot_RGRCurve <- function(fitObj, title="RGR")
{
   dat   <- fitObj$data
   fsig  <- fitObj$fsig
   dfsig <- fitObj$dfsig
   par1  <- fitObj$params
   x <- dat$x # 1:max(dat$x)
   RGR <- dfsig(x,par1)/fsig(x,par1)
   m <- spline(x, RGR)
   x <- m$x
   RGR <- m$y
   plot( x, RGR/max(RGR), type="l", col="red", main=title)
}

print_fittedParams <- function(fitObj)
{
   par1  <- fitObj$params
   model <- fitObj$model
   if (model==1) {
       names(par1) <- c('a','b','c','d')
       modlabel <- 'Single Sigmoid'
   }
   if (model==2) {
       modlabel <- 'Sum of two Sigmoid'
       names(par1) <- c('a','b','c','d','f','g','h')
   }
   params <- as.data.frame(t(par1))
   cat(modlabel,"\n")
   print(params)
   cat("\n")
}


Linear modeling for Relative Growth Rate (RGR)

# Linear Modeling with cv.glmnet
# TODO : Need to add a validation test - cf https://cran.r-project.org/web/packages/hdi/hdi.pdf
run_cvglmnet <- function(fitObj, dh, Tname, setNameList, alpha=1, info=FALSE)
{
   if (info) cat(" Getting & Merging ... \n")
      options(stringsAsFactors=FALSE)
      set.seed(1674)

    # Get X variables by merging the data subsets (setNameList)
      dsMerged <- dh$getSubsetByName(setNameList)
      Nvars <- length(dsMerged$varnames)
      X <- dsMerged$data[, dsMerged$varnames ]

    # NA imputation
      resNIPALS <- pcaMethods::pca(as.matrix(X), method = "nipals", center = FALSE)
      X <- resNIPALS@completeObs

    # Aggregation according to 'Tname'
      t <- as.numeric(gsub("[^0-9\\.]+", "",dsMerged$data[, Tname]))
      M <- as.data.frame(cbind(t, X))
      M <- aggregate(. ~ t, M, mean)

    # Increase the data size (3 fold) by data interpolating using the cubic spline method
      X.new <- NULL
      for( k in 1:Nvars ) {
          m <- spline(M[,1], M[,k+1])
          X.new <- cbind(X.new, m$y)
          t.new <- m$x
      }
      colnames(X.new) <- dsMerged$varnames
      X <- X.new

    # Scaling of the variables
      X <- scale(X, center=T, scale=T)

    # Get the sigmoid model
      dat   <- fitObj$data
      par1  <- fitObj$params
      fsig  <- fitObj$fsig
      dfsig <- fitObj$dfsig

    # Calculate the RGR corresponding to the same points as for X
      Ymodel <- fsig(t.new,par1)
      dYmodel <- dfsig(t.new,par1)
      RGR <- dYmodel/Ymodel
      RGR <- RGR/max(RGR)
      Y <- RGR
      dataList <- list(X=X, Y=Y, Tx3=t.new)

   if (info) cat("Fitting ... ")
      X   <- dataList$X
      Y   <- dataList$Y
      Tx3 <- dataList$Tx3
      Vars <- colnames(X)

    # Application of linear modeling using the 'cv.glmnet' function of the R glmnet package
      nfolds <- min(20, length(Vars))
      cvfit = glmnet::cv.glmnet(X, Y, type.measure = "mse", nfolds = nfolds, alpha=alpha)

      M <- coef(cvfit, s = "lambda.min")
      cvfitList <- M
      V <- as.vector(M[rownames(M) %in% Vars, ])
      Yest <- X%*%V + M[1,1]
      corVal <- cor( Y, Yest)
      E <- cbind(Tx3, Y, Yest)
      colnames(E) <- c("T", "RGR", "Estimated")
      cvfitOut <- E

   if (info) cat("R2 =",corVal," OK\n")

   list(cvfitList=cvfitList, cvfitOut=cvfitOut, corVal=corVal, alpha=alpha)

}

# Plot the fitted 'RGR' curve
plot_fittedRGRCurve <- function(cvfitOut, title="Fitted RGR") {
   plot(cvfitOut[,1],cvfitOut[,2], col="blue", main=title)
   lines(cvfitOut[,1], cvfitOut[,3],  col='red')
}



Session Information

  • Collect Information About The Current R Session
options(width=80)
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
## [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
## [5] LC_TIME=French_France.1252    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] UpSetR_1.3.3        Rodam_0.1.6         RCurl_1.95-4.11    
##  [4] bitops_1.0-6        glmnet_2.0-16       foreach_1.4.4      
##  [7] Matrix_1.2-14       minpack.lm_1.2-1    limma_3.36.3       
## [10] gplots_3.0.1        pcaMethods_1.72.0   Biobase_2.40.0     
## [13] BiocGenerics_0.26.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3             lattice_0.20-35        tidyr_0.8.3           
##  [4] visNetwork_2.0.4       gtools_3.8.1           assertthat_0.2.1      
##  [7] digest_0.6.22          R6_2.3.0               plyr_1.8.4            
## [10] evaluate_0.14          ggplot2_3.2.1          pillar_1.4.2          
## [13] rlang_0.4.1            lazyeval_0.2.1         rstudioapi_0.10       
## [16] gdata_2.18.0           rmarkdown_1.11         DiagrammeR_1.0.0      
## [19] downloader_0.4         readr_1.1.1            stringr_1.4.0         
## [22] htmlwidgets_1.5.1.9000 igraph_1.2.2           munsell_0.5.0         
## [25] compiler_3.5.1         influenceR_0.1.0       rgexf_0.15.3          
## [28] xfun_0.3               pkgconfig_2.0.2        htmltools_0.4.0.9000  
## [31] tidyselect_0.2.5       tibble_2.0.1           gridExtra_2.3         
## [34] codetools_0.2-15       XML_3.98-1.16          viridisLite_0.3.0     
## [37] crayon_1.3.4           dplyr_0.8.0.1          grid_3.5.1            
## [40] jsonlite_1.6           gtable_0.2.0           magrittr_1.5          
## [43] scales_1.0.0           KernSmooth_2.23-15     stringi_1.4.3         
## [46] viridis_0.5.1          brew_1.0-6             data.tree_0.7.8       
## [49] RColorBrewer_1.1-2     iterators_1.0.10       tools_3.5.1           
## [52] glue_1.3.1             purrr_0.2.5            hms_0.4.2             
## [55] Rook_1.1-1             yaml_2.2.0             colorspace_1.3-2      
## [58] caTools_1.17.1.1       knitr_1.21



About this document

This document was automatically generated within R Studio from a R markdown file (Rmd) helped with the R packages knitr & rmarkdown and the pandoc tool (‘knit to HTML’ in RStudio)

"C:/Program Files/RStudio/bin/pandoc/pandoc" +RTS -K512m -RTS Report1.utf8.md 
          --to html4 
          --from markdown+autolink_bare_uris+ascii_identifiers+tex_math_single_backslash+smart 
          --output Report1.html 
          --email-obfuscation none 
          --self-contained 
          --standalone 
          --section-divs 
          --table-of-contents 
          --toc-depth 4 
          --template "C:\Users\GAIA\Documents\R\win-library\3.5\rmarkdown\rmd\h\default.html" 
          --no-highlight 
          --variable highlightjs=1 
          --variable "theme:bootstrap" 
          --include-in-header "C:\Users\GAIA\AppData\Local\Temp\Rtmpiud8P8\rmarkdown-str4e5c5b7a59f5.html" 
          --mathjax 
          --variable "mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"