Introduction

This report serves to describe all the process involved in our project workflow.

The objective is to run the various analysis scripts automatically, and, ideally, quickly.

The project starts to create a database in MongoDB with the selection of the samples to analyze from americangut database.

The samples was provided properly treated and clean, ready to be analyzed.

Download and Store requirements:

mongodb version 4.0.3 or more

python 3.5 or more

please make sure that build-essential package is installed. (sudo apt install build-essential)

Python Libraries used:

pymongo
BeautifulSoup
Pandas
os
pathlib
requests

Use the following program run sequence for sample database download and store: (please make sure to have the utilities.py with the import and download datasets.

create_project_folders.sh make it executable and run it (chmod +x)

create_mongo_collections.sh make it executable and run it (chmod +x)

Import_data_AG.py

xml_AG_tocollection.py

Import_data_T2D.py

xml_T2D_tocollection.py

Import_data_IBD.py

xml_IBD_tocollection.py

Use the following program run sequence for kraken database download and store:

create_kraken_greengenes.sh make it executable and run it (chmod +x)

Taxonomy

In this section, we will describe how we extract the taxonomy information from our samples.

In total we have 16709 samples acquired and stored in our mongoDB database.

The assign taxonomy was made with Kraken and greengenes database.

For the phylogenitics analysis with trees, we merge all the Kraken reports in one single file.

It was done with ete3tree running in python version 3.7 .

Diversity Analysis

To analyse the composition of our microbiome we proceed with the Alpha and Beta analysis.

To start we create a small dataset with 50 samples, 10 from each cavity : eye, nasopharyngeal, oral, skin and vaginal. This dataset only have 50 samples due to computacional issues.

After the selection, we transforme the kraken reports in biom tables. Then we merge all of them in one single otu table with qiime 1.9 and create a mapfile with the metadata provided.

Download for create and merge biom tables

Qiime version 1.9 (http://qiime.org/install/install.html)
Kraken-biom version 1.0.1 (https://bioconda.github.io/recipes/kraken-biom/README.html)

Script:

Convert kraken reports in biom tables:

  for file in /path_that_contais_reports/*; do kraken-biom $file --output_fp $file.biom; done

  optional: --fmt json or tsv

  Merge otus tables in Qiime 1.9:

 merge_otu_tables.py -i input_table1.biom , input_table2.biom -o merged_table.biom

Download for Statistical Analysis

 R version 3.5.2 (https://www.r-project.org/)

 R Studio (optional) https://www.rstudio.com/

Packages for R:

    Biomformat version 3.8 (https://www.bioconductor.org/packages/release/bioc/html/biomformat.html)
    Phyloseq version 3.5 (https://www.bioconductor.org/packages/release/bioc/html/phyloseq.html)

For Beta Diversity: Bray-curtis distances

Question: Are there relevant differences between microbial species between samples collected from different body sites?

Community Composition

To analyse the composition of our samples we use the microbiome 1.2.1 package for R (in this case we use Rstudio):

https://microbiome.github.io/microbiome/ and Phyloseq 3.5 .

This component will resort to bar plots and heatmaps to provide attractive results.

Analysis

We start to import the packges that we will need.

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2018 Leo Lahti et al. <microbiome.github.io>

## 
## Attaching package: 'microbiome'

## The following object is masked from 'package:base':
## 
##     transform

We will start to show the phylogenetic composition of our samples with a tree. So we import in biom format all the samples in our database but merged.

#Import Dataset | All Samples | Used for plot the general tree

tree_path <- "/home/ray2g/kraken.biom"
tree_biom <- read_biom(tree_path)

## Warning in strsplit(conditionMessage(e), "\n"): input string 1 is invalid
## in this locale

samples <- import_biom(tree_biom, parseFunction = parse_taxonomy_greengenes)

Then we do the same but only for the small dataset with the 50 samples to analyse by sample location and merged with the created mapfile.

# Import Dataset | Sample Location | 50 samples
path <- "/home/ray2g/merged_table.biom"

# Read biom table from path
table <- read_biom(path)

## Warning in strsplit(conditionMessage(e), "\n"): input string 1 is invalid
## in this locale

names(table)

##  [1] "id"                  "format"              "format_url"         
##  [4] "type"                "generated_by"        "date"               
##  [7] "matrix_type"         "matrix_element_type" "rows"               
## [10] "columns"             "shape"               "data"

# Import biom table in phyloseq
physeq <- import_biom(table, parseFunction = parse_taxonomy_greengenes)
tax_table <- tax_table(physeq)[, 1:7]
rank_names(physeq)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

sample_names(physeq)

##  [1] "SAMEA104163469" "SAMEA4495937"   "SAMEA104163465" "SAMEA104163316"
##  [5] "SAMEA104163279" "SAMEA104163534" "SAMEA104163525" "SAMEA104187573"
##  [9] "SAMEA104163396" "SAMEA104163589" "SAMEA104163327" "SAMEA96014668" 
## [13] "SAMEA94507918"  "SAMEA104206961" "SAMEA104207015" "SAMEA94527418" 
## [17] "SAMEA104207091" "SAMEA104206884" "SAMEA94417918"  "SAMEA4353811"  
## [21] "SAMEA4353790"   "SAMEA4353772"   "SAMEA3608419"   "SAMEA104206404"
## [25] "SAMEA3614155"   "SAMEA104207003" "SAMEA4353845"   "SAMEA4353787"  
## [29] "SAMEA3614284"   "SAMEA3609566"   "SAMEA3612262"   "SAMEA104206872"
## [33] "SAMEA104163596" "SAMEA3611028"   "SAMEA3614295"   "SAMEA94378168" 
## [37] "SAMEA104206940" "SAMEA3612994"   "SAMEA104163382" "SAMEA104206951"
## [41] "SAMEA104187432" "SAMEA3613042"   "SAMEA104187428" "SAMEA3613103"  
## [45] "SAMEA3628058"   "SAMEA104187430" "SAMEA94514668"  "SAMEA4353846"  
## [49] "SAMEA104187431" "SAMEA3632950"

# Import map file 
mapfile <- read.table('/home/ray2g/mapfile.txt', sep=',', comment='', head=T, row.names=1, strip.white = TRUE)
rownames(mapfile)

##  [1] "SAMEA104163469" "SAMEA4495937"   "SAMEA104163465" "SAMEA104163316"
##  [5] "SAMEA104163279" "SAMEA104163534" "SAMEA104163525" "SAMEA104187573"
##  [9] "SAMEA104163396" "SAMEA104163589" "SAMEA104163327" "SAMEA96014668" 
## [13] "SAMEA94507918"  "SAMEA104206961" "SAMEA104207015" "SAMEA94527418" 
## [17] "SAMEA104207091" "SAMEA104206884" "SAMEA94417918"  "SAMEA4353811"  
## [21] "SAMEA4353790"   "SAMEA4353772"   "SAMEA3608419"   "SAMEA104206404"
## [25] "SAMEA3614155"   "SAMEA104207003" "SAMEA4353845"   "SAMEA4353787"  
## [29] "SAMEA3614284"   "SAMEA3609566"   "SAMEA3612262"   "SAMEA104206872"
## [33] "SAMEA104163596" "SAMEA3611028"   "SAMEA3614295"   "SAMEA94378168" 
## [37] "SAMEA104206940" "SAMEA3612994"   "SAMEA104163382" "SAMEA104206951"
## [41] "SAMEA104187432" "SAMEA3613042"   "SAMEA104187428" "SAMEA3613103"  
## [45] "SAMEA3628058"   "SAMEA104187430" "SAMEA94514668"  "SAMEA4353846"  
## [49] "SAMEA104187431" "SAMEA3632950"

# Check if the sample names are the same as the file file
all(rownames(mapfile) %in% sample_names(physeq))

## [1] TRUE

# Import map file for phyloseq
sample_data <- sample_data(mapfile)

# Random Tree | Samples Location
tree_random <- rtree(ntaxa(physeq), rooted=TRUE, tip.label=taxa_names(physeq))

# Merged biom table and map file
physeq2 <- merge_phyloseq(sample_data, physeq, tree_random )
#sample_variables(physeq2)

physeq2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1735 taxa and 50 samples ]
## sample_data() Sample Data:       [ 50 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 1735 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1735 tips and 1734 internal nodes ]

# Permutacion
pseq <- physeq2

In this point we proceed with the sumarization of our samples.

#Sumarization

summarize_phyloseq(pseq)

## Compositional = NO
## 1] Min. number of reads = 62 
## 2] Max. number of reads = 261891 
## 3] Total number of reads = 1680494 
## 4] Average number of reads = 33609.88 
## 5] Median number of reads = 23319.5 
## 7] Sparsity = 0.835170028818444 
## 6] Any OTU sum to 1 or less? YES 
## 8] Number of singletons = 275 
## 9] Percent of OTUs that are singletons 15.850144092219 
## 10] Number of sample variables are: 1 
## Sample_location

nsamples(pseq)

## [1] 50

ntaxa(pseq)

## [1] 1735

#top_taxa(pseq, n = 10)
meta <- meta(pseq)
meta

The creation of the trees for all samples and the for each sample location.

# Tree code | All samples

OTU <- otu_table(samples)           
TAX <- tax_table(samples)

tree2 <- (phyloseq(OTU,TAX))

# Create random tree
random_tree = rtree(ntaxa(tree2), rooted=TRUE, tip.label=taxa_names(tree2))

# Merge phyloseq objects
tree3 <- merge_phyloseq(tree2, random_tree) # merge phyloseq objects

# Top 50
treetop50 <- names(sort(taxa_sums(tree3), decreasing = TRUE))[1:50]
tree_top50 <- transform_sample_counts(tree3, function(OTU) OTU/sum(OTU))
tree.top50 <- prune_taxa(treetop50, tree_top50)

# Plot tree | All Samples
plot_tree_all <- plot_tree(tree.top50, label.tips="taxa_names", ladderize="left", plot.margin=0.3)

# Random Tree | Samples Location
tree_random <- rtree(ntaxa(physeq), rooted=TRUE, tip.label=taxa_names(physeq))

# Top 20
top20 <- names(sort(taxa_sums(pseq), decreasing = TRUE))[1:20]
pseq_top20 <- transform_sample_counts(pseq, function(OTU) OTU/sum(OTU))
pseq.top20 <- prune_taxa(top20, pseq_top20)

# Plot tree for sample location
plot_tree_sl_genus <- plot_tree(pseq.top20, color="Genus", shape="Sample_location", label.tips="taxa_names", ladderize="right", plot.margin=0.3)

plot_tree_all

plot_tree_sl_genus

Input ete3tree plot

Alpha Diversity with the 50 samples, 10 for each cavity.

# Count
samplesum <- data.frame(sum= sample_sums(physeq2))

#Rarefation samples(optional)
ps.rarefied = rarefy_even_depth(physeq2, rngseed=1, sample.size=0.9*min(sample_sums(physeq2)), replace=F)

## `set.seed(1)` was used to initialize repeatable random subsampling.

## Please record this for your records so others can reproduce.

## Try `set.seed(1); .Random.seed` for the full vector

## ...

## 1462OTUs were removed because they are no longer 
## present in any sample after random subsampling

## ...

#Abundance
otu_tab <- t(abundances(pseq)) 

ps_abundance <- vegan::rarecurve(otu_tab, step = 50, label= FALSE, sample = min(rowSums(otu_tab), col = "blue", cex = 0.6))

# Distribution
distribution_sample <- ggplot(samplesum, aes(x=sum)) + geom_histogram(color="black", fill="indianred", binwidth = 2500) +
ggtitle("Distribution of sample depth") + xlab("Read counts") + theme(axis.title.y = element_blank())

distribution_sample

#Taxa Prevalence
taxa_prev_phylum <- plot_taxa_prevalence(ps.rarefied, "Phylum")
taxa_prev_class <- plot_taxa_prevalence(ps.rarefied, "Class")

taxa_prev_phylum

taxa_prev_class 

Beta Analysys for Samples Location.

#Set the variables of the samples to use in the colors of the plots
Sample_location <- factor(sample_data(physeq2)$Sample_location,levels = c("eye", "nasopharyngeal", "oral", "skin", "vaginal"))

# Calculate bray-curtis distance for beta diversity
distbc <- distance(physeq2, method = "bray")
ordBC <- ordinate(physeq = physeq2, method = "PCoA", distance = distbc)

# Rarefation samples(optional)
ps.rarefied = rarefy_even_depth(physeq2, rngseed=1, sample.size=0.9*min(sample_sums(physeq2)), replace=F)

## `set.seed(1)` was used to initialize repeatable random subsampling.

## Please record this for your records so others can reproduce.

## Try `set.seed(1); .Random.seed` for the full vector

## ...

## 1462OTUs were removed because they are no longer 
## present in any sample after random subsampling

## ...

#Plot for beta diversity (with Bray-Curtis distance)
plot_ordination(physeq = physeq2,ordination = ordBC,color = Sample_location) + geom_point(aes(color = Sample_location), alpha = 0.7, size = 4)

## Warning in plot_ordination(physeq = physeq2, ordination = ordBC, color =
## Sample_location): The `color` variable argument should have length equal to
## 1.Taking first value.

## Warning in plot_ordination(physeq = physeq2, ordination = ordBC, color =
## Sample_location): Color variable was not found in the available data you
## provided.No color mapped.

## No available covariate data to map on the points for this plot `type`

population_density <- plot_landscape(ps.rarefied, "NMDS", "bray", col = "Sample_location") + labs(title = paste("NMDS / Bray-Curtis"))

population_density + scale_color_brewer(palette= "Dark2") + scale_fill_gradient(low = "#e0ecf4", high = "#6e016b") + theme_classic()

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

Definitions of Tops

# Tops

top30<- names(sort(taxa_sums(pseq), decreasing = TRUE))[1:30]
pseq_top30 <- transform_sample_counts(pseq, function(OTU) OTU/sum(OTU))
pseq.top30 <- prune_taxa(top30, pseq_top30)

top50<- names(sort(taxa_sums(pseq), decreasing = TRUE))[1:50]
pseq_top50 <- transform_sample_counts(pseq, function(OTU) OTU/sum(OTU))
pseq.top50 <- prune_taxa(top50, pseq_top50)

top100<- names(sort(taxa_sums(pseq), decreasing = TRUE))[1:100]
pseq_top100 <- transform_sample_counts(pseq, function(OTU) OTU/sum(OTU))
pseq.top100 <- prune_taxa(top100, pseq_top100)

# Head Taxa

head(get_taxa_unique(pseq.top20, "Kingdom"))

## [1] "Bacteria"

head(get_taxa_unique(pseq.top20, "Phylum"))

## [1] "Firmicutes"     "Actinobacteria" "Proteobacteria" "Cyanobacteria" 
## [5] "Bacteroidetes"

head(get_taxa_unique(pseq.top20, "Class"))

## [1] "Clostridia"          "Actinobacteria"      "Gammaproteobacteria"
## [4] "Chloroplast"         "Bacteroidia"         "Bacilli"

head(get_taxa_unique(pseq.top20, "Order"))

## [1] "Clostridiales"   "Actinomycetales" "Pasteurellales"  "Streptophyta"   
## [5] "Bacteroidales"   "Bacillales"

head(get_taxa_unique(pseq.top20, "Family"))

## [1] "Veillonellaceae"    "Micrococcaceae"     "Corynebacteriaceae"
## [4] "[Tissierellaceae]"  "Pasteurellaceae"    NA

head(get_taxa_unique(pseq.top20, "Genus"))

## [1] "Veillonella"     "Rothia"          "Corynebacterium" "Finegoldia"     
## [5] NA                "Staphylococcus"

head(get_taxa_unique(pseq.top20, "Species"))

## [1] NA             "copri"        "mucilaginosa"

Creation and Vizualization of barplots for each location

# BarPlot for each location

plot_bar(pseq.top20, fill = "Kingdom") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill = "Phylum") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill = "Class") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill = "Order") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill = "Family") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill = "Genus") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

plot_bar(pseq.top20, fill= "Species") + facet_wrap(~Sample_location, scales = "free_x", nrow = 1)

Creation and Vizualization of barplots for all the samples with our dataset, the 50 samples one.

# BarPlot 

plot_bar(pseq.top20, fill="Kingdom")

plot_bar(pseq.top20, fill="Phylum")

plot_bar(pseq.top20, fill="Class")

plot_bar(pseq.top20, fill="Order")

plot_bar(pseq.top20, fill="Family")

plot_bar(pseq.top20, fill="Genus")

plot_bar(pseq.top20, fill="Species")

At the start of this analysis composition report we used the plot_composition function from microbiome package but beyond show the otu names by ids, it dont assume the taxonomic level for each domain.

So i let them in the report to analyse which may not have worked correctly.

Chatting in the forums it was said that

" The plot_composition returns a ggplot object. My suggestion is to use standard ggplot utilities for more fine-tuned figure manipulation. The plot_composition and other similar functions in the microbiome pkg are designed to facilitate fast data exploration. Fine-tuning needs are so variable and user dependent that it is not efficient to try to implement most of those solutions in the function itself since manipulation tools for ggplot objects exist already."

ggplot did not recognize the sample type provided so i dont proceed.

# Plot Composition for sample location

psl_composition_Kingdom <- plot_composition(pseq.top20, taxonomic.level = "Kingdom", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Kingdom") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Phylum <- plot_composition(pseq.top20, taxonomic.level = "Phylum", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Phylum") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Class <- plot_composition(pseq.top20, taxonomic.level = "Class", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Class") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Order <- plot_composition(pseq.top20, taxonomic.level = "Order", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Order") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Family <- plot_composition(pseq.top20, taxonomic.level = "Family", average_by = "Sample_location", transform = "compositional") +
     labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Family") +
     geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Genus <- plot_composition(pseq.top20, taxonomic.level = "genus", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Genus") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Species <- plot_composition(pseq.top20, taxonomic.level = "genus", average_by = "Sample_location", transform = "compositional") +
  labs(x = "Sample Location", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Species") +
  geom_point(color="black", position="jitter", size=1) + theme_classic()

psl_composition_Kingdom

psl_composition_Phylum

psl_composition_Class

psl_composition_Order

psl_composition_Family

psl_composition_Genus

psl_composition_Species

# BarPlot

bp_Kingdom <-plot_composition(pseq.top20,
                             taxonomic.level = "Kigdom", plot.type = "barplot") + 
  labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Kingdom") + 
  theme_ipsum(grid="Y")+ theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))


bp_Phylum <-plot_composition(pseq.top20,
                             taxonomic.level = "Phylum", plot.type = "barplot") +
  labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Phylum") + 
  theme_ipsum(grid="Y") + theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))


bp_Class <-plot_composition(pseq.top20,
                             taxonomic.level = "Class", plot.type = "barplot") +
  labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Class") + 
  theme_ipsum(grid="Y") + theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

bp_Order <-plot_composition(pseq.top20,
                             taxonomic.level = "Order", plot.type = "barplot") +
  labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Order") + 
  theme_ipsum(grid="Y")+ theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))


bp_Family <-plot_composition(pseq.top20,
                             taxonomic.level = "Family", plot.type = "barplot") + 
                             labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Family") + 
                             theme_ipsum(grid="Y")+ theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))



bp_Genus <- plot_composition(pseq.top20,
                      taxonomic.level = "Genus", plot.type = "barplot") +
                      labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Genus") + 
                      theme_ipsum(grid="Y")+ theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))


bp_Species <- plot_composition(pseq.top20,
                                taxonomic.level = "Species", plot.type = "barplot") + 
                                labs(x = "Samples", y = "Relative abundance (%)", title = "Relative abundance data", subtitle = "Species") + 
                                theme_ipsum(grid="Y")+ theme_grey() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

bp_Kingdom

bp_Phylum

bp_Class

bp_Order

bp_Family

bp_Genus

bp_Species

HeatMaps

# Heatmap

#k_bacteria <- subset_taxa(pseq, Kingdom =="Bacteria")
#heatmap(otu_table(k_bacteria), scale="column", col = terrain.colors(256), main="Heatmap | Kingdom |Bacteria")

# Agregate by Domain
pseq_kingdom <- aggregate_taxa(pseq, "Kingdom")
pseq_phylum <- aggregate_taxa(pseq, "Phylum")
pseq_class <- aggregate_taxa(pseq, "Class")
pseq_order <- aggregate_taxa(pseq, "Order")
pseq_family <- aggregate_taxa(pseq, "Family") 
pseq_genus <- aggregate_taxa(pseq, "Genus")
pseq_Species <- aggregate_taxa(pseq, "Species")

plot_heatmap(pseq_kingdom, taxa.label="Kingdom", title = " HeatMap | Kingdom")

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_phylum, taxa.label="Phylum", title = " HeatMap | Phylum")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_class, taxa.label="Class", title = " HeatMap | Class")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_order, taxa.label="Order", title = " HeatMap | Order")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_family, taxa.label="Family", title = " HeatMap | Family")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_genus, taxa.label="Genus", title = " HeatMap | Genus")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(pseq_Species, taxa.label="Species", title = " HeatMap | Species")

## Warning: Transformation introduced infinite values in discrete y-axis

Session Info

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
## 
## Matrix products: default
## BLAS: /home/ray2g/miniconda3/envs/rstudio/lib/R/lib/libRblas.so
## LAPACK: /home/ray2g/miniconda3/envs/rstudio/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] bindrcpp_0.2.2     ape_5.2            tinytex_0.9       
## [4] RColorBrewer_1.1-2 microbiome_1.2.1   hrbrthemes_0.5.0.1
## [7] ggplot2_3.1.0      phyloseq_1.24.2    biomformat_1.8.0  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.19        lattice_0.20-35     tidyr_0.8.1        
##  [4] Biostrings_2.48.0   assertthat_0.2.0    digest_0.6.18      
##  [7] foreach_1.4.4       R6_2.3.0            plyr_1.8.4         
## [10] stats4_3.5.1        evaluate_0.12       pillar_1.3.0       
## [13] zlibbioc_1.26.0     rlang_0.2.2         lazyeval_0.2.1     
## [16] data.table_1.11.8   extrafontdb_1.0     vegan_2.5-2        
## [19] S4Vectors_0.18.3    Matrix_1.2-14       rmarkdown_1.11     
## [22] labeling_0.3        splines_3.5.1       extrafont_0.17     
## [25] stringr_1.3.1       igraph_1.2.2        munsell_0.5.0      
## [28] compiler_3.5.1      xfun_0.4            pkgconfig_2.0.2    
## [31] BiocGenerics_0.26.0 multtest_2.36.0     mgcv_1.8-24        
## [34] htmltools_0.3.6     tidyselect_0.2.5    tibble_1.4.2       
## [37] IRanges_2.14.12     codetools_0.2-15    permute_0.9-4      
## [40] crayon_1.3.4        dplyr_0.7.7         withr_2.1.2        
## [43] MASS_7.3-51         grid_3.5.1          nlme_3.1-137       
## [46] jsonlite_1.5        Rttf2pt1_1.3.7      gtable_0.2.0       
## [49] magrittr_1.5        scales_1.0.0        stringi_1.2.4      
## [52] XVector_0.20.0      reshape2_1.4.3      Rhdf5lib_1.2.1     
## [55] iterators_1.0.10    tools_3.5.1         ade4_1.7-13        
## [58] Biobase_2.40.0      glue_1.3.0          purrr_0.2.5        
## [61] parallel_3.5.1      survival_2.42-6     yaml_2.2.0         
## [64] colorspace_1.3-2    rhdf5_2.24.0        cluster_2.0.7-1    
## [67] knitr_1.21          bindr_0.1.1