.. _`ris-r-workshop`:

==============
RIS R Workshop
==============

.. contents::
   :depth: 1
   :local:

.. Note::

    Connecting to get command line access: ``ssh wustlkey@compute1-client-1.ris.wustl.edu``

    Queue to use: ``workshop, workshop-interactive``

    Group to use: ``compute-workshop`` (if part of multiple groups)


.. _`ris-r-workshop-video`:

R Workshop Video
----------------

.. raw:: confluence

    <ac:structured-macro ac:name="html">
        <ac:plain-text-body><![CDATA[

.. raw:: confluence html

    <iframe width="560" height="315" src="https://www.youtube.com/embed/z6upbR0cvZA" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

.. raw:: confluence

    ]]></ac:plain-text-body>
        </ac:structured-macro>

.. _`ris-r-workshop-provides`:

What will this documentation provide?
-------------------------------------

- An introduction to using R on the RIS Compute Platform.
- Expanding the base r-base Docker image to include new packages.
- Use case example of data analysis with R within the RIS Compute Platform.

.. _`ris-r-workshop-needed`:

What is needed?
---------------

- Access to the RIS Compute Platform.
- The example data. Found here: https://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1000/matrix/
- An account on Docker Hub.

.. _`ris-r-workshop-images`:

Available R Docker Images
-------------------------

- There are already basic R Docker images available from public sources. Information about what sources RIS recommends, :ref:`can be found here. <ris-vetted-app-containers>`
- This workshop will be using the r-base publicly available official image as a base for building upon. https://hub.docker.com/_/r-base

.. _`ris-r-workshop-start`:

Getting Started
---------------

- You can run the latest r-base R Docker image on the RIS Compute Platform with the following command.

.. code::

    export LSF_DOCKER_VOLUMES='/path/to/storage:/path/to/storage /path/to/scratch:/path/to/scratch'

    bsub -Is -G compute-workshop -q workshop-interactive -a 'docker(r-base:latest)' /bin/bash


- This example shows running an R script that prints out X terms of the fibonacci sequence where X is an integer input by the user in the script.

.. image:: images/compute104/R.example.job.start.gif

.. image:: images/compute104/R.example.script.gif

.. image:: images/compute104/R.example.run.gif

.. image:: images/compute104/R.example.run2.gif

.. image:: images/compute104/R.example.run3.gif

- This type of job can also be ran as a non interactive job by changing up the bsub command slightly.

.. code::

    export LSF_DOCKER_VOLUMES='/path/to/storage:/path/to/storage /path/to/scratch:/path/to/scratch'

    bsub -G compute-workshop -q workshop -a 'docker(r-base:latest)' Rscript /path/to/script

.. image:: images/compute104/R.example.run.ni.gif

.. image:: images/compute104/R.example.email.png

.. _`ris-r-workshop-packages`:

Adding Packages to the Docker Image
-----------------------------------

.. Note::

    - The example of Docker image building contained here assumes we are building on a local machine.
    - You can build Docker images using the Compute Platform as well. :ref:`That documentation can be found here. <build-image-on-compute>`

- The first thing we need to do is determine what version of R we what to use and use that specific tag in the base docker build.
- In this process, we will be creating a Dockerfile, which can then be used by Docker to build a Docker image.
- In this case we're going to use ``R 4.0.3``.

.. code::

    #Set base image from r-base using R version 4.0.3
    FROM r-base:4.0.3

- Once we have the base image decided, we need to tell the Docker file to install the packages we want.
- In the case of our example we'll be using DeSeq2.

.. code::

    #Set base image from r-base using R version 4.0.3
    FROM r-base:4.0.3

    #Install needed R packages
    RUN R -e "install.packages('DeSeq2',dependencies=TRUE, repos='http://cran.rstudio.com/')"

.. Note::

    - You can install multiple packages at the same time by adjusting the install command to the following.

    .. code::

        RUN R -e "install.packages(c('EIAdata', 'gdata', 'ggmap', 'ggplot2'), dependencies=TRUE, repos='http://cran.studio.com/')"

- Once we save the file, we then make sure we are in the working directory that contains the Docker directory.
- Then, we can tell Docker to build this image using the following command.

.. code::

    docker build -t username/Dockerfile:tag directory-of-docker

.. Note::

    - ``username`` is the username of your Docker hub account.
    - ``Dockerfile`` is the name you wish to give the Docker image
    - ``tag`` is the tag you want to give this particular image, you can use latest or something like ``1.0``
        - If no tag is given, Docker will default to the ``latest`` tag, overwriting the previous image with the same name and the ``latest`` tag.
    - ``directory-of-docker`` is the directory where you have saved the particular Dockerfile we are using to build the image.

