(1) Login to your AWS-hosted RStudio Instance (URL/username in the Google Sheet)

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).

(2) Generate barplot(s) using subset of metagenomic samples

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 how we make this plot in R:

library(ggplot2)   # These lines load additional libraries into our environment

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  # Read our taxonomic table into a computational object from a file

taxa = taxa[taxa$rank == 'phylum',]  # This filters our taxonomic table to a specific taxonomic rank. One of Kingdom, Phylum, Class, Order, Family, Genus, Species. Play around with a few other ranks.

ggplot(taxa, aes(x=sample, y=percent_abundance, fill=taxon)) + # this creates a ggplot object. Can you figure out what the aes(...) section is doing?
  geom_bar(stat="identity") +  # this tells ggplot how we want our data to be displayed
  xlab('Sample') +  # These lines tell ggplot what our axis labels should be
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) # this rotates the x-axis text 90 degrees

Here is a video of this step being performed:

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

Chunk 1:

library(ggplot2)   # These lines load additional libraries into our environment

Chunk 1+2:

library(ggplot2)   # These lines load additional libraries into our environment

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  # Read our taxonomic table into a computational object from a file

Chunk 1+2+3:

library(ggplot2)   # These lines load additional libraries into our environment

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  # Read our taxonomic table into a computational object from a file

taxa = taxa[taxa$rank == 'phylum',]  # This filters our taxonomic table to a specific taxonomic rank. One of Kingdom, Phylum, Class, Order, Family, Genus, Species. Play around with a few other ranks.

Chunk 1+2+3+4:

library(ggplot2)   # These lines load additional libraries into our environment

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  # Read our taxonomic table into a computational object from a file

taxa = taxa[taxa$rank == 'phylum',]  # This filters our taxonomic table to a specific taxonomic rank. One of Kingdom, Phylum, Class, Order, Family, Genus, Species. Play around with a few other ranks.

ggplot(taxa, aes(x=sample, y=percent_abundance, fill=taxon)) + # this creates a ggplot object. Can you figure out what the aes(...) section is doing?
  geom_bar(stat="identity") +  # this tells ggplot how we want our data to be displayed
  xlab('Sample') +  # These lines tell ggplot what our axis labels should be
  ylab('Abundance') +
  labs(fill='Phylum') +
  theme(axis.text.x = element_text(angle = 90)) # this rotates the x-axis text 90 degrees

(4) Generate barplot at the family level

2. Modify code: edit your script to generate a barplot at the “family” level, by modifying this line of code:

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

Checkpoint (CP) 1 Please upload the following to the Google Sheet:

1. Barplot at the phylum level. (To download this plot as a file, click the (+) Zoom button, right click on the plot and select “Save Image As”.)

2. Upload barplot at the family level. (To download this plot as a file, click the (+) Zoom button, right click on the plot and select “Save Image As”.)

3. Answer to these questions:

• 1. Which line of code ingests data from the taxa_table.csv file into the taxa object (copy paste the exact line of code)?
• 2. Which line of code allows you to filter down to a specific taxonomic level?
• 3. What is the $ operator doing in this line of code: taxa = taxa[taxa$rank == 'phylum',] ?
• 4. Why wouldn’t you want to plot a barplot at the species level? (If you’re not sure, try it!)
• 5. What can the barplots tell you about the similarities and dissimilarities of microbiome samples?
• 6. What questions do barlots leave unanswered about microbiome samples?

(5.) Access your Linux terminal via your AWS-hosted RStudio instance

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:

(6.) 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 behaviour of the command. There are many online resources that tell you which options a particular command can take, and how each option modifies the behaviour 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.

Checkpoint (CP) 2 Please take a screenshot of the last 15 lines of your history output and upload it to the Google Sheet.

(7) Look at your assigned dataset (time estimate: 10 min)

