1.1 RStudio, Assigned Dataset, Checkpoint Information

Here is a Google Sheet that contains each participant:

• AWS-hosted RStudio URL, Username, and Password (General Info tab)
• Sample Assignment Information (General Info tab)
• Worksheet where you will upload your Checkpoint (CP) Results (Checkpoint (CP) Results tab)

1.2 Slack

Please join the Slack channel for this course

1.3 Troubleshooting

1. If you’re having trouble with a step analyzing your assigned sample, compare the code to the that of the example sample code to see if you can determine where you’ve gone wrong. Each coding step requires you to first reproduce the results using the Example Sample (mini_SL342389). You should always be able to copy the code using the icon in the upper righthand corner of the screen, paste it into the R Console or Linux Terminal (depending on the step) and and run it as is. You should not ask a question about the analysis of your sample if you haven’t first reproduced the results with the Example Sample.
2. If you have a question about the Example Sample Analysis, or if you’ve successfully completed the Example Sample Analysis for a given step and have a question about analyzing your sample, go to Slack and search and see if someone has asked this question.
3. If no one has previously asked your question on Slack, please post the question with text and screenshot(s) in a single message on the main channel. A TA or MILRD Assistant Mentor will reply in a thread to to that message to assist you.
4. If a threaded conversation in Slack does not resolve your question, the person helping you will join you in a Zoom room, have you share your screen, and help you troubleshoot.

5.1 Login to your RStudio Instance

You are provided an instance of the R editor RStudio on AWS already waiting for you to use. Login to your AWS-hosted RStudio instance using the URL, username, and password assigned to you in the Google Sheet.

Please make sure the RStudio user interface dashboard looks as follows: Console Pane (Lower Left Quadrant), Scripting Pane (Upper Left Quadrant), File Manager/Plots Display/Help-Tab Pane (Lower Right Quadrant), Environment/History-Tab Pane (Upper Right Quadrant).

5.2 Barplots

In this section, you will learn how to use R to generate two plots that help visualize the similarities and differences between a subset of the Metasub samples and compare it to metagenomics samples from the Human Microbiome project.

We have provided a subset of of metagnomics data (taxa_table.csv) and one with, for your convenience. This file contains taxonomically characterized and quantified metagenomics data from nine microbiome samples: three from the human microbiome project, three from the Metasub project, and three from a mystery source.

First we will plot our samples as a stacked-barplot at the phylum level. A stacked-barplot shows a set of numbers as a series of columns one on top of each other colored by a label. In our case, a taxonomic profile is a set of numerical abundances labeled by the microbial species it belongs to.

Here’s the template code to make this plot in R:

library(ggplot2)  

taxa_yourSample = read.csv('yourSample.csv', header=TRUE, sep=',') 

family_yourSample = taxa_yourSample[taxa_yourSample$rank == 'phylum',] 

family_yourSample = family_yourSample[family_yourSample$percent_abundance >= 2,] 

ggplot(family_yourSample, aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

1. These line loads the ggplot2 library into our environment. The library() function is used to load additional libraries that contain functions to be utilized by the script that follows.
3. Reads tabular taxonomic data into a computational object from a .csv file using the read.csv() function.
5. This filters our taxonomic table to a specific taxonomic rank (i.e. Kingdom, Phylum, Class, Order, Family, Genus, Species).
7. This filters out and discards data for a specified variable below the specified threshold.
   
9. this creates a ggplot() object.
10. this tells ggplot() how we want our data to be displayed
11. These lines tell ggplot() what our axis labels should be
14. this rotates the x-axis text 90 degrees

Execute with Example Sample mini_SL342389

Create: Phylum-level Barplot

Run Script: all at once

This code generates a phylum-level barplot by modifying taxa_yourSample, mini_yourSample.csv, and phyla_yourSample per Example sample.

library(ggplot2)  

taxa_mini_SL342389 = read.csv('mini_SL342389_taxa.csv', header=TRUE, sep=',') 

phyla_mini_SL342389 = taxa_mini_SL342389[taxa_mini_SL342389$rank == 'phylum',] 

phyla_mini_SL342389 = phyla_mini_SL342389[phyla_mini_SL342389$percent_abundance >= 2,] 

ggplot(phyla_mini_SL342389 , aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

Output

Step-by-Step Video

Run Script: Chunk-by-chunk

Clear your Console, Environment Tab, and Plots Tab using the 🧹 button and re-run the script chunk-by-chunk (noting what each chunk appears to be doing):

library(ggplot2)   

Note 1

Code should be sent to the Console. Library will be loaded behing the scenes. This command worked if no error messages were printed in the console)

+

taxa_mini_SL342389 = read.csv('mini_SL342389_taxa.csv', header=TRUE, sep=',')

Note 1

After execution, the taxa_mini_SL342389 object should be listed in the Environment Tab. Click on the name of this object in the environment tab and it should load as a table in a new tab in the scripting pane (as long as it’s not too big). Alternatively, you can use the head() function to return the first “parts of a vector, matrix, table, data frame or function”; to learn more about this function, execute ?head() in the Console. Go ahead and execute head(taxa_mini_SL342389) in the Console now. It will show the first 6 lines of the taxa_mini_SL342389 by default.)