- When we try to build the image using ``R 4.0.3``, we get an error for DeSeq2.

.. image:: images/r-workshop/deseq2.error.png

- This error is because ``DESeq2`` is not available on ``CRAN``, the default R packages location. It is available via Bioconductor. http://bioconductor.org/packages/release/bioc/html/DESeq2.html

.. Note::

    - There are a lot of bioinformatics packages that are not available through ``CRAN`` but are available through Bioconductor.

- To install ``DESeq2`` we first need to install the Bioconductor manager package as this is required for Bioconductor packages.
- To do that, our updated Dockerfile code would look like the following.

.. code::

    #Set base image from r-base using R version 4.0.3
    FROM r-base:4.0.3

    #Install needed R packages
    RUN R -e "install.packages('BiocManager',dependencies=TRUE, repos='http://cran.rstudio.com/')"
    RUN R -e "BiocManager::install('DESeq2')"

- Now when we run the ``docker build`` command, it will attempt to install the packages we've told it to.
- In this case that is ``BiocManager`` and all of it's dependencies, followed by ``DESeq2``.
- As we look at the R output the Docker gives us when attempting to build this package, we can see that an error occurs.

.. image:: images/r-workshop/biocmanager.error.png

- If we look further back in the install output, we can see where this error originates.

.. image:: images/r-workshop/curl.error.png

- When the image tries to install ``DESeq2`` we get a few errors.

.. image:: images/r-workshop/deseq2.install.errors.png

- A lot of these errors look like they're dependency errors.
- By looking at the ``BiocManager`` error, we can see that we need to install a base OS library  that's needed for the R package.
- We can do that by adding the following code to install said library through the appropriate install method for the OS.
- The r-base Docker images use Debian as their OS, so we can install the ``libcurl4-openssl-dev`` library with ``apt-get``.

.. Note::

    - You can find the OS version of Linux by using the ``cat /etc/os-release`` command.

.. code::

    #Set base image from r-base using R version 4.0.3
    FROM r-base:4.0.3

    #Install library dependencies
    RUN apt-get update \
        && apt-get install -y --no-install-recommends libcurl4-openssl-dev \
        && apt-get clean

    #Install needed R packages
    RUN R -e "install.packages('BiocManager',dependencies=TRUE, repos='http://cran.rstudio.com/')"
    RUN R -e "BiocManager::install('DESeq2')"

- We can rerun the ``docker build`` command with the library installed and see if that solves our issues.
- This time we will use the ``--no-cache`` option of the ``docker build`` command so that we can see if we have any install issues.

.. code::

    docker build --no-cache -t username/Dockerfile:tag directory-of-docker

.. Note::

    - Using the ``--no-cache`` option means that Docker will not use the cache when building the image, which means that everything will compile fresh.
    - This option is not recommended all the time, only when attempting to debug an image.

- This time when the image tries to install the ``Curl`` package, it is successful.

.. image:: images/r-workshop/curl.success.png

- We can also see the message at the end of the ``BiocManager`` contains no warnings.

.. image:: images/r-workshop/biocmanager.success.png

- However, ``DESeq2`` still contains errors when it tries to install.

.. image:: images/r-workshop/deseq2.install.errors2.png

- If we delve into the install code, we can see that we need two more libraries based on the errors given.

.. image:: images/r-workshop/openssl.install.error.png

.. image:: images/r-workshop/xml.install.error.png

- For the first package, the error gives us the library needed, which is: ``libssl-dev`` for Debian.
- For the second package, we need to to a Google search of ``xml2-config`` as that is what the error says the package needs.
- The search results lead us to the ``libxml2-dev`` library for Debian.
- Now that we have identified which libraries we need, we can add them to the Docker file and run the build command again.
- Our Dockerfile now looks like the following.

.. code::

    #Set base image from r-base using R version 4.0.3
    FROM r-base:4.0.3

    #Install library dependencies
    RUN apt-get update \
        && apt-get install -y --no-install-recommends libcurl4-openssl-dev libssl-dev libxml2-dev \
        && apt-get clean

    #Install needed R packages
    RUN R -e "install.packages('BiocManager',dependencies=TRUE, repos='http://cran.rstudio.com/')"
    RUN R -e "BiocManager::install('DESeq2')"

- After this build attempt, we can see that we no longer have any errors the ``DESeq2`` install either.

