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 You should have received an email to do this. If you didn’t, please send your email address to [email protected].

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: Indiviual Samples

In this section, you will learn how to use R to generate barplots of metagenomic data..

We have provided processed metagenomic data broken down by sample (e.g. mini_SL342389_taxa.csv and yourSample.csv), for your convenience. Each file contains taxonomically characterized and quantified metagenomics data.

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.

Next we will plot our samples as a stacked-barplot at the family level.

library(ggplot2)  

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

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

family_yourSample <- family_yourSample[family_yourSample$percent_abundance >= 0.3,] 

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

5.2.2 Run Step

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, yourSample_taxa.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 >= 0.3,] 

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

This video shows how to perform this step with a generic dataset, but the process of code execution is ths same:

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 >= 0.3,] 

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, yourSample_taxa.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 >= 0.3,] 

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

Execute with Your Sample

Create: Phylum-level Barplot

Modify and run template script

Generate a phylum-level barplot by modifying taxa_yourSample, yourSample_taxa.csv, and phyla_yourSample in the template scriptper your assigned sample.

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 >= 0.3,] 

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

Re-run: Modified Template 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)   

+

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 >= 0.3,]

+

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

Create: Family-level Barplot

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

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

5.3 Configure Linux Terminal

5.3.1 Workflow Context

Current step Linux Steps R Steps

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: Aapter Removal, Filtering

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 Metagenomics Analysis: Taxonomic Characterization & Quantification

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

5.8.2 Run Step

Launch Kraken

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:

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.

5.8.2.1 Execute Kraken

Execute with Example Sample mini_SL342389

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

Output Screen

Output Confirm the mini_SL342389_taxonomic_report.csv file is present using the ls command.

head mini_SL342389_taxonomic_report.csv

Confirm the mini_SL342389_taxonomic_report.csv file is present using the ls command.

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.

Step-by-step Guide

Execute with Your Sample

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:

head mini_my_reads_taxonomic_report.csv

5.8.2.2 Convert Kraken Output to User-friendly Format

The mini_SL342389_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

Execute with Example Sample mini_SL342389

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

Look at the first 20 lines by executing:

head -20 mini_SL342389_final_taxonomic_results.mpa

Look at the last 20 lines by executing:

tail -20 mini_SL342389_final_taxonomic_results.mpa

Execute with Your Sample

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

5.8.2.3 Extract Species Level Taxonomic Assignments

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

Execute with Example Sample mini_SL342389

grep 's__' mini_SL342389_final_taxonomic_results.mpa > mini_SL342389_species_only.tsv

Look at the first 10 lines:

head mini_SL342389_species_only.tsv

Visually confirm that all taxonmomic assignments are at the species level.

See how many species were identified:

wc -l mini_SL342389_species_only.tsv

Execute with Your Sample

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

Visually confirm that all taxonmomic assignments are at the species level.

See how many species were identified:

wc -l mini_my_reads_species_only.tsv

5.9 Identify: Abundant Species

5.9.1 Workflow Context

Current step Linux Steps R Steps
Note: because this step is completed in the R prorgamming language (not Linux Terminal), you should organize the RStudio dashboard to optimally utilize R. It should look like this:

5.9.2 Run Step

Execute with Example Sample mini_SL342389

taxa_mini_SL342389_species = read.csv('~/reads/mini_SL342389_species_only.tsv', header=FALSE, sep='\t')

taxa_mini_SL342389_species_ordered <- taxa_mini_SL342389_species[order(taxa_mini_SL342389_species$V2, decreasing = TRUE), ]  

head(taxa_mini_SL342389_species_ordered, 5)

Output

Step-by-step Guide

Insert Video

Execute with Your Sample

taxa_yourSample_species <- read.csv('~/reads/mini_my_reads_species_only.tsv', header=FALSE, sep='\t')

taxa_yourSample_species_ordered <- taxa_yourSample_species[order(taxa_yourSample_species$V2, decreasing = TRUE), ]  

