Graphical Analysis of Amplicon Metagenomics Data

Just some housekeeping

Load libraries for RStudio

Open the phyloseq object

ps <- readRDS("C:/Users/faysmith/Desktop/seq/ps.rds")
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 555 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 555 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 555 tips and 553 internal nodes ]

Filter out everything but the kingdom of fungi, since it was not our intent to amplify anything else

ps_sub <- subset_taxa(ps, Kingdom == "k__Fungi")
ps_sub

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 549 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 549 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 549 tips and 547 internal nodes ]

Looks like we only lost 6 taxa from the OTU table, meaning our amplification seems to have worked well to only target fungi.

Normalize the OTUs to account for sampling/sequence differences

There is some randomization at work here, so better set some seed (so the randomized numbers are reproducable) rarefy_even_depth is a function that will randomly sample OTUs from each group until it reaches the same number of OTUs that is in the smallest sample size. There are some that claim that rarefaction should not be used on this kind of data, but everyone does it anyway. Some fancier methods of normalizing the data exist, but they are more complicated and not as well published on.

set.seed(100)
#Randomly sample to the lowest denomination
ps_sub <- rarefy_even_depth(ps_sub, rngseed = TRUE)
ps_sub

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 503 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 503 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 503 tips and 501 internal nodes ]

We went from 549 OTUs in our dataset to 503 after normalizing for sample size.

#Optional: Limit all analyses to the top 10 most abundant taxa I’m not doing this one just yet, but it can make a first-look anlyses a little easier to understand.

# TopNOTUs = names(sort(taxa_sums(ps_sub), TRUE)[1:50])
# ps_10 <- prune_taxa(TopNOTUs, ps_sub)
# plot_bar(ps_10, "Depth", fill = "Genus", facet_grid = ~Area)
# 
# ps_sub <- ps_10

Stacked barplot by taxa

Plot a stacked barchart that shows the distribution of taxa (at all levels). Extremely small groups (<2%) are removed for simplification and all groups are plotted as a % of the total taxa in each sample.

#Plot the Phylum composition of Soil Fungi
ps_phylum <- ps_sub %>%
  tax_glom(taxrank = "Phylum") %>%                     # agglomerate at phylum level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > 0.02) %>%                         # Filter out low abundance taxa
  arrange(Phylum)                                      # Sort data frame alphabetically by phylum
ggplot(ps_phylum, aes(x = Depth, y = Abundance, fill = Phylum)) + 
  facet_grid(Area~Replication) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels=c("0-10 cm", "10-20 cm"))+
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance (Phyla > 2%) \n") +
  ggtitle("Phylum Composition of Soil Fungi \n Fungal Communities by Sampling Site and Depth") 

#Plot the Class composition of Soil Fungi
ps_Class <- ps_sub %>%
  tax_glom(taxrank = "Class") %>%                     # agglomerate at Class level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > 0.02) %>%                         # Filter out low abundance taxa
  arrange(Class)                                      # Sort data frame alphabetically by Class
ggplot(ps_Class, aes(x = Depth, y = Abundance, fill = Class)) + 
  facet_grid(Area~Replication) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels=c("0-10 cm", "10-20 cm"))+
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance (Class > 2%) \n") +
  ggtitle("Class Composition of Soil Fungi \n Fungal Communities by Sampling Site and Depth") 

#Plot the Order composition of Soil Fungi
ps_Order <- ps_sub %>%
  tax_glom(taxrank = "Order") %>%                     # agglomerate at Order level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > 0.02) %>%                         # Filter out low abundance taxa
  arrange(Order)                                      # Sort data frame alphabetically by Order
ggplot(ps_Order, aes(x = Depth, y = Abundance, fill = Order)) + 
  facet_grid(Area~Replication) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels=c("0-10 cm", "10-20 cm"))+
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance (Order > 2%) \n") +
  ggtitle("Order Composition of Soil Fungi \n Fungal Communities by Sampling Site and Depth")