.. image:: images/r-workshop/deseq2.success.png

- Now that we have the R package we want to use successfully installed in a Docker image, we can push said image to our Docker Hub.

.. code::

    docker push username/dockerfile:tag

- Once you push the image it will be publicly available on Docker hub, unless you have changed settings.

.. image:: images/r-workshop/docker.push.png

.. image:: images/r-workshop/docker.hub.example.png

.. _`ris-r-workshop-using-image`:

Using the R Docker Image
------------------------

- Now that we have our Docker image upgraded and hosted on Docker hub, we can use it like we did the r-base Docker image.

.. code::

    export LSF_DOCKER_VOLUMES='/path/to/storage:/path/to/storage /path/to/scratch:/path/to/scratch'

    bsub -Is -G compute-workshop -q workshop-interactive -a 'docker(username/dockerfile:tag)' /bin/bash

.. Note::

    - ``username/dockerfile:tag`` is the name we gave the image when we built it and pushed it to Docker hub.

.. image:: images/r-workshop/running.docker.image.png

.. image:: images/r-workshop/deseq2.installed.png

.. _`ris-r-workshop-setup`:

Analysis Setup
--------------

- We can now use ``DESeq2`` to do analysis on the sample data we previously gathered.
- If we look at our data file using the linux ``head`` command, we can see that it has a lot of information at the beginning of the file.
- Each line of information starts with a ``!``. We can remove those lines so that the file plays nice with R.

.. image:: images/r-workshop/example.info.head.png

- We can use the following command to do this.

.. code::

    grep -v "\!" GSE1000_series_matrix.txt > GSE1000_series_matrix_fixed.txt

.. Note::

    - This file should be in the directory you will be using as the working directory for R.
    - This command should be run from within the working directory. If not, you will need the full path to the file and the output.

.. image:: images/r-workshop/example.without.info.png

- We can now see that the file has only a single header and then a row for each sample and columns of expression data.
- The first thing we want to do is set our working directory and to load the data into R.

.. code::

    R> setwd("/path/to/storage/")
    R> data <- read.delim("GSE1000_series_matrix_fixed.txt", header = TRUE, sep = "\t")

.. image:: images/r-workshop/r.data.load.png

- Now that we have the data loaded, we'll want to manipulate the data just a little bit to make it more friendly to use.
- Then we'll set the row names to the Gene IDs and then remove that column.

.. code::

    R> rownames(data) <- data[,1]
    R> data <- data[,-1]

.. image:: images/r-workshop/r.rownames.png

- For analysis with ``DESeq2`` there needs to be phenotype data for the samples. Since we don't have that available, we'll generate some example data.
- In this dataset, there are 10 samples, so we'll create a 10 by 4 matrix.

.. code::

    R> pheno.data <- as.data.frame(matrix(rep(rep(0, 4), 10), nrow=10, ncol=4))

.. image:: images/r-workshop/r.create.matrix.png

- Now that we have the dataframe created, we need to add column and row names.

.. code::

    R> rownames(pheno.data) <- colnames(data)
    R> colnames(pheno.data) <- c("group","treatment","disease","weight")

.. image:: images/r-workshop/r.pheno.names.png

- Now we need to populate the data.
- The group data we will create 2 groups for the group column, the first group will simply be the first 5 samples and the remaining 5 will be the second group.

.. code::

    R> pheno.data[,1] <- c(rep(1,5), rep(2,5))

- Next we will create 2 treatments for the treatment column, we will alternate the treatments bewtween 1 and 2 regardless of group.

.. code::

    R> pheno.data[,2] <- c(rep(c(1,2), 5))

- Next we need to create a binary disease state of either 0 or 1. We will do this randomly.

.. code::

    R> pheno.data[,3] <- c(round(abs(rnorm(10, mean=0.5, sd=0.15)), digits = 0))

- Finally we will need to generate a random set of weights. We can generate those very similarly to the disease state.

.. code::

    R> pheno.data[,4] <- c(round(rnorm(10, mean=175, sd=25), digits = 0))

.. image:: images/r-workshop/r.pheno.data.png

.. Note::

    - You can read more about the functions used to create the example data via the R documentation.
    - ``rep``: https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/rep
    - ``rnorm``: https://www.rdocumentation.org/packages/compositions/versions/2.0-0/topics/rnorm
    - ``abs``: https://www.rdocumentation.org/packages/SparkR/versions/2.1.2/topics/abs
    - ``round``: https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/Round

.. _`ris-r-workshop-deseq2-analysis`:

