intestinimonas.Rmd

---
title: "Physiology of Intestinimonas butyriciproducens"
runtime: shiny
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: https://github.com/microsud/TeachingMIB20306
    theme: lumen
---


```{r setup, include=FALSE}


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

#BiocManager::install("DEP")
#BiocManager::install("SummarizedExperiment")
#BiocManager::install("clusterProfiler")
#BiocManager::install("EnrichmentBrowser")
#BiocManager::install("gage")


#BiocManager::install("ReactomePA")
library(flexdashboard)
library(DEP)
library(SummarizedExperiment)
# install.packages("picante",repos="http://R-Forge.R-project.org") 
# library(picante)
library(clusterProfiler)
library(data.table)
library(DT)
library(RColorBrewer)
library(shiny)
library(shinythemes)
library(pheatmap)
library(plotly)
#library(metagenomeSeq)
library(Biobase)
library("EnrichmentBrowser")
library(plotly)
library(gage)
library(tidyr)
library(dplyr)
library(png)
library(knitr)
library(enrichplot)
library("EnhancedVolcano")
library(reshape2)
library("pathview")
#library(ReactomePA)
library("filesstrings")
#library("shinycssloaders")
options(shiny.maxRequestSize = 30*1024^2) # 30mb limit

av = function(x){
  if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){
    return(NULL)
  }
  return(x)
}

source("codes/make_data_se_file.R")
source("codes/plot_volcano_custom.R")
```

DATA ANALYSIS  
=======================================================================
Sidebar {.sidebar}
-----------------------------------------------------------------------
```{r upload}

# Proteome table
h3("MIB20306 - PRACTICAL 2019")
hr()
img(src="www/WUR_LOGO.png",height=100,width=200)
h6('image source google images')
hr()
h5("This is an interactive Shiny application to visualize and analyze proteomics data")
hr()

tags$hr()

# Example data
actionButton("load_data", label = 'Click here to initiate analysis', width = "225px")

h6('')
tags$strong()

h4('Experiment 3 ')
tags$strong()

h5('Genome-guided elucidation of the physiology of the human gut bacterium Intestinimonas butyriciproducens')
tags$strong()

h6('contact: sudarshan.shetty@wur.nl')
tags$strong()

h5('For educational purpose only')  
tags$strong()

```


```{r reactive}

# Reactive values to store imported data
pseq <- reactive({
  
  # Observe for example data
  if(input$load_data > 0){
    showNotification("Using example dataset...", type = "message")
    data_prot <- readRDS("data/data_se_intestinimonas.rds")
    ps0 <- data_prot
    return(ps0)
  } 
  

  # Err if null
  if(is.null(input$upload_ProteomeData) || is.null(input$upload_SampleData)){
    
    return(NULL)

  }
  

  showNotification("Importing data please be patient...", 
                   duration = NULL, type = "message")
   showNotification("Filtering Reverse != "+", Potential.contaminant != "+" please be patient...", 
                   duration = NULL, type = "message")
    
   
   my_data_se <- make_data_sefile(proteomeData = input$upload_ProteomeData$datapath, 
                              SampleData = input$upload_SampleData$datapath)
  
  return(my_data_se)

})


```


Column {data-width=1000} {.tabset .tabset-fade}
-----------------------------------------------------------------------
### Overview of Experiment 3  

*I. butyriciproducens* cells were inoculated in three anaerobic bottles containing Glucose plus Acetate and in another three anaerobic bottles containing Lysine as sole carbon source. Cells were harvested at late exponential phase, centrifuged to obtain a cell pellet. The cell pellet was subjected to protein extraction procedure.  

![Overview](www/microbiome_sm.jpg)

### Protein Overlap   

This barplot shows the number of different samples in which individual proteins were detected. We can see that most proteins are detected in all six samples.    
 

