Update Practical 3 Processing 16S rRNA amplicon data authored by Taylor Priest's avatar Taylor Priest
...@@ -160,48 +160,48 @@ Now we are going to use the produced files in order to make some plots and visua ...@@ -160,48 +160,48 @@ Now we are going to use the produced files in order to make some plots and visua
First, let's load the necessary libraries: First, let's load the necessary libraries:
library(phyloseq) > library(phyloseq)
library(ggplot2) > library(ggplot2)
library(reshape2) > library(reshape2)
library(RColorBrewer) > library(RColorBrewer)
library(vegan) > library(vegan)
library(plyr) > library(plyr)
Before we create our phyloseq object, we can import some metadata about our samples, e.g. sample type, fraction etc. Before we create our phyloseq object, we can import some metadata about our samples, e.g. sample type, fraction etc.
sample_metadata <- read.table("sample_metadata.csv", header=TRUE, sep="\t", stringsAsFactors=FALSE) > sample_metadata <- read.table("sample_metadata.csv", header=TRUE, sep="\t", stringsAsFactors=FALSE)
rownames(sample_metadata) <- sample_metadata$Sample > rownames(sample_metadata) <- sample_metadata$Sample
View(sample_metadata) > View(sample_metadata)
Now let's create our phyloseq object Now let's create our phyloseq object
physeq_marmic <- phyloseq(otu_table(asv_count_tab, taxa_are_rows=TRUE), > physeq_marmic <- phyloseq(otu_table(asv_count_tab, taxa_are_rows=TRUE),
tax_table(taxa_with_species), sample_data(sample_metadata)) tax_table(taxa_with_species), sample_data(sample_metadata))
Once you have created your phyloseq object, the datasets can be accessed by using the following syntax Once you have created your phyloseq object, the datasets can be accessed by using the following syntax
otu_table(physeq_marmic) > otu_table(physeq_marmic)
tax_table(physeq_marmic) > tax_table(physeq_marmic)
sample_data(physeq_marmic) > sample_data(physeq_marmic)
Make some basic diversity plots and rarefaction curves Make some basic diversity plots and rarefaction curves
## Diversity plots ## Diversity plots
alpha_plot_marmic <- plot_richness(physeq_marmic, measures=c("Chao1", "Shannon", "Simpson"), x="sample", nrow=3) > alpha_plot_marmic <- plot_richness(physeq_marmic, measures=c("Chao1", "Shannon", "Simpson"), x="sample", nrow=3)
alpha_plot_marmic > alpha_plot_marmic
## Rarefaction curves to provide indication of the coverage of our community ## Rarefaction curves to provide indication of the coverage of our community
rarefaction_marmic <- rarecurve(t(otu_table(physeq_marmic)), step=50, cex=0.5) > rarefaction_marmic <- rarecurve(t(otu_table(physeq_marmic)), step=50, cex=0.5)
rarefaction_marmic > rarefaction_marmic
Now let's visualise the composition of our sample communities using basic barplots. For this, we first need to calculate relative abundance Now let's visualise the composition of our sample communities using basic barplots. For this, we first need to calculate relative abundance
physeq_marmic_rel <- transform_sample_counts(physeq_marmic, function(otu) otu/sum(otu)) > physeq_marmic_rel <- transform_sample_counts(physeq_marmic, function(otu) otu/sum(otu))
physeq_marmic_rel > physeq_marmic_rel
First let's plot phylum-level composition First let's plot phylum-level composition
marmic_phylum_barplot <- plot_bar(physeq_marmic_rel, fill="Phylum") + > marmic_phylum_barplot <- plot_bar(physeq_marmic_rel, fill="Phylum") +
facet_grid(~Sample_group, scales="free_x") + facet_grid(~Sample_group, scales="free_x") +
guides(fill=guide_legend(ncol=1)) + guides(fill=guide_legend(ncol=1)) +
labs(y="Relative abundance", x="Sample") + labs(y="Relative abundance", x="Sample") +
...@@ -212,11 +212,11 @@ First let's plot phylum-level composition ...@@ -212,11 +212,11 @@ First let's plot phylum-level composition
legend.title = element_text(size=12,face="bold",colour="black"), legend.title = element_text(size=12,face="bold",colour="black"),
strip.background.x = element_rect(fill="white"), strip.background.x = element_rect(fill="white"),
strip.text.x = element_text(size=11,face="bold",colour="black")) strip.text.x = element_text(size=11,face="bold",colour="black"))
herschel_phylum_barplot > marmic_phylum_barplot
Now at the Family level Now at the Family level
marmic_family_barplot <- plot_bar(physeq_marmic_rel, fill="Family") + > marmic_family_barplot <- plot_bar(physeq_marmic_rel, fill="Family") +
facet_grid(~Sample_group, scales="free_x") + facet_grid(~Sample_group, scales="free_x") +
guides(fill=guide_legend(ncol=1)) + guides(fill=guide_legend(ncol=1)) +
labs(y="Relative abundance", x="Sample") + labs(y="Relative abundance", x="Sample") +
... ...
......