DESeq2 Analysis
---------------

- We now have our expression data ``data`` and our phenotypic data ``pheno.data``, now we're ready to analyze the data via ``DESeq2``
- The first step is loading the ``DESeq2`` library.

.. code::

    R> library(DESeq2)

.. Note::

    - The manual for DESeq2 can be found here: http://bioconductor.org/packages/release/bioc/manuals/DESeq2/man/DESeq2.pdf

- With the library loaded, we need to create a ``DESeq2`` object that contains the data we've generated.

.. code::

    R> deseq2.data <- DESeqDataSetFromMatrix(data, pheno.data, ~ group + treatment + disease + weight)

.. image:: images/r-workshop/r.integer.error.png

- As we can see from the error we receive when attempting this command, our count data needs to be in integers, ie should not include decimals.
- We can fix this via the following method.

.. code::

    R> data <- round(data, digits = 0)

- When we do and rerun the ``DESeqDataSetFromMatrix`` command we now get a warning about our data and that certain columns of data should be designated as factors.
- We can adjust that using the following method.

.. code::

    R> pheno.data$group <- factor(pheno.data$group, levels = c("1","2"))
    R> pheno.data$treatment <- factor(pheno.data$treatment, levels = c("1","2"))
    R> pheno.data$disease <- factor(pheno.data$disease, levels = c("0","1"))

- Though this does not appear to get rid of the warning, we can continue doing analysis fine.
- The next step is to run the actual analysis using the ``DESeq()`` command.

.. code::

    R> dds <- DESeq(deseq2.data)

.. image:: images/r-workshop/r.analysis.png

- Now we want to put the results into an object so that we can see what they are.

.. code::

    R> res <- results(dds)

- Now if we want to look at a head display of the results we can simply type the variable name ``res`` and it will give us a nice display.

.. image:: images/r-workshop/r.results.png

- We can see here that there are several columns for each gene that was analyzed. The column of interest for us at this point is the p-value.
- We can write a line of code using a ``which`` statement combined with the use of the dataframe signifiers to get a list of the genes that are below a certain p-value.
- First we need an object with the gene names. We can do that by using ``rownames`` and our dataset.

.. code::

    R> genes <- rownames(data)

- Then we create our statement to get which genes are significant at p-value < 0.01.

.. Note::

    - In this case we are selecting p < 0.01 as our signifcant threshold, but other values can be chosen.
    - There are statistical methods to determine what p-value you should use, but those are not a part of this workshop.

.. code::

    R> significant.genes <- genes[which(res$pvalue < 0.01)]

.. image:: images/r-workshop/r.sig.genes.png

- There are a lot of genes here that are significant. We can find out how many by using the ``length`` command.

.. code::

    R> length(significant.genes)

.. image:: images/r-workshop/r.number.genes.png

- As we can see, there are 446 genes that are significant (for this example data) based on the values we gave. We can use that data for displaying this information via graphs and charts.

.. Note::

    - You can find more information about the ``which`` and ``length`` commands at the following links.
        - ``which``: https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/which
        - ``length``: https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/length

.. Note::

    - There will be differences in the number of genes you find significant based upon the phenotype data. Since this data is randomly generated, it will change the results of the output.
    - You can test this yourself by re-running the data generation steps and doing the analysis up to the ``length`` step and seeing the differences.
    - You can do this until you have a workable number (~500) of significant genes to do further analysis with.

.. _`ris-r-workshop-further-analysis`:

Further Analysis
----------------

- With our now identified significant genes, we can do further analysis and create graphs.
- The first thing we'll want to do is create a data object with just the significant genes.

.. code::

    R> sig.data <- data[which(res$pvalue < 0.01),]

- Now we are going to look at clustering the data based off of the significant gene counts to see what that looks like.
- We will do a simple K-means clustering analysis using the ``kmeans`` command.
- But first we must transpose the dataframe for the kmeans analysis to work correctly.

.. code::

    R> t.sig.data <- as.data.frame(t(sig.data))
    R> rownames(t.sig.data) <- colnames(sig.data)
    R> colnames(t.sig.data) <- rownames(sig.data)

- We have two groups and two treatments. The first step will be to cluster into 2 groups.

.. code::

    R> fit2 <- kmeans(t.sig.data, 2)

- We can look at what groups the samples got clustered into by looking at ``fit$cluster``.

.. image:: images/r-workshop/r.fit.cluster.png

- Now we can compare that back to the original groupings by putting this information into a dataframe.