```{r}

output$LibSizeHist <- renderPlotly({
  
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  
  ps0 <- pseq()
  overlap <- plot_frequency(ps0)
  overlap.fil <- overlap$data
  overlap.fil.sub <- subset(overlap.fil, Var1!=0)
  
  p <- ggplot(overlap.fil.sub, aes(x = Var1, y = Freq, fill = Var1)) + 
  geom_col() + scale_fill_grey(start = 0.8, end = 0.2) + theme_DEP2() + theme(legend.position = "none")
  
  plotly::ggplotly(p+ labs(title = "Protein identifications overlap", 
                           x = "Identified in number of samples", 
                           y = "Number of proteins"), height = 400, width=700)

}
#, width = 400, height = 400
)

fillPage(plotlyOutput("LibSizeHist"))
#fillPage(plotOutput("LibSizeHist"))
```


### Pre-processing and QC  

Quality control (QC) is a critical step for any high-throughput data-driven approach, such as proteomics or transcriptomics. Here, we can filter the data to keep proteins that are detected in any specific  number of replicates per condition (maximum 3). Due to the technical limitations* in measuring, there can be proteins detected in only one of the three replicates, which can bias statistical comparisons between two groups. 
Based on experience, it is normally found that proteins missing in two of the replicates usually has low detection in the third replicate. If the protein is detected in high amounts in one and missing in remaining two replicates, this proteins requires further investigation for potential causes of descrepency.

The two groups in this example represent two different growth conditions; a culture grown with glucose and acetate added, and another culture grown with lysine added. We can change the filter thresholds to test the influence they have on the number of proteins that can be kept and used for further analysis.    