#Plot the Family composition of Soil Fungi
ps_Family <- ps_sub %>%
  tax_glom(taxrank = "Family") %>%                     # agglomerate at Family level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > 0.02) %>%                         # Filter out low abundance taxa
  arrange(Family)                                      # Sort data frame alphabetically by Family
ggplot(ps_Family, aes(x = Depth, y = Abundance, fill = Family)) + 
  facet_grid(Area~Replication) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels=c("0-10 cm", "10-20 cm"))+
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance (Family > 2%) \n") +
  ggtitle("Family Composition of Soil Fungi \n Fungal Communities by Sampling Site and Depth")

#Plot the Genus composition of Soil Fungi
ps_Genus <- ps_sub %>%
  tax_glom(taxrank = "Genus") %>%                     # agglomerate at Genus level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > 0.02) %>%                         # Filter out low abundance taxa
  arrange(Genus)                                      # Sort data frame alphabetically by Genus
ggplot(ps_Genus, aes(x = Depth, y = Abundance, fill = Genus)) + 
  facet_grid(Area~Replication) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels=c("0-10 cm", "10-20 cm"))+
  guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance (Genus > 2%) \n") +
  ggtitle("Genus Composition of Soil Fungi \n Fungal Communities by Sampling Site and Depth")

Visually, it seems that there are differences between the samples, even the ones collected from the same area. It might be worth exploring within each class represented in this bargraph to more clearly see what differences in genus abundance there are among the samples. It would also be great to seperate these into functional groups, something that seems possible but a new program may have to be written for it to happen.

NMDS of Bray-Curtis Distances

This is a non-metric multidimensional scaling plot that uses the dissimilarities among each sample data set to map out distances. Samples are analyzed as a sum of all the sequences within and Taxa are scaled seperatly by unique sequences.

Square root transformation
Wisconsin double standardization
Run 0 stress 0.08932695 
Run 1 stress 0.08607235 
... New best solution
... Procrustes: rmse 0.1994656  max resid 0.3850024 
Run 2 stress 0.07945461 
... New best solution
... Procrustes: rmse 0.1334179  max resid 0.3312774 
Run 3 stress 0.08607371 
Run 4 stress 0.07945426 
... New best solution
... Procrustes: rmse 0.0003676254  max resid 0.0006650386 
... Similar to previous best
Run 5 stress 0.08888285 
Run 6 stress 0.08932516 
Run 7 stress 0.1318923 
Run 8 stress 0.08888249 
Run 9 stress 0.08888243 
Run 10 stress 0.08674561 
Run 11 stress 0.08888329 
Run 12 stress 0.08607378 
Run 13 stress 0.1079878 
Run 14 stress 0.1318941 
Run 15 stress 0.07945433 
... Procrustes: rmse 0.0001557373  max resid 0.000284411 
... Similar to previous best
Run 16 stress 0.07945426 
... Procrustes: rmse 6.085333e-05  max resid 0.0001132677 
... Similar to previous best
Run 17 stress 0.0893256 
Run 18 stress 0.1079886 
Run 19 stress 0.1607972 
Run 20 stress 0.07945433 
... Procrustes: rmse 0.0001964965  max resid 0.0003594765 
... Similar to previous best
*** Solution reached

This ordination plot seems to indicate clustering of samples according to the area from which they were taken. There also might be some clustering of taxa within classes around these points, indicating a possible difference in community structure between mounds, depressions, and the adjacent pasture.

Plot alpha diversity levels

I wanted to see these done for the groups of sample area and by depth (not single samples). To create this average, I merged the replications together using “merge_samples” function in phyloseq. This causes the other values to be currupted, so they are renamed in the following code as categorical variables.

[1] "SampleID"    "Depth"       "Area"        "Replication"

It is important to note that, for each of these measurements of diversity, the lower numbers indicate higher levels of diversity. For instance, the Simpson diversity index of “1” indicates no diversity and an index of “0” indicates infinite diversity. Clearly, the top 10 cm of the mound area has the highest level of diversity, as to be expected.

Plot the phylogenetic tree for all samples