.. code::

    R> compare.data <- data.frame(pheno.data$group, pheno.data$treatment, pheno.data$disease, fit2$cluster)

.. image:: images/r-workshop/r.compare.cluster.png

- As we can see, it does not appear to group the data based off any of our phenotypes.
- We can see if using 4 groups changes that or not.

.. code::

    R> fit4 <- kmeans(t.sig.data, 4)
    R> compare.data2 <- data.frame(pheno.data$group, pheno.data$treatment, pheno.data$disease, fit4$cluster)

.. image:: images/r-workshop/r.compare.cluster2.png

- As we can see here, using the significant genes to group the samples into 4 groups does not explain the data any better than the 2 groups.
- If we want to save the significant genes to a file, we can do that as well use the ``write`` command.

.. code::

    R> write(significant.genes, "significant.genes.txt")

.. _`ris-r-workshop-visualization`:

Visualization
-------------

- You can visualize data using many different methods in R and we'll go through a couple of examples.

.. Note::

    Because of the way jobs run on the server and given the visualization limitations of ssh connections,
    most methods will have to output the visualization directly to a file.

- The first thing we're going to start off with is a simple bar chart. We'll create one based off the first significant gene.

.. code::

    R> barplot(t.sig.data[,1], main="Gene 1", xlab="Counts", names.arg=rownames(t.sig.data))

- This is simply the command to create a bar plot, now we have to make it so that it prints the plot to a file

.. code::

    R> pdf("gene1.pdf")
    R> barplot(t.sig.data[,1], main="Gene 1 Counts", names.arg=rownames(t.sig.data))
    R> dev.off()

- As we can see, it generates a pdf that has our graph.

.. image:: images/r-workshop/r.barplot.png

- If we look at the graph, we can see that some of the names aren't displayed as they're rather long, we can adjust this with a couple of options.
    - The first is to rotate the names with the ``las`` option.
    - The second is to adjust the margin using the parameter command ``par()``.
    - The third is to adjust the height on the pdf. The default is 7.

.. code::

    R> pdf("gene1.pdf", height=10)
    R> par(mar=c(11,4,4,4))
    R> barplot(t.sig.data[,1], main="Gene 1 Counts", names.arg=rownames(t.sig.data), las=2)
    R> dev.off()

.. image:: images/r-workshop/r.barplot.fixed.png

- Given that we can have multiple, and possibly hundreds of genes that we may want to create a graph for, we can use a for loop in R to loop through the genes.
- The first thing we'll need to do is adjust the commands to use the names of the genes for the pdf file and for the graph title.

.. code::

    R> pdf(paste(significant.genes[1], ".pdf", sep=""), height=10)
    R> par(mar=c(11,4,4,4))
    R> barplot(t.sig.data[,1], main=paste(significant.genes[1], "Counts", sep=" "), names.arg=rownames(t.sig.data), las=2)
    R> dev.off()

.. image:: images/r-workshop/r.barplot.with.name.png

- Now that we have that setup, we can introduce the for loop so that it loops through all the significant genes.

.. code::

    R> for (i in 1:length(significant.genes)) {
    R> pdf(paste(significant.genes[i], ".pdf", sep=""), height=10)
    R> par(mar=c(11,4,4,4))
    R> barplot(t.sig.data[,i], main=paste(significant.genes[i], "Counts", sep=" "), names.arg=rownames(t.sig.data), las=2)
    R> dev.off()
    R> }

- When we run this, we get an error because there are special characters in the gene name that Linux does not like in file names.

.. image:: images/r-workshop/r.file.name.error.png

- To correct this, we can use the ``str_replace_all`` command combined with a regular expression to replace all non alphanumberic characters with and underscore.
- This command is part of the ``stringr`` library.

.. code::

    R> library(stringr)
    R> str_replace_all(significant.genes[i], regex("\\W+"), "_")

- Once we add that to our for loop, we see that we no longer get the error.

.. code::

    R> for (i in 1:length(significant.genes)) {
    R> pdf(paste(str_replace_all(significant.genes[i], regex("\\W+"), "_"), ".pdf", sep=""), height=10)
    R> par(mar=c(11,4,4,4))
    R> barplot(t.sig.data[,i], main=paste(significant.genes[i], "Counts", sep=" "), names.arg=rownames(t.sig.data), las=2)
    R> dev.off()
    R> }

.. image:: images/r-workshop/r.str.replace.png

- If we look in the working directory, we can see that we have a list of graphs of all the genes that are significant.

.. image:: images/r-workshop/r.all.graphs.png

- You can put all of these commands into one R script.

