... | @@ -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") +
|
... | | ... | |