N = 50
GPN = prune_taxa(names(sort(taxa_sums(ps_sub), decreasing = TRUE)[1:N]), 
    ps_sub)
p <- plot_tree(GPN, nodelabf = nodeplotblank, ladderize="left", color="Area", shape = "Depth", label.tips = "Genus", text.size = 2.5)
# p + coord_polar(theta = "y") #Opttion for doing a radial grap
p

There are clear groups that fall within the samples in the “depression” area, although they are “NA” for genus, meaning that they are identified as fungi, but do not have a clear enough match to describe them further than Phylum in some cases. Interesting that the NA seems to be more common in the “Depression”" areas.

Plot heatmap

Shows the prevalence of OTUs for each sample time

N = 200
GPN = prune_taxa(names(sort(taxa_sums(ps_sub), decreasing = TRUE)[1:N]), 
    ps_sub)
GPN1 = subset_samples(ps_sub, Depth=="1")
GPN1 = prune_taxa(names(sort(taxa_sums(GPN1), decreasing = TRUE)[1:N]), 
    GPN1)
GPN2 = subset_samples(GPN, Depth=="2")
GPN2 = prune_taxa(names(sort(taxa_sums(GPN2), decreasing = TRUE)[1:N]), 
    GPN2)
plot_heatmap(GPN, sample.label = "Area", taxa.label = "Genus")

plot_heatmap(GPN1, sample.label = "Area", taxa.label = "Genus", title="0-10 cm Depth")

plot_heatmap(GPN2, sample.label = "Area", taxa.label = "Genus", title="10-20 cm Depth")

There is clearly something going on in the 0-10 cm depth microbial community. It is odd that one of the mound replications contained an abundance in taxa that the other samples did not have.

In order to seperate the fungi by functional guilds, we will call in the whole FUNGuild database through the following code:

parse_funguild <- function(url = 'http://www.stbates.org/funguild_db.php', tax_name = TRUE){
  # require(XML)
  # require(jsonlite)
  # require(RCurl)
  ## Parse data
  tmp <- XML::htmlParse(url)
  tmp <- XML::xpathSApply(doc = tmp, path = "//body", fun = XML::xmlValue)
  ## Read url and convert to data.frame
  db <- jsonlite::fromJSON(txt=tmp)
  ## Remove IDs
  db$`_id` <- NULL
  if(tax_name == TRUE){
    ## Code legend
    ## Taxon Level: A numeral corresponding the correct taxonomic level for the taxon
    taxons <- c(
      "keyword",                                                       # 0
      "Phylum", "Subphylum", "Class", "Subclass", "Order", "Suborder", # 3:8
      "Family", "Subfamily", "Tribe", "Subtribe", "Genus",             # 9:13
      "Subgenus", "Section", "Subsection", "Series", "Subseries",      # 15:19
      "Species", "Subspecies", "Variety", "Subvariety", "Form",        # 20:24
      "Subform", "Form Species")
    ## Table with coding
    taxmatch <- data.frame(
      TaxID = c(0, 3:13, 15:26),
      Taxon = factor(taxons, levels = taxons))
    ## Match taxon codes
    db$taxonomicLevel <- taxmatch[match(x = db$taxonomicLevel, table = taxmatch$TaxID), "Taxon"]
  }
  # remove rows with missing data
  # which(
  #     with(db, trophicMode == "NULL" & guild == "NULL" & growthForm == "NULL" & trait == "NULL" & notes == "NULL")
  #     )
  ## Add database dump date as attributes to the result
  attr(db, "DownloadDate") <- date()
  return(db)
}
FUNGuild <- parse_funguild()
FUNGuild

COMING SOON: Using the FUNGuild database to apply functional guild distinctions to OTU tables in phyloseq.

Currently combining files outside of R:

write.table(FUNGuild, “C:/Users/faysmith/Desktop/R metagenomics output/fun_table.txt”, sep=“”) write.table(ps_Genus, “C:/Users/faysmith/Desktop/R metagenomics output/genus_table.txt”, sep=“”)