*The missing values can result from very low abundance i.e. below or boderline of the limit of detection for the instrument, misidentification/ambiguous matching where the same peptide matches to multiple proteins sequences and several other reasons related to sampple processing, protein size, location (cell wall or cytoplasmic), etc. [Further reading Lazar et al., J. Proteome Res.2016](https://pubs.acs.org/doi/10.1021/acs.jproteome.5b00981).  

```{r}

fluidPage(
  titlePanel("Filter missing values"),
  
  sidebarLayout(
    
  
sidebarPanel(
helpText("Below set the threshold for allowing 
                          missing values in replicates within the conditions being compared."),
  
  width = 4, column = 3,

 numericInput(
        "Filter_missval_Thres",
        label = "Missing value threshold",0,
       min = 0, 
        max = 10)
),

output$OTUHist_1 <- renderPlot({
  
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  
  ps0 <- pseq()
  data_filt <- filter_missval(ps0, thr = input$Filter_missval_Thres)
  p.num.prot <- plot_numbers(data_filt)
  
  p.num.prot <- p.num.prot + geom_text(aes(label = sprintf("%.1f", sum), y= sum),  
                                     vjust = 3) 
  #p.num.prot <- p.num.prot #+ scale_fill_manual("Carbon Source", values = c("#fbb4ae", "#b3cde3", #ccebc5"))

  #print(p.num.prot)
  return(p.num.prot)
  
  
  }, width = 600, height = 600)


))

``` 

Note 1:
In this step, Filter missing values filters the proteomics dataset based on missing values. The dataset is filtered to keep proteins that have a specified maximum of missing values in at least one condition. In other words, this filter removes proteins that have been missed in any specified number of replicates in at least one condition. This can help make the analysis more stringent by focusing on only those proteins that have been detected with high confidence. In the next tab, we can check effect of this parameter on identifying similarities and differences between replicates and between two conditions.       


### Principal Components Analysis    

The PCA plot is a dimension-reduction method used to identify general properties of data. It is widely used to identify similarities and differences between samples e.g. between replicates or between conditions. How dissimilar are the replicates of the glucose+acetate or the lysine condition? The two axes explain the variance between points.

**Number of top proteins to view:** This sets the number of the most abundant proteins to be used for PCA.   

**Missing value threshold** check the previous section **Pre-processing and QC** for information. Check how changing this value affects replicate as well as between group similarity. Note: Here, y-axis shows variation between replicates and x-axis between two conditions.     
 
```{r}

fluidPage(
  titlePanel("PCA"),
  
  sidebarLayout(
    sidebarPanel(width = 2,  column = 3,
                 ## link the JS file
                 tags$head(tags$script(src="script.js")),
                 ## link the CSS file
                 tags$head(tags$link(rel="stylesheet", 
                                     type="text/css",
                                     href="style.css")),
                 numericInput(
                   "top_prot",
                   label = "Number of most abundant proteins to view",100,
                   min = 10, 
                   max = 100000),
                 helpText("Below set the threshold for allowing 
                          missing values in replicates within the conditions being compared."),
                 numericInput(
                   "Filter_missval",
                   label = "Missing value threshold",0,
                   min = 0, 
                   max = 10)
    ),
    mainPanel(
      plotlyOutput("PCA"))))

output$PCA <- renderPlotly({
  
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  
  ps0 <- pseq()
  data_filt <- filter_missval(ps0, thr = input$Filter_missval)
  data_norm <- normalize_vsn(data_filt)
  set.seed(2156)
  # All possible imputation methods are printed in an error, if an invalid function name is given.
  # Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
  data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)
  
  # Differential enrichment analysis  based on linear models and empherical Bayes statistics
  
  # Test all possible comparisons of samples
  data_diff_all_contrasts <- test_diff(data_imp, type = "all")
  
  dep <- add_rejections(data_diff_all_contrasts, alpha = 0.05, lfc = 1.5)
  
  p.pca <- plot_pca(dep, x = 1, y = 2, n = input$top_prot, point_size = 4,
                    label=F)
  p.pca<- p.pca + theme(legend.title = element_blank()) + ggtitle(paste0("PCA plot - ", input$top_prot, " most abundant proteins"))
  ggplotly(p.pca)
  #return(p.pca)
  #mainPanel(plotlyOutput("p.pca"))
  
}
#, width = 800, height = 600
)


#))
``` 



### Volcano Plot      

Volcano plots are routinely used to visualise results of different -omics experiments. We can visualise statistical significance (P-value) versus magnitude of change (Log2 fold change) in protein expression. The genes or proteins with the largest fold changes may be the most interesting and biologically relevant ones. Again, the reference condition used to define fold-changes is Glucose+Acetate. Consequently, those proteins that have a higher expression in the lysine condition compared to the glucose+acetate condition are shown at the right of the plot, and vice versa. To avoid messy overlap of locus tags, only few of the locus tags with 1.5 fold change between conditions are labelled in the plot below. Complete list can be obtained from the next tab **Differential Expression**.   

**Threshold for the adjusted P value:** This is the threshold for false-discovery rate. It helps prevent false positives from being considered significant. [Further reading](https://link.springer.com/referenceworkentry/10.1007%2F978-1-4419-9863-7_223).   

**Threshold for the log2 fold change:** The log2-fold change is the standard measure of change between conditions. e.g. growth in lysine media compared to growth in glucose media,
Protein A has 1000 and 10 counts in lysin and glucose condition, respectively. Lysine compared to Glucose+Acetate: 1000/10 = 100.log2(100) = 6.6. This means that the expression of Protein A is ca. 6.6-fold higher in the Lysine condition. A positive value indicates that the protein is present at higher abundance in the lysine condition, while a negative value indicates a higher abundance in the glucose plus acetate condition.
Here, Glucose+Acetate is the initial condition. We compare expression in lysine to expression in glucose+acetate. The value here is absolute i.e when 1.5 is set as cut-off, only those prtoeins that have a fold change of -1.5 or below and 1.5 or above are used.  

**Missing value threshold** check the previous section **Pre-processing and QC** for information.    

Note: The cut-offs used here will be the ones used for the table produced in the **Differential Expression tab** later in the workflow. Make note of the values that you define here.  

```{r}

fluidPage(titlePanel("Volcano plot showing differential expression"),
  
  sidebarLayout(
    sidebarPanel(width = 2,  column = 3,

                 numericInput(
                   "alpha",
                   label = "Threshold for the adjusted P value",0,
                   min = 0, 
                   max = 10,
                   step = 0.01),
                 numericInput(
                   "Lfc",
                   label = "Threshold for the log2 fold change",0,
                   min = 0, 
                   max = 50,
                   step= 0.1),
                 
                 helpText("Below set the threshold for allowing 
                          missing values in replicates within the conditions being compared."),
                 numericInput(
                   "Filter_missval2",
                   label = "Missing value threshold",0,
                   min = 0, 
                   max = 6)
                 ),
mainPanel(
            plotOutput('Volcano')
        )))

output$Volcano <- renderPlot({
  
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  
  ps0 <- pseq()
  data_filt <- filter_missval(ps0, thr = input$Filter_missval2)
  data_norm <- normalize_vsn(data_filt)
  set.seed(2156)
  # All possible imputation methods are printed in an error, if an invalid function name is given.
  # Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
  data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)
  
  # Differential enrichment analysis  based on linear models and empherical Bayes statistics
  
  # Test all possible comparisons of samples
  data_diff_all_contrasts <- test_diff(data_imp, type = "all")
  
  dep <- add_rejections(data_diff_all_contrasts, alpha = input$alpha, lfc = input$Lfc)
  
  data_results <- get_results(dep)
  colnames(data_results) <- c("name", "ID", "p.val", "padj", "significant",
                              "Glucose.plus.Acetate_vs_Lysine_significant", 
                              "ratio", "log2FoldChange", "lfc2")
  #p.glu_vs_lysine <- plot_volcano_custom(dep, contrast = "Glucose.plus.Acetate_vs_Lysine", 
  #
  #                                label_size = 2, 
  #                               add_names = TRUE) 
  #return(p.glu_vs_lysine)
  
  
  p.glu_vs_lysine <- EnhancedVolcano::EnhancedVolcano(data_results,
                                                      lab = data_results$name,  
                                                      x = 'log2FoldChange',
                                                      y = 'p.val',
                                                      FCcutoff = input$Lfc,
                                                      cutoffLineType = 'twodash',
                                                      cutoffLineWidth = 0.8,
                                                      transcriptPointSize = 1.5,
                                                      transcriptLabSize = 4.0,
                                                      colAlpha = 1,
                                                      #DrawConnectors= T,
                                                      legend=c('NS','Log (base 2) fold-change','adj.P value',
                                                               'adj.P value & Log (base 2) fold-change'),
                                                      legendPosition = 'right',
                                                      legendLabSize = 12,
                                                      legendIconSize = 5.0)
  
  return(p.glu_vs_lysine)
  
  }, width = 1100, height = 700)


``` 

  

### Differential Expression  

Here, we test the the null hypothesis that there is no difference in the protien expression between glucose plus acetate and lysine condition. 

Note: The exact calulation steps are out of scope for this experiment, if interested students are encouraged to read more about [limma](https://academic.oup.com/nar/article/43/7/e47/2414268). For this Bachelor course, we are interested in understanding the simple biological questions using these standard approaches.  


```{r}


output$Table <- DT::renderDataTable({
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  ps0 <- pseq()
    data_filt <- filter_missval(ps0, thr = input$Filter_missval2)
    data_norm <- normalize_vsn(data_filt)
    set.seed(2156)
  # All possible imputation methods are printed in an error, if an invalid function name is given.
  # Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
    data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)
  
  # Differential enrichment analysis  based on linear models and empherical Bayes statistics
  
  # Test all possible comparisons of samples
  data_diff_all_contrasts <- test_diff(data_imp, type = "all")
  
  dep <- add_rejections(data_diff_all_contrasts, alpha = input$alpha, lfc = input$Lfc)
  data_results <- get_results(dep)
  df_wide <- get_df_wide(dep)
  write.csv(data_results, "data_results.csv")
  data_results.sub <- data_results[,-7]
  colnames(data_results.sub) <-
    c(
    "Locus Tag",
    "Name",
    "Glucose plus Acetate vs Lysine p.val",
    "Glucose plus Acetate vs Lysine p.adj",
    "significant",
    "Glucose.plus Acetate vs Lysine significant",
    "Glucose plus Acetate log2 FC",
    "Lysine log2 FC"
    )
    data_results.sub2 <- data_results.sub
    
    #data_results.sub2$Name <- str_replace(data_results.sub2$Name, "^[>AF_)", "")
    #head(data_results.sub2)
    
    
    data_results.sub2$Name = colsplit(data_results.sub2$Name, " ", names = c("name1", "Name"))[, 2]
    data_results.sub2$Name = colsplit(data_results.sub2$Name, " ", names = c("name1", "Name"))[, 2]
    head(data_results.sub2)
    data_results.sub2$Name <-
    gsub(" reverse ", " ", data_results.sub2$Name)
    data_results.sub2$Name <-
    gsub(" forward ", " ", data_results.sub2$Name)
    return(data_results.sub2)
  
  
},
    server = F, 
    extensions = 'Buttons',
    options = list(lengthChange = TRUE, 
                   dom = 'Bfrtip',
                   buttons = c('copy'),
                   #autoWidth = F,
                   #pageLength = 100,
                   paging = FALSE,
                   searchHighlight = TRUE#,
                   #scrollX = T
                   )
)
      

#fillPage(dataTableOutput('Table'))
fluidPage(dataTableOutput("Table"),style = "height:500px; overflow-y: scroll;overflow-x: scroll;"
          )

```



### Heatmap  

Heatmaps can be useful to investigate, highly variable genes, their abundance, co-expression and similarities between replicates. 

Note: The hierarchy of clusters is depecied by a dendrogram on top and right. On the left proteins are grouped into a cluster with similar abundances and on the top, samples with similar protein expression are grouped together.  

```{r}

fluidPage(titlePanel("Heatmap"),
  
  sidebarLayout(
    
    sidebarPanel(width = 2,  
  
  numericInput(
  "TopVar",
  label = "Most differentially expressed proteins",
  10,
  min = 5,
  max = 1000
  )),
  
mainPanel( 
imageOutput('Heatmap')
        )))


output$Heatmap <- renderPlot({
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
  
  ps0 <- pseq()
  data_filt <- filter_missval(ps0, thr = input$Filter_missval2)
  data_norm <- normalize_vsn(data_filt)
  set.seed(2156)
  # All possible imputation methods are printed in an error, if an invalid function name is given.
  # Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
  data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)
  
  # Differential enrichment analysis  based on linear models and empherical Bayes statistics
  
  # Test all possible comparisons of samples
  data_diff_all_contrasts <- test_diff(data_imp, type = "all")
  
  dep <-
  add_rejections(data_diff_all_contrasts,
  alpha = input$alpha,
  lfc = input$Lfc)
  
  topVarGenes <- head(order(-rowVars(assay(dep))), input$TopVar)
  mat <- assay(dep)[topVarGenes,]
  mat <- mat - rowMeans(mat)
  df <- as.data.frame(colData(dep)[, c("Condition", "ID")])
  #pheatmap(mat, annotation_col = df)
  var_prot <- paste0("Top ", input$TopVar, " differentially expressed proteins")
  heatmap_prot <- pheatmap(mat, annotation_col = df, main = var_prot)
  return(heatmap_prot)
  
  }, width = 900, height = 1200
  )
  
```

### KEGG Pathway IDs-Names  

Use this table to find pathway IDs for metabolic pathways of interest. for. e.g check the KEGG Pathway IDs-Name for Butanoate metabolism. Then use these IDs to search for the corresponding pathways on the next page **KEGG Pathway Analysis**.  

```{r}

choices =   data <- read.csv("data/kegg_pathway_list.csv", 
                   header = T, 
                   row.names = 1)

output$Table2 <- renderDT({
  
  df <- choices},
  options = list(  # options
      scrollX = TRUE, # allow user to scroll wide tables horizontally,
      scrollY = "400px",
      stateSave = FALSE
))

# Reactive values to store imported data

fillPage(dataTableOutput('Table2'))
  
```


### KEGG Pathway Anlaysis  

KEGG pathway analysis is a powerful approach to identify gene/protien expression and thier differences between conditions.   

Note 1: The boxes colored in green are present in the genome. The boxes coloured in red are differentially expresse. Check for thier locus tags and search in **Differential Expression** to see in which condition they are high or low in expression.  

Note 2: Top right of the image is key showing the gradient of color from <b> <font color="red">RED</font></b> to <b><font color="#ffa500">YELLOW</font></b> to <b><font color="steelblue">BLUE</font></b>. 
<b><font color="red">RED</font></b> indicates high expression (fold induction) in Glucose plus acetate condition and <b><font color="steelblue">BLUE</font></b> indicates high expression (fold induction) in Lysine condition. The <b><font color="#ffa500">YELLOW</font></b> color signifies no difference in expression between the two growth conditions.    

For your reference, the KEGG pathway based on genomic data only can be found here: [*Intestinimonas butyriciproducens* AF211](https://www.genome.jp/dbget-bin/www_bget?gn:T04182)  

Please wait for the image to be mapped and viewed in the browser.  


```{r}


choices =   data <- read.csv("data/kegg_pathway_list.csv", 
                   header = T, 
                   row.names = 1)
# List of choices for selectInput
mylist <- rownames(choices)
# Name it
#kk <- enrichKEGG(gene         = gene,
#                 organism     = 'ibu',
#                pvalueCutoff = 0.05)

#browseKEGG(kk, 'ibu00650')
fluidPage(
        headerPanel('Pathway Analysis'),
        sidebarPanel(
            selectInput('Select', label = h3("Pathway"), choices = mylist, selected = NA)
        ),
        mainPanel(
            imageOutput('plot1')
        ),

output$plot1 <- renderImage({
  
  validate(need(expr = !is.null(pseq()), message = "Please upload data..."))  
  
  ps0 <- pseq()
  data_filt <- filter_missval(ps0, thr = input$Filter_missval_Thres)
  data_norm <- normalize_vsn(data_filt)
  set.seed(2156)
  # All possible imputation methods are printed in an error, if an invalid function name is given.
  # Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
  data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)
  
  # Differential enrichment analysis  based on linear models and empherical Bayes statistics
  
  # Test all possible comparisons of samples
  data_diff_all_contrasts <- test_diff(data_imp, type = "all")
  
  dep <- add_rejections(data_diff_all_contrasts, alpha = input$alpha, lfc = input$Lfc)
  data_results <- get_results(dep)
  #df_wide <- get_df_wide(dep)
  
  foldchanges.view = data_results$Lysine_centered
  names(foldchanges.view) = data_results$name
  #gene <- names(foldchanges.view)[abs(foldchanges.view) > 1.5]
  #kk <- enrichKEGG(gene         = gene,
  #                 organism     = 'ibu',
  #                 pvalueCutoff = 0.05,
  #                 pAdjustMethod = "none",
  #                 minGSSize = 5)
  
  id.select <- gsub("ibu","",input$Select) 
  
  outfilepng <- paste0("kegg_maps_out/","ibu", id.select, '.pathview.png')
  
  p <- pathview(gene.data = foldchanges.view, 
                pathway.id = id.select, 
                gene.idtype = "KEGG",
                species = "ibu", 
                #out.suffix="kegg",
                kegg.native=T,
                limit      = list(gene=max(abs(foldchanges.view)), cpd=1),
                same.layer=T,
                low = "#d73027", mid = "#fee090", high = "#4575b4",
                kegg.dir = "kegg_maps")
  
  showNotification("Rendering the pathway map, please wait...", type = "message")
  #browseKEGG(kk, input$Select)
  #pathview::pathview(gene.data = gene, 
   #      pathway.id = input$Select, 
    #     gene.idtype = "REFSEQ",
     #    species = "ibu", 
      #   out.suffix="kegg",
       #           kegg.native=T,
        #          same.layer=T
         #)
  
 file.move(paste0("ibu", id.select, '.pathview.png', sep=''), "kegg_maps_out", 
           overwrite = TRUE )
 
 #outfilepng <- paste0("kegg_maps_out/","ibu", id.select, '.pathview.png')
 
  #pngpath <- paste0("ibu", id.select, '.pathview.png', sep='')
 #pngpath
  list(src =outfilepng)
  
        }, deleteFile = FALSE))


```