Changes
Page history
Update Practical 3 Processing 16S rRNA amplicon data
authored
Jan 14, 2022
by
Taylor Priest
Show whitespace changes
Inline
Side-by-side
Practical-3-Processing-16S-rRNA-amplicon-data.md
View page @
a56d69e6
...
@@ -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") +
...
...
...
...