.. code::

    #Set the working directory
    setwd("/path/to/working/directory")

    #Load the Libraries
    library(DESeq2)
    library(stringr)

    #Read the Expression Data into an Object
    data <- read.delim("GSE1000_series_matrix_fixed.txt", header = TRUE, sep = "\t")

    #Rename the row names with the Gene IDs and remove the Gene ID row
    rownames(data) <- data[,1]
    data <- data[,-1]
    data <- round(data, digits = 0)

    #Create a data frame for the create phenotype data
    pheno.data <- as.data.frame(matrix(rep(rep(0, 4), 10), nrow=10, ncol=4))

    #Give the row and column names the appropriate sample and phenotype names
    rownames(pheno.data) <- colnames(data)
    colnames(pheno.data) <- c("group","treatment","disease","weight")

    #Populate the phenotype data frame
    pheno.data[,1] <- c(rep(1,5), rep(2,5))
    pheno.data[,2] <- c(rep(c(1,2), 5))
    pheno.data[,3] <- c(round(abs(rnorm(10, mean=0.5, sd=0.15)), digits = 0))
    pheno.data[,4] <- c(round(rnorm(10, mean=175, sd=25), digits = 0))

    #Fixing the phenotype data
    pheno.data$group <- factor(pheno.data$group, levels = c("1","2"))
    pheno.data$treatment <- factor(pheno.data$treatment, levels = c("1","2"))
    pheno.data$disease <- factor(pheno.data$disease, levels = c("0","1"))

    #Creating the DESeq2 object
    deseq2.data <- DESeqDataSetFromMatrix(data, pheno.data, ~ group + treatment + disease + weight)

    #Running the DESeq2 analysis
    dds <- DESeq(deseq2.data)

    #Getting Results
    res <- results(dds)
    genes <- rownames(data)
    significant.genes <- genes[which(res$pvalue < 0.01)]

    #Getting the significant genes
    sig.data <- data[which(res$pvalue < 0.01),]

    #Doing K-means clustering analysis
    t.sig.data <- as.data.frame(t(sig.data))
    rownames(t.sig.data) <- colnames(sig.data)
    colnames(t.sig.data) <- rownames(sig.data)
    fit2 <- kmeans(t.sig.data, 2)
    compare.data <- data.frame(pheno.data$group, pheno.data$treatment, pheno.data$disease, fit2$cluster)
    fit4 <- kmeans(t.sig.data, 4)
    compare.data2 <- data.frame(pheno.data$group, pheno.data$treatment, pheno.data$disease, fit4$cluster)

    #Writing the list of significant genes to a file
    write(significant.genes, "significant.genes.txt")

    #Visualization of data

    #Bar Plot
    pdf("gene1.pdf", height=10)
    par(mar=c(11,4,4,4))
    barplot(t.sig.data[,1], main="Gene 1 Counts", names.arg=rownames(t.sig.data), las=2)
    dev.off()

    for (i in 1:length(significant.genes)) {
    pdf(paste(str_replace_all(significant.genes[i], regex("\\W+"), "_"), ".pdf", sep=""), height=10)
    par(mar=c(11,4,4,4))
    barplot(t.sig.data[,i], main=paste(significant.genes[i], "Counts", sep=" "), names.arg=rownames(t.sig.data), las=2)
    dev.off()
    }

.. _`ris-r-workshop-rscript`:

Running an R-script
-------------------

- If we put the previous into a file and call it example.R, we have created an R-script that we can run without having to go into R.
- The R-script can be run either interactively or a non-interactive job.

.. code::

    export LSF_DOCKER_VOLUMES='/path/to/storage:/path/to/storage /path/to/scratch:/path/to/scratch'
    bsub -Is -G compute-workshop -q workshop-interactive -a 'docker(username/dockerfile:tag)' Rscript /path/to/example.R

    export LSF_DOCKER_VOLUMES='/path/to/storage:/path/to/storage /path/to/scratch:/path/to/scratch'
    bsub -G compute-workshop -q workshop -a 'docker(username/dockerfile:tag)' Rscript /path/to/example.R

- We start with just the sample data file in our working directory.

.. image:: images/r-workshop/r.working.dir.png

- Then we run the R-script with an interactive command.

.. image:: images/r-workshop/rscript.example.gif

- And now we can see that we have graphs of the significant genes as pdf files, and the list of significant genes in the file: ``significant.genes.txt``.

.. image:: images/r-workshop/rscript.output.png

.. Note::

    - We will be covering Rstudio in another workshop.