To make sure you undertstand the workflow, we’d like to first have you complete the analysis with a small dataset. We’ve taken 50,000 reads from the read1 and read2 files of each assigned sample. They are labeled with the prefix mini_ and located in the reads directory (files without the prefix mini are the original, full-size files, which each contain millions of reads.)

Data from each Metasub sample is contained in two files—read1 (i.e. R1) and read2 (i.e.R2). This is because the researchers paired-end sequenced each sample. When samples are sequenced on an Illumina genomic sequencer researchers get to decide: (1) the number of reads to generate for each sample library; (2) the read length (within a range of ~50 - ~300 nucleotides); and (3) whether to generate a single read or a pair of reads from each sequenced DNA library fragment.

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

A FASTQ record has the following format:

1. A line starting with @, containing the sequence ID.
2. One or more lines that contain the sequence.
3. A new line starting with the character +, and being either empty or repeating the sequence ID.
4. One or more lines that contain the quality scores.

Here’s an example of a FASTQ file with two records:

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

Navigate to the reads directory (which contains the “mini” fastq data):

$ cd ~/reads/

PLEASE NOTE: make sure to complete all steps in this section from the reads directory.

Take a look at the first 10 lines of your assigned fastq files:

$ zcat mini_my_reads.R1.fastq.gz | head

$ zcat mini_my_reads.R2.fastq.gz | head

(Please note: mini_my_reads is a placeholder for your assigned miniature reads dataset.)

Checkpoint (CP) 3

In the Google Sheet, please:

1. Upload a screenshot of your zcat results.

2. Answer these questions:

• 1. How are fastq files are structured? (i.e. What does line 1, line 2, line 3, and line 4, of each read signify?)
• 2. Why are there two reads files for each sample?

(8) Remove adapters and filter-out low quality reads

Check AdapterRemoval is working by opening your terminal and typing:

$ AdapterRemoval -h

You should see a long list of all the arguments that can be given to this tool.

Run AdapterRemoval on your mini sample (executable version):

$ 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

Please note: mini_my_reads is a placeholder for your assigned sample ID. The file extension .fq is a widely used abbreviation for .fastq.

Run AdapterRemoval on your sample (explanatory version):

$ AdapterRemoval
    --file1 mini_my_reads.R1.fastq.gz --file2 mini_my_reads.R2.fastq.gz  # this tells the program where to find our input data.
    --trimns  # this tells the program to remove reads with ambiguous bases or 'N's
    --trimqualities --minquality 40  # this tells the program to remove reads with low-quality bases
    --threads 2  # this tells the computer to multithread the program with 2 cores
    --gzip  # this tells the program that both the input and output files are compressed with a program call gzip
    --output1 clean.mini_my_reads.R1.fastq.gz  --output2 clean.mini_my_reads.R2.fastq.gz  # this tells the program where to put its output

Now we have read files that are clean. Look at one of the output files by running:

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

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

For comparison, look again at the first 10 lines of your raw reads:

$ zcat mini_my_reads.R1.fastq.gz | head

$ zcat mini_my_reads.R2.fastq.gz | head

Notice any differences between the Clean and Raw reads? If so, write them down.

Count the number of reads in each of your cleaned fastq files by running:

$ zcat clean.mini_reads.R1.fastq.gz | wc -l and manually dividing by four.

$ zcat clean.mini_reads.R2.fastq.gz | wc -l and manually dividing by four.

Count the number of reads in each of your raw (i.e. initial) fastq files by running:

$ zcat mini_reads.R1.fastq.gz | wc -l and manually dividing by four.

$ zcat mini_reads.R2.fastq.gz | wc -l and manually dividing by four.

CP 4 In the Google Sheet, please: 1. Upload a screenshot of your zcat results. 2. Answer these questions:

• 1. Do you notice any differences in the read lengths between your Raw and Cleaned reads? If so, what?
• 2. Report the numbers of your Raw and Cleaned reads.

(9) 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.

(10) Execute Metagenomic Analysis

Launch Kraken (executable version):

$ conda install --channel bioconda --yes kraken 