Note 2

Take a look at the rank column. Even from the first few lines of the taxa_mini_SL34238 object, it’s evident there are taxonomic classifications made a several levels (e.g. phylum, class order, genus, etc). Keep this in mind as you execute the next step.

+

phyla_mini_SL342389 = taxa_mini_SL342389[taxa_mini_SL342389$rank == 'phylum',] 

Note 1

After Execution, run head(taxa_mini_SL342389) in the Console. How did this line of code transform the data?

+

phyla_mini_SL342389 = phyla_mini_SL342389[phyla_mini_SL342389$percent_abundance >= 2,] 

Note 1

In the R parlance, an object’s rows are referred to as oberservations (obs. for short); an object’s column headers are referred to as variables. How did this line of code transform the data? (Hint: look at the change in observations.)

+

ggplot(phyla_mini_SL342389 , aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

Create: Family-level Barplot

This code generates family-level barplot by modifying taxa_yourSample, mini_yourSample.csv, phyla_yourSample, and labs() per Example sample.

library(ggplot2)  

taxa_mini_SL342389 = read.csv('mini_SL342389_taxa.csv', header=TRUE, sep=',') 

family_mini_SL342389 = taxa_mini_SL342389[taxa_mini_SL342389$rank == 'family',] 

family_mini_SL342389 = family_mini_SL342389[family_mini_SL342389$percent_abundance >= 2,] 

ggplot(family_mini_SL342389 , aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Family') +
  theme(axis.text.x = element_text(angle = 90)) 

Output

library(ggplot2)  

taxa_yourSample = read.csv('yourSample_taxa.csv', header=TRUE, sep=',') 

phyla_yourSample = taxa_yourSample[taxa_yourSample$rank == 'phylum',] 

phyla_yourSample = phyla_yourSample[phyla_yourSample$percent_abundance >= 2,] 

ggplot(phyla_yourSample, aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

library(ggplot2)   

taxa_yourSample = read.csv('yourSample_taxa.csv', header=TRUE, sep=',') 

phyla_yourSample = taxa_yourSample[taxa_yourSample$rank == 'phylum',] 

phyla_yourSample = phyla_yourSample[phyla_yourSample$percent_abundance >= 2,]

ggplot(taxa, aes(x=sample, y=percent_abundance, fill=taxon)) + 
  geom_bar(stat="identity") + 
  xlab('Sample') +  
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

phyla_yourSample = taxa_yourSample[taxa_yourSample$rank == 'phylum',] 

5.3 Access your Linux terminal via RStudio

Access your Linux terminal by logging into RStudio via the web browser (can be found in the the Google Sheet.

Setup to your Linux environment:

It should look like this once you’re finished setting up your Linux environment:

5.4 Navigating your Linux terminal.

Linux is a command-line based operating system. For our purposes Linux = Unix.

Let’s do a brief review of Linux commands on our AWS-hosted Linux instance in RStudio. This the language and structure of this tutorial is heavily borrowed from http://www.ee.surrey.ac.uk/Teaching/Unix/. Go to this website and read the “Introduction to the UNIX Operating System”. For our purposes, you can treat Unix and Linux as interchangeable.

ls (list) When you first login, your current working directory is your home directory by default. Your home directory in every Linux system is denoted by the shorthand ~.

To find out what is in your home directory, type:

$ ls

The ls command ( lowercase L and lowercase S ) lists the contents of your current working directory.

You should see these Directories and files: R databases miniconda3 pass.txt reads rstudio taxa_table.csv

Note that the directories and files are colored differently.

ls does not cause all the files in your home directory to be listed—only those ones whose name does not begin with a dot (.). Files beginning with a dot (.) are known as hidden files and usually contain important program configuration information. They are hidden because you should not change them unless you are very familiar with Linux!!!

To list all files in your home directory including those whose names begin with a dot, type:

$ ls -a

As you likely guessed, the a in -a stands for all.

As you should now see, ls -a lists files that are normally hidden.

ls is an example of a command which can take options: -a is an example of an option. The options change the behavour of the command. There are many online resources that tell you which options a particular command can take, and how each option modifies the behavour of the command. We’ve posted a Linux command “Cheat Sheet” in the channel with Cody.

Making Directories

mkdir (make directory)

We will now make a subdirectory in your home directory to hold the files you will be creating and using in the course of this tutorial. To make a subdirectory calledunixstuff in your current working directory type:

$ mkdir unixstuff

To see the directory you have just created, type:

$ ls

Changing to a different directory

cd (change directory)

The command cd directory means change the current working directory to ‘directory’. The current working directory may be thought of as the directory you are in, i.e. your current position in the file-system tree.

To change to the directory you have just made, type:

$ cd unixstuff

Type ls and press Enter to see the contents (which should be empty)

Make another directory inside the unixstuff directory called backups

Still in the unixstuff directory, type:

$ ls -a

As you can see, in the unixstuff directory (and in all other directories), there are two special directories called (.) and (..)

The current directory (.) In UNIX, (.) means the current directory, so typing

$ cd .

NOTE: there is a space between cd and the dot

means stay where you are (the unixstuff directory).

This may not seem very useful at first, but using (.) as the name of the current directory will save a lot of typing, as we shall see later in the tutorial.

The parent directory (..) (..) means the parent of the current directory, so typing

$ cd ..

will take you one directory up the hierarchy (back to your home directory). Try it now.

Note: typing cd with no argument always returns you to your home directory. This is very useful if you are lost in the file system.

Pathnames pwd (print working directory) Pathnames (also called filepaths) enable you to work out where you are in relation to the whole file-system. For example, to find out the absolute pathname of your home-directory, type cd to get back to your home-directory and then type:

$ pwd

The full pathname will look something like this:

/home/user#

which means that user# (your home directory) is in the sub-directory home (the group directory), which is in the top-level root directory called ” / ” .

Unix File Structure

Use the commands cd, ls and pwd to explore the file system.

More about home directories and pathnames

Understanding pathnames First type cd to get back to your home-directory, then type:

$ ls unixstuff

to list the conents of your unixstuff directory.

Now type:

$ ls backups

You will get a message like this: backups: No such file or directory

This is because backups is not in your current working directory. To use a command on a file (or directory) not in the current working directory (the directory you are currently in), you must either cd to the correct directory, or specify its full pathname. To list the contents of your backups directory, you must type:

$ ls unixstuff/backups

~ (your home directory)

Home directories can also be referred to by the tilde ~ character. It can be used to specify paths starting at your home directory. So typing

$ ls ~/unixstuff

will list the contents of your unixstuff directory, no matter where you currently are in the file system.

What do you think ls ~ lists?

At this point, you have completed through Introduction and Tutorial One on http://www.ee.surrey.ac.uk/Teaching/Unix/. If you’re interested, you can continue with Tutorial Two (Copying Files, Moving Files, Removing Files and directories, Displaying the contents of a file on the screen, Searching the contents of a file) and subsequent Tutorials 3-8, but they are outside of the scope of the Linux skills you’ll need to complete the VTP.

Once you have finished all of the above exercises, please execute:

$ history 

The history command shows the recently executed commands in your linux terminal.

Side Note: you can save your Linux command history to a file using this command:

$ history > yourName_linux_work.txt

yourName_linux_work.txt is a file that should contain all of the commands you’ve executed in your Linux terminal.

5.5 Look: assigned dataset

5.5.4 Paired-end vs. Single-end

Here is an diagram that summarizes paired-end versus single-end sequencing:

@071112_SLXA-EAS1_s_7:5:1:817:345
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC
@071112_SLXA-EAS1_s_7:5:1:801:338
GTTCAGGGATACGACGTTTGTATTTTAAGAATCTGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII6IBI

cd ~/reads/

  zcat mini_SL342389.R1.fastq.gz | head

  zcat mini_SL342389.R2.fastq.gz | head

zcat mini_my_reads.R1.fastq.gz | head

zcat mini_my_reads.R2.fastq.gz | head

5.6 Quality Control

AdapterRemoval --file1 mini_my_reads.R1.fastq.gz --file2 mini_my_reads.R2.fastq.gz --trimns --trimqualities --minquality 40 --threads 2 --gzip --output1 clean.mini_my_reads.R1.fastq.gz  --output2 clean.mini_my_reads.R2.fastq.gz

AdapterRemoval -h

AdapterRemoval --file1 mini_SL342389.R1.fastq.gz --file2 mini_SL342389.R2.fastq.gz --trimns --trimqualities --minquality 40 --threads 2 --gzip --output1 clean.mini_SL342389.R1.fastq.gz  --output2 clean.mini_SL342389.R2.fastq.gz

AdapterRemoval --file1 mini_my_reads.R1.fastq.gz --file2 mini_my_reads.R2.fastq.gz --trimns --trimqualities --minquality 40 --threads 2 --gzip --output1 clean.mini_my_reads.R1.fastq.gz  --output2 clean.mini_my_reads.R2.fastq.gz

  zcat clean.mini_SL342389.R1.fastq.gz | head

  zcat clean.mini_SL342389.R2.fastq.gz | head

  zcat mini_SL342389.R1.fastq.gz | head

  zcat mini_SL342389.R2.fastq.gz | head

  zcat clean.mini_SL342389.R1.fastq.gz | wc -l

  zcat clean.mini_SL342389.R2.fastq.gz | wc -l

  zcat mini_SL342389.R1.fastq.gz | wc -l

  zcat mini_SL342389.R2.fastq.gz | wc -l

  zcat clean.mini_my_reads.R1.fastq.gz | head

  zcat clean.mini_my_reads.R2.fastq.gz | head

  zcat mini_my_reads.R1.fastq.gz | head

  zcat mini_my_reads.R2.fastq.gz | head

  zcat clean.mini_reads.R1.fastq.gz | wc -l

  zcat clean.mini_reads.R2.fastq.gz | wc -l

  zcat mini_reads.R1.fastq.gz | wc -l

  zcat mini_reads.R2.fastq.gz | wc -l

5.7 Not shown: Remove contaminating human DNA from your sample

Researchers often decide at this point to remove DNA from unwanted species, such as human. DNA that we don’t want is called “Contaminating DNA”. We are not going to perform this step due to time constraints.

5.8 Taxonomic Characterization & Quantification

Launch Kraken (executable version):

(base) user#@ip-#:~/reads $

conda install --channel bioconda --yes kraken 

conda is a package manager
install is a subcommand ofconda
--yes gives conda permission to install software
kraken is the program we are installing

You should now be able to test Kraken by running:

(base) user#@ip-#:~/reads $

kraken --help

--help is an option that displays the options for the tool

If a list of kraken options and explanations is displayed, it’s confirmation the tool is running. Provided this happnes, we can now perform taxonomic classification on our data using this command:

(base) user#@ip-#:~/reads $

kraken --gzip-compressed --fastq-input --threads 2 --paired --preload --db ~/databases/minikraken_20171101_8GB_dustmasked clean.mini_my_reads.R1.fastq.gz clean.mini_my_reads.R2.fastq.gz > mini_my_reads_taxonomic_report.csv

kraken command invocation
--gzip-compressed indicates our reads are gzip compressed
--fastq-input indicates our reads are formatted as fastq files
--threads tells the compute instance how many threads to use
--paired indicates we’re using paired end data
--preload specifies to load the database into RAM before running
--db specifies the path to the database
clean.mini_my_reads.R1.fastq.gz filepath to Read1 file
clean.mini_my_reads.R2.fastq.gz filepath to Read2 file
mini_my_reads_taxonomic_report.csv filepath to file where kraken will deposit results

5.8.3 Execute metagonomics analysis: See one, do one

text

See One (Example)
Sample: mini_SL342389

(base) user#@ip-#:~/reads $

kraken --gzip-compressed --fastq-input --threads 2 --paired --preload --db ~/databases/minikraken_20171101_8GB_dustmasked clean.mini_SL342389.R1.fastq.gz clean.mini_SL342389.R2.fastq.gz > mini_SL342389_taxonomic_report.csv

Confirm the file is present using the ls command.

Look at Kraken output:

(base) user#@ip-#:~/reads $

head mini_SL342389_taxonomic_report.csv

(the head command prints the first 10 lines of any file to the terminal)

It should look similar to this:

Note Column 1 indicates whether it was classified by Kraken or not. Column 2 denotes the Read ID (line 1 of the fastq read record). Subsequent columns indicated taxonomic mappings, if any.

Execute with your sample

(base) user#@ip-#:~/reads $

kraken --gzip-compressed --fastq-input --threads 2 --paired --preload --db ~/databases/minikraken_20171101_8GB_dustmasked clean.mini_my_reads.R1.fastq.gz clean.mini_my_reads.R2.fastq.gz > mini_my_reads_taxonomic_report.csv

Confirm the file is present using the ls command.

Once this command runs, take a look at the output by running: (base) user#@ip-#:~/reads $

head mini_my_reads_taxonomic_report.csv

kraken-mpa-report --db ~/databases/minikraken_20171101_8GB_dustmasked mini_SL342389_taxonomic_report.csv > mini_SL342389_final_taxonomic_results.mpa

head -20 mini_SL342389_final_taxonomic_results.mpa

tail -20 mini_SL342389_final_taxonomic_results.mpa

kraken-mpa-report --db ~/databases/minikraken_20171101_8GB_dustmasked mini_my_reads_taxonomic_report.csv > mini_my_reads_final_taxonomic_results.mpa

head -20 mini_my_reads_final_taxonomic_results.mpa

tail -20 mini_my_reads_final_taxonomic_results.mpa

grep 's__' mini_SL342389_final_taxonomic_results.mpa > mini_SL342389_species_only.tsv

head mini_SL342389_species_only.tsv

wc -l mini_SL342389_species_only.tsv

grep 's__' mini_my_reads_final_taxonomic_results.mpa > mini_my_reads_species_only.tsv

head mini_my_reads_species_only.tsv

wc -l mini_my_reads_species_only.tsv

5.9 Barplots

Now create a barplot with all of the metasub samples:

library(ggplot2)  

taxa_all_barplot = read.csv('kraken_all_metasub.csv', header=TRUE, sep=',') 

taxa_all_barplot = taxa_all_barplot[taxa_all_barplot$rank == 'phylum',] 

taxa_all_barplot = taxa_all_barplot[taxa_all_barplot$percent_abundance >= 2,] 

ggplot(taxa_all_barplot, aes(x=sample, y=percent_abundance, fill=taxon)) +
  geom_bar(stat="identity") +  
  xlab('Sample') + 
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) 

5.10 ID Abundant species

Read-in all of the Kraken data from the samples you and the other participants analyzed:

taxa_all = read.csv('kraken_all_metasub.csv', header=TRUE, sep=',')

taxa_all_spp = taxa_all[taxa_all$rank == 'species',] 

taxa_all_spp_mySampleID = taxa_all_spp[taxa_all_spp$sample == 'SampleID',]  

taxa_all_spp_mySampleID_ordered <- taxa_all_spp_mySampleID[order(taxa_all_spp_mySampleID$percent_abundance, decreasing = TRUE), ]  

head(taxa_all_spp_mySampleID_ordered, 5)

1. Read-in all metasub data
3. Filter down to species level
5. Sort to your sample. Be sure to specify your SampleID (e.g. SL342389) where indicated.
7. Order object from highest-to-lowest percent abundance
9. Print top 5 most abundant species in your sample to the Console

For your information, you can open taxa_all_spp_mySampleID as a table in the Scripting Pane by clicking on the name of the object in the Environment Pane. You can order rows based on the incresing and decreasing values of a specific column by clicking the ▲ and ▼ arrows of that column.

5.11 ID AMR Loci

Navigate to the Bacterial and Viral Bioinformatics Resource Center (BV-BRC): https://www.bv-brc.org/

Use the tools in BV-BRC to ID at least one AMR locus/gene for the top 5 most abundant species in your sample.

1. Select Taxa option from the “All Data Types” dropdown menu
2. Search for a given species via the search bar
3. Select the “species” result (usually the first one) by checking the box next to that row
4. Click TAXON OVERVIEW from the green column
5. Click the Specialty Genes tab
6. Click the blue Filter button
7. For Property, select Antibiotic Resitance
8. For Source, select PATRIC
9. For Classification, search for efflux pump conferring antibiotic resistance using the magnifying glass feature
10. Select a few of the results rows by checking their boxes
11. Click the FASTA option in the green column and select View FASTA DNA This will show you the nucleotide sequences of the selected resistance loci.

Here’s a screencast using the species Staphylococcus aureus (Methicillin-resistant Staphylococcus aureus (MRSA) is a widely known antibiotic resistant hospital acquired infection):