head(taxa_yourSample_species_ordered, 5)

5.10 Barplot: All 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 >= 0.3,] 

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.11 Lookup: Abund. Spp. Resistance Loci

5.12 Use Species-Level Information to Determine How Similar Samples are to Each Other

5.12.2 Explanation

In this section, you will learn how to use R to generate 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. You will first see how two do this at the phylum and family levels using barplots (which you have done before) and then at the species level using a technique called Principal Component Analysis.

As humans, it’s easiest for us to see plots in two dimensions, but microbiome data is not inherently two dimensional. In a formal mathematical sense, microbiome data has a dimensionality equal to the number of microbes which can be identified – often thousands. However, most of these dimensions are highly correlated meaning, roughly, that most of the information about the abundance of some particular microbe can probably be well represented by using the abundance of a different microbe and a constant scaling factor.

In fact, it is possible to find sets of scaling factors that (likely) represent as much information as possible about a set of samples using a technique called Principal Component Analysis (or PCA for short). PCA is mathematically complex and we won’t explain it here but we will do an exercise to show how it can be used in R to plot your data in 2D.

The main idea of PCA is to “reduce the dimensionality of a data set consisting of many variables correlated with each other, either heavily or lightly, while retaining the variation present in the dataset, up to the maximum extent” 1.

What constitutes a dimension in our specific dataset? As mentioned in the instructions file: >In a formal mathematical sense, microbiome data has a dimensionality equal to the number of microbial species which can be identified – often thousands. However, most of these dimensions are highly correlated meaning, roughly, that most of the information about the abundance of some particular microbe can probably be well represented by using the abundance of a different microbe and a constant scaling factor.

As you can see, the when the PCA plot is colored by sample, you can see how well each samples type clusters and how similar they are to the other sample types.

PC = Principal Component. PC1 is the first principal component, PC 2 is the second principal component. Each principal component comprises a percentage of variance in the data. Here, we plot two principal components but many more can be computed. It’s often diminishing returns, however, and it’s not possible to (easily) visualize more than 2-3 (e.g. on X, Y, Z, axes), without doing more complicated analysis. A good rule of thumb: if a lot of variation is captured by the first two PCs, just plot those and be done:)

Sometimes you see each PCs contribution to variance in parentheses on each access; an example of this would be in the Human Microbiome Project Paper:

An interesting exercise is to figure out what other PCs are computed using the script we provided and how much variance each accounts for.

I would interpret your PCA results as follows: the mystery samples don’t cluster nicely via PC1/PC2 like the Metasub (i.e. “environmental”) and GI HMP samples do. One sample seems to be similar on PC2 but not PC1, and vice versa. A third mystery sample seems to split the difference.

library(ggplot2)   

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  

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

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

library(ggplot2)   

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')  

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

ggplot(taxa, 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)) 

library(ggplot2)
library(reshape2)

taxa = read.csv('taxa_table.csv', header=TRUE, sep=',')
taxa = taxa[taxa$rank == 'species',]
taxa = acast(taxa, sample~taxon)  # add some print statements and see if you can understand what this line does

pca = prcomp(taxa)  # run PCA on our samples. look at the pca object and get a sense of its features

principal_components = pca$x[,1:2]  # grab the first two most important components from PCA. Play with some other sets.

principal_components = data.frame(principal_components)

print(principal_components)

principal_components$sample_type = unlist(
  lapply(rownames(principal_components), FUN=function(x){strsplit(x, '\\.')[[1]][1]})
)

principal_components$sample_type = unlist(lapply(principal_components$sample_type, FUN=function(x){
  if(x == 'gastrointestinal' || x == 'mystery'){
    return(x)
  }
  return('environmental')
}))

p = ggplot(principal_components, aes(x=PC1, y=PC2, color=sample_type)) +
  geom_point(size=3) +
  labs(fill='Samples')

print(p)

6.1 Run Step

library(ggplot2)  

taxa_all_barplot <- read.csv('~/Downloads/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 >= 0.3,] 

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