Launch Kraken (explanatory version):

$ conda  # the name of the package manager
install  # this is a subcommand, miniconda has several other functions
--channel bioconda
--yes  # give conda permission to install software
kraken  # the name of the program we are installing

You should now be able to test Kraken by running:

$ kraken --help

If everything is working, we can now perform taxonomic classification on our data using this command (executable version):

$ 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

BE SURE to replace mysampleID with your assigned sample ID.

If everything is working, we can now perform taxonomic classification on our data using this command (explanatory version):

$ kraken
--gzip-compressed  # tells kraken that our reads are compressed
--fastq-input  # tells kraken that our reads are formatted as fastq files
--threads 2  # uses 2 threads
--paired  # tells kraken that we have paired end data
--preload  # loads the database into RAM before running
--db ~/databases/minikraken_20171101_8GB_dustmasked  # tells kraken where to find 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

Once this command runs, take a look at the output by running:

$ head mini_my_reads_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.

The mini_my_reads_taxonomic_report.csv is formatted for storage efficiency and is not very useful for biological interpretation. Convert this file to an .mpa format so that it can be interpreted more easily:

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

Look at the first 20 lines by executing:

$ head -20 mini_my_reads_final_taxonomic_results.mpa

Look at the last 20 lines by executing:

$ tail -20 mini_my_reads_final_taxonomic_results.mpa

Pull out the species level assignments in your .mpa file and store them in a .txt file:

grep 's__' mini_my_reads_final_taxonomic_results.mpa > mini_my_reads_species_only.tsv

Look at the first 10 lines:

head mini_my_reads_species_only.txt

See how many species were identified:

wc -l mini_my_reads_species_only.txt

CP 5

In the Google Sheet, please:

1. Upload a screenshot of your kraken command output. (i.e. the one that states the % sequences classified and unclassified)

2. Upload a screenshot of your head -20 mini_my_reads_final_taxonomic_results.mpa command.

3. Upload a screenshot of your tail -20 mini_my_reads_final_taxonomic_results.mpa command.

4. Answer these questions:

• 1. What do d, p, and other abbreviations stand for?
• 2. What do you the numbers at the end of each line signify?
• 3. Why do you think Kraken can’t characterize some reads at all?
• 4. Why do you think Kraken can’t characterize some reads down to the species level?
• 5. When identifying only the species-level taxonomic assignments, what “search term” do you use in your code?

(11) Generate barplot using all kraken samples that you and your classmates analyzed

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)) 

CP 6

Please upload the following to the Google Sheet:

1. Barplot at the phylum level for all samples. (To download this plot as a file, click the (+) Zoom button, right click on the plot and select “Save Image As”.)

(12) Sort species abundance from highest to lowest

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

# Read-in all metasub data
taxa_all = read.csv('kraken_all_metasub.csv', header=TRUE, sep=',')

# Filter down to species level
taxa_all_spp = taxa_all[taxa_all$rank == 'species',] 

# Sort to your sample. Be sure to specify your SampleID (e.g. SL342389) where indicated.
taxa_all_spp_mySampleID = taxa_all_spp[taxa_all_spp$sample == 'SampleID',]  

# Order object from highest-to-lowest percent abundance 
taxa_all_spp_mySampleID_ordered <- taxa_all_spp_mySampleID[order(taxa_all_spp_mySampleID$percent_abundance, decreasing = TRUE), ]  

# Print top 5 most abundant species in your sample to the Console
head(taxa_all_spp_mySampleID_ordered, 5)

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.

CP 7

In the Google Sheet:

1. List the top 5 most abundant species in your sample.

(13) ID Potential Antibiotic Resistance Loci in the SPP in Your Assigned Sample

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):

CP 8

In the Google Sheet:

1. List three efflux pump conferring antibiotic resistance locus IDs for the species you queried using the BV-BRC.
2. Hypothesize how you would figure out if a resistant strain of one of the abundant species in your sample was present.