Data

Import raw data and assign sample key:

map_corrected <- read.csv("../map_corrected.txt", sep = "\t", row.names = 1)
knitr::kable(map_corrected) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
BarcodeSequence LinkerPrimerSequence FileInput SampleType RA BVAS REMISSION_BVAS ENTAS Description
kg001 NA NA kg001-n-fischer-16s_S31_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg002 NA NA kg002-n-fischer-16s_S32_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg003 NA NA kg003-n-fischer-16s_S33_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg004 NA NA kg004-n-fischer-16s_S34_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg005 NA NA kg005-n-fischer-16s_S35_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg007 NA NA kg007-n-fischer-16s_S36_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg008 NA NA kg008-n-fischer-16s_S37_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg009 NA NA kg009-n-fischer-16s_S38_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg010 NA NA kg010-n-fischer-16s_S39_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg011 NA NA kg011-n-fischer-16s_S40_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg014 NA NA kg014-n-fischer-16s_S41_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg015 NA NA kg015-n-fischer-16s_S42_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg016 NA NA kg016-n-fischer-16s_S43_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg019 NA NA kg019-n-fischer-16s_S44_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg020 NA NA kg020-n-fischer-16s_S45_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg021 NA NA kg021-n-fischer-16s_S46_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg022 NA NA kg022-n-fischer-16s_S47_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg023 NA NA kg023-n-fischer-16s_S48_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg025 NA NA kg025-n-fischer-16s_S49_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg026 NA NA kg026-n-fischer-16s_S50_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg027 NA NA kg027-n-fischer-16s_S51_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg028 NA NA kg028-n-fischer-16s_S52_L001.assembled_filtered.nonchimera.fasta control control control control control kg
kg029 NA NA kg029-n-fischer-16s_S53_L001.assembled_filtered.nonchimera.fasta control control control control control kg
mw001 NA NA mw001-n-fischer-16s_S54_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw003 NA NA mw003-n-fischer-16s_S55_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw004 NA NA mw004-n-fischer-16s_S56_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw005 NA NA mw005-n-fischer-16s_S57_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw006 NA NA mw006-n-fischer-16s_S58_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2+ mw
mw007 NA NA mw007-n-fischer-16s_S59_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw008 NA NA mw008-n-fischer-16s_S60_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2+ mw
mw009 NA NA mw009-n-fischer-16s_S61_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS+ unknown unknown mw
mw010 NA NA mw010-n-fischer-16s_S62_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw011 NA NA mw011-n-fischer-16s_S63_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2+ 2+ mw
mw013 NA NA mw014-n-fischer-16s_S64_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- unknown unknown mv
mw014 NA NA mw014-n-fischer-16s_S65_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2+ 2+ mw
mw015 NA NA mw015-n-fischer-16s_S66_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS+ unknown unknown mw
mw016 NA NA mw016-n-fischer-16s_S67_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2- mw
mw017 NA NA mw017-n-fischer-16s_S68_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2+ 2+ mw
mw018 NA NA mw018-n-fischer-16s_S69_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS+ unknown 2+ mw
mw019 NA NA mw019-n-fischer-16s_S70_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown mw
mw021 NA NA mw021-n-fischer-16s_S71_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2- mw
mw022 NA NA mw022-n-fischer-16s_S72_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS+ unknown 2- mw
micro1 NA NA micro1-lamprecht_S1_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- unknown micro
micro3 NA NA micro3-lamprecht_S3_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2- micro
micro4 NA NA micro4-lamprecht_S4_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2- micro
micro6 NA NA micro6-lamprecht_S6_L001.assembled_filtered.nonchimera.fasta GPA GPA unknown unknown unknown micro
micro7 NA NA micro7-lamprecht_S7_L001.assembled_filtered.nonchimera.fasta GPA GPA unknown unknown unknown micro
micro12 NA NA micro12-lamprecht_S12_L001.assembled_filtered.nonchimera.fasta GPA GPA unknown unknown unknown micro
micro13 NA NA micro13-lamprecht_S13_L001.assembled_filtered.nonchimera.fasta GPA GPA unknown unknown unknown micro
micro16 NA NA micro16-lamprecht_S16_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2+ unknown micro
micro17 NA NA micro17-lamprecht_S17_L001.assembled_filtered.nonchimera.fasta GPA GPA BVAS- 2- 2- micro
ra002 NA NA gpa-ra002-n-fischer-16s_S10_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra003 NA NA ra003-n-fischer-16s_S12_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra004 NA NA ra004-n-fischer-16s_S11_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra005 NA NA ra005-n-fischer-16s_S13_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra006 NA NA ra006-n-fischer-16s_S14_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra007 NA NA ra007-n-fischer-16s_S15_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra008 NA NA ra008-n-fischer-16s_S16_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra009 NA NA ra009-n-fischer-16s_S17_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra010 NA NA ra010-n-fischer-16s_S18_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra013 NA NA ra013-n-fischer-16s_S19_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra014 NA NA ra014-n-fischer-16s_S20_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra015 NA NA ra015-n-fischer-16s_S21_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra017 NA NA ra017-n-fischer-16s_S22_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra018 NA NA ra018-n-fischer-16s_S23_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra019 NA NA ra019-n-fischer-16s_S24_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra020 NA NA ra020-n-fischer-16s_S25_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra021 NA NA ra021-n-fischer-16s_S26_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra022 NA NA ra022-n-fischer-16s_S27_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra023 NA NA ra023-n-fischer-16s_S28_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra
ra024 NA NA ra024-n-fischer-16s_S29_L001.assembled_filtered.nonchimera.fasta RA RA- RA RA RA ra
ra025 NA NA ra025-n-fischer-16s_S30_L001.assembled_filtered.nonchimera.fasta RA RA+ RA RA RA ra

Prerequisites to be installed

install.packages("dplyr")     # To manipulate dataframes
install.packages("readxl")    # To read Excel files into R
install.packages("ggplot2")   # for high quality graphics
install.packages("heatmaply")
source("https://bioconductor.org/biocLite.R")
biocLite("phyloseq")
library("readxl")  # necessary to import the data from Excel file
library("ggplot2")  # graphics
library("picante")
library("microbiome")  # data analysis and visualisation
library("phyloseq")  # also the basis of data object. Data analysis and visualisation
library("ggpubr")  # publication quality figures, based on ggplot2
library("dplyr")  # data handling, filter and reformat data frames
library("RColorBrewer")  # nice color options
library("heatmaply")
library(vegan)
library(gplots)

Read the data and create phyloseq objects

Three tables are needed * OTU * Taxonomy * Samples

# Change your working directory to where the files are located
ps.ng.tax <- import_biom("./table_mc1300.biom", "../clustering/rep_set.tre")
sample <- read.csv("../map_corrected.txt", sep = "\t", row.names = 1)
SAM = sample_data(sample, errorIfNULL = T)
# rownames(SAM) <-
# c('1','2','3','5','6','7','8','9','10','12','13','14','15','16','17','18','19','20','21','22','23','24','25','26','27','28','29','30','31','32','33','34','35','36','37','38','39','40','41','42','43','44','46','47','48','49','50','51','52','53','55')
ps.ng.tax <- merge_phyloseq(ps.ng.tax, SAM)
print(ps.ng.tax)
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 631 taxa and 60 samples ]
sample_data() Sample Data:       [ 60 samples by 9 sample variables ]
tax_table()   Taxonomy Table:    [ 631 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 631 tips and 630 internal nodes ]
colnames(tax_table(ps.ng.tax)) <- c("Domain", "Phylum", "Class", "Order", "Family",
    "Genus", "Species")
saveRDS(ps.ng.tax, "./ps.ng.tax.rds")

Visualize data

sample_names(ps.ng.tax)
 [1] "kg001"   "kg002"   "kg003"   "kg004"   "kg005"   "kg007"   "kg009"  
 [8] "kg015"   "kg016"   "kg019"   "kg020"   "kg021"   "kg022"   "kg023"  
[15] "kg025"   "kg026"   "kg027"   "kg028"   "kg029"   "micro1"  "micro3" 
[22] "micro4"  "micro6"  "micro7"  "micro12" "micro13" "micro16" "micro17"
[29] "mw001"   "mw004"   "mw005"   "mw006"   "mw007"   "mw009"   "mw010"  
[36] "mw013"   "mw014"   "mw015"   "mw017"   "mw021"   "ra002"   "ra003"  
[43] "ra004"   "ra005"   "ra006"   "ra007"   "ra008"   "ra009"   "ra010"  
[50] "ra013"   "ra014"   "ra015"   "ra017"   "ra018"   "ra019"   "ra020"  
[57] "ra022"   "ra023"   "ra024"   "ra025"  
rank_names(ps.ng.tax)
[1] "Domain"  "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(ps.ng.tax)
[1] "BarcodeSequence"      "LinkerPrimerSequence" "FileInput"           
[4] "SampleType"           "RA"                   "BVAS"                
[7] "REMISSION_BVAS"       "ENTAS"                "Description"

Normalize number of reads in each sample using median sequencing depth.

# RAREFACTION set.seed(9242) # This will help in reproducing the filtering and
# nomalisation. ps.ng.tax <- rarefy_even_depth(ps.ng.tax, sample.size = 42369)
# total <- 42369 NORMALIZE number of reads in each sample using median
# sequencing depth.
total = median(sample_sums(ps.ng.tax))
# > total [1] 42369
standf = function(x, t = total) round(t * (x/sum(x)))
ps.ng.tax = transform_sample_counts(ps.ng.tax, standf)
ps.ng.tax_rel <- microbiome::transform(ps.ng.tax, "compositional")
saveRDS(ps.ng.tax, "./ps.ng.tax.rds")
hmp.meta <- meta(ps.ng.tax)
hmp.meta$sam_name <- rownames(hmp.meta)

Heatmaps

We consider the most abundant OTUs for heatmaps. For example one can only take OTUs that represent at least 1% of reads in at least one sample. Remember we normalized all the sampples to median number of reads (total). We are left with only 168 OTUS which makes the reading much more easy.

# Custom function to plot a heatmap with the specified sample order
ps.ng.tax_abund <- phyloseq::filter_taxa(ps.ng.tax, function(x) sum(x > total * 0.01) >
    0, TRUE)
kable(otu_table(ps.ng.tax_abund)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
kg001 kg002 kg003 kg004 kg005 kg007 kg009 kg015 kg016 kg019 kg020 kg021 kg022 kg023 kg025 kg026 kg027 kg028 kg029 micro1 micro3 micro4 micro6 micro7 micro12 micro13 micro16 micro17 mw001 mw004 mw005 mw006 mw007 mw009 mw010 mw013 mw014 mw015 mw017 mw021 ra002 ra003 ra004 ra005 ra006 ra007 ra008 ra009 ra010 ra013 ra014 ra015 ra017 ra018 ra019 ra020 ra022 ra023 ra024 ra025
80113 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 8 0 0 0 0 0 0 0 0 0 49 13 5 80 59 8 4 18 112 17 0 5 15 11 22 21 9 6 111 14 0 0 0 0 0 0 0 0 0 0 0 0
521073 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 81 29 48 0 18 0 21 18 1 0 4 68 81 2 122 49 11 20 22 296 13 17 14 26 31 38 31 35 35 336 87 3 0 0 2 0 0 1 0 0 0 1 0
242293 0 0 3 0 0 0 0 0 12 0 0 0 0 0 0 0 0 12 163 0 0 1 16 0 1 1 6 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 2 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4
812921 0 0 0 0 0 0 0 0 0 0 6 2 0 0 0 0 0 3 7 0 15 0 0 1 0 0 0 12 0 0 0 0 0 0 2 1 6 0 0 0 0 0 0 150 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0
4405869 0 2 0 0 1 0 0 0 0 0 0 0 4 0 1 0 0 3 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 2 0 3 2 0 0 0 87 0 5 0 1 0 0 39 0 0 0
341460 0 2 228 0 37 0 0 7 5 0 1 1 9 0 2 10 80 0 0 2 0 1 0 1 6 0 16 0 0 0 0 0 3 3 4 18 0 0 0 4 0 0 2 0 9 0 0 0 117 35 0 18 0 102 19 0 50 0 8 0
240755 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 3498 0 2 1 2 1 0 0 0 0 0 0 6 0 1 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0
3125352 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 147 0 0 0 0 0 0 0 0 0 0 286 2 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
237745 0 0 0 0 28 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 115 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
637049 418 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
656889 332 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
538000 0 0 0 27 0 0 0 384 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
345362 63 0 0 0 4 0 14 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 444 0 0 0 9 0 16 0 2 0 230 9 0 10 0 0 0 2 18 175 0 1 101 3 3 9 0 98 15 6 0
768553 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 239 0 0 0 0 0 0 0 63 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 239 3 0 0
202829 0 0 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 60 0 0 0 0 0 0 0 0 0 0 0 0 635 0 0 0 0 0 0 0 28 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 396 28 0 0
512485 0 0 0 0 82 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
143362 0 0 6 0 16 5 0 0 10 0 0 4 0 0 0 0 0 12 3 0 0 0 0 0 3 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 13 34 1 0 0 2 0 0 120 0 0 0 3 0 0 15 0 0 0
861881 0 0 0 0 0 547 0 0 2 5 6248 29 1 2 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 33 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
574102 0 7 12 0 2 0 0 0 73 0 0 0 5 0 13 0 0 27 125 18 0 0 24 0 14 1 7 6 0 0 0 0 0 3 0 0 0 0 0 0 2 4 1 6 24 0 2 0 0 38 1 0 0 0 185 7 15 0 3 0
2360704 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 88 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
68621 6 0 0 8 4 0 3 0 31 23 0 17 15 1 6 0 0 27 11 100 0 35 0 79 10 2 1 16 83 9 5 7 26 9 6 47 95 2 7 12 27 10 15 30 20 21 305 11 14 79 4 0 0 10 22 13 24 10 1 0
4318352 0 0 37 0 8 0 0 0 0 0 0 0 0 0 0 2 30 0 0 2 0 0 0 0 0 0 2 4 0 0 0 0 0 0 2 52 0 0 0 4 0 0 1 0 0 0 0 0 116 0 0 0 0 0 8 0 0 0 0 0
4366522 0 3 5 0 92 0 0 3 0 0 0 0 0 1 1 0 0 0 0 2 0 0 19 0 1 0 44 3 0 0 0 0 0 5 0 6 0 0 0 2 0 0 0 0 0 0 0 0 0 25 0 0 0 133 0 0 52 0 1 0
900973 0 105 203 125 67 5 0 0 2 190 0 0 181 0 0 1 23 3 77 3 24 0 99 23 15 46 0 0 0 0 0 0 0 0 60 18 0 71 0 91 1 0 14 2 106 31 378 180 0 0 0 0 0 0 3 16 102 18 294 94
933546 0 57 115 7 0 0 0 0 0 92 30 0 0 0 0 0 0 0 0 0 265 0 66 84 207 1 0 33 0 0 0 0 0 0 0 1 0 0 49 0 0 197 0 0 486 0 890 0 0 0 0 0 0 11 1148 0 113 0 0 0
3393186 0 10 30 0 0 0 0 10 0 0 0 0 0 0 6 0 0 0 55 0 102 0 0 2 0 0 0 0 0 0 0 0 0 18 0 219 0 0 7 0 53 446 0 0 0 0 0 0 1 0 8 0 0 0 172 0 0 0 479 0
386273 0 0 0 18 0 0 79 0 0 0 0 227 0 0 19 2 69 3 105 0 0 0 50 0 7 0 0 0 0 0 0 0 3 62 226 0 0 0 12 675 21 94 0 52 30 0 0 2 0 0 181 0 0 33 6 0 65 5 0 16
4368261 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 85 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
495017 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 113 0 0 0
979707 0 0 0 7 0 0 0 0 0 0 0 81 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 1 0 0 45 0 0 0 0 0 0 0 0 0 0 0 0
760967 23 0 0 15 0 0 0 0 2 0 4 0 0 0 25 2 81 0 321 0 0 0 0 0 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 156 35 0 0 0 0 0 333 0 0 0 0 0 0 0 0 0
1077373 40 0 0 7 18 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 3 0 0 299 0 0 2 0 840 4 0 0 17 1 0 0 3 0 0 0 0 0 0 1 0 0 0 0 34
292041 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 0 0 83 0 0 0 0 0 0
4426163 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 193 0 0 2 0 7 0 0 0 0 0 1 0 0 0 0 2 0 0 0 0 0 0 17 0 0 0 0 0 0
4302472 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 330 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 2 0 0 0
4451252 0 8 18 2 3 0 0 0 0 0 0 0 0 0 0 2 11 0 0 0 0 0 14 0 0 2 3 0 10 282 0 0 0 0 0 121 0 0 0 0 0 6 0 0 1 0 9 0 59 0 0 0 0 13 3 0 0 0 81 0
4307391 0 0 1 0 10 0 0 0 2 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 1 0 0 120 0 0 0 0 0 34 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 4 0 0 0
851668 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 147 0 1 0 0 0 0 0 0 0 0 0 0 0 210 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
970138 0 0 15 0 6 0 0 0 0 0 0 0 0 0 0 0 69 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2 7 0 0 0 2 0 0 0 0 0 0 0 0 0 87 0 0 0 17 0 0 11 0 0 0
4420570 0 0 136 0 8 0 17 4 521 0 1 0 0 0 1 1 0 0 0 0 0 4 0 0 0 26 2 1 0 0 0 10 0 2 9 0 0 0 1 0 13 0 19 0 18 0 31 6 17 161 12 6 0 60 19 0 152 3 14 0
503315 74 362 167 29 158 342 624 458 0 599 8 609 1 0 364 286 208 419 365 0 0 0 54 29 322 19 12 25 0 21 83 10 0 659 897 44 0 38 84 952 72 146 1 818 176 0 2 21 30 101 246 6 44 141 180 0 687 15 206 534
504674 74 84 18 4 112 14 271 0 0 169 0 99 0 0 155 61 21 105 28 16 0 0 6 2 190 0 1 17 29 0 0 0 0 59 66 6 89 31 9 13 3 37 17 117 0 3 0 0 73 163 13 3 0 148 180 23 28 8 10 586
937813 11 12 5 0 23 9 14 0 0 52 0 19 0 0 29 1 10 33 7 8 0 0 2 1 52 0 0 0 5 0 0 0 0 26 24 0 11 0 2 2 0 8 5 25 1 0 0 0 19 0 4 0 0 28 37 3 7 0 5 124
224222 0 0 0 0 5 0 3 2 183 0 0 3 0 0 1 2 4 0 3 0 0 0 1 4 3 0 1 0 0 0 0 0 0 3 2 0 0 0 0 0 20 1 0 1 16 0 62 0 1 0 0 0 0 2 1 0 7 0 1 0
451449 0 22 270 11 16 0 10 10 0 2 0 144 1 1 146 5 0 152 140 1 25 0 4 0 12 0 2 12 0 0 0 0 0 12 3 0 0 0 4 233 1 1 2 381 1 0 4 0 25 0 0 0 0 3 251 0 0 0 42 8
2751958 269 156 0 5 192 139 679 534 1 26 3 334 1 0 116 340 8 102 274 14 22 0 287 0 600 1 2 4 0 4 0 7 0 929 370 10 67 15 62 8 13 167 37 122 0 0 0 27 17 204 22 3 0 257 228 0 96 25 143 514
403258 0 23 0 1 52 39 0 0 0 0 0 371 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1029165 0 0 39 20 11 2 117 0 0 0 0 0 0 0 33 0 13 33 38 0 0 0 29 0 37 0 0 0 0 0 0 0 0 0 193 0 0 0 4 24 0 3 0 78 0 0 0 26 0 0 72 0 0 0 0 0 54 0 8 142
207936 0 0 0 0 8 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 138
494906 92 280 358 60 205 247 1009 429 415 871 10 1929 0 9 593 172 331 885 1262 0 0 0 189 0 560 0 0 1 0 0 0 21 10 807 1059 38 0 4 61 490 1 305 20 1251 3 0 0 126 0 368 223 0 20 0 393 0 511 109 0 728
183158 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 4 0 0 1 0 0 0 0 1 0 0 1 10 0 273 0 0 0 0 0 0 0 17 0 0 43 3 0 4 92 0 16 0 0 0 0 0 0 13 4 0 0 0 0 0
4374753 0 0 3 0 1 0 0 1 69 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 107 0 0 0 0 2 37 0 0 0 0 1 1 0 0 0 0 9 0 23 0 0 1 0 0 1 3 0 0 7 0
4475758 0 5 9 0 7 0 0 0 7 0 0 0 0 0 0 9 6 0 0 0 0 0 0 0 1 0 1 0 0 9 0 0 0 3 0 7 0 0 0 10 0 0 0 0 0 0 2 0 31 27 414 3 0 97 5 0 26 0 8 0
4296244 0 0 9 0 1 0 0 4 171 0 0 0 0 0 1 21 0 0 0 0 0 2 0 0 0 0 10 0 0 141 0 0 0 0 3 90 0 0 0 0 2 2 0 0 0 0 18 0 43 0 5 0 0 0 0 3 0 0 12 0
1143045 0 8 13 3 24 0 0 0 11 0 0 1 0 0 0 27 10 0 0 0 0 0 0 0 9 1 25 2 0 17 0 0 0 5 0 15 0 0 0 18 0 1 0 0 0 0 2 0 115 57 835 7 0 161 4 0 39 0 15 0
4310208 0 2 7 0 8 0 0 6 7 0 0 0 0 0 0 11 1 0 0 0 0 0 0 0 2 0 8 1 0 0 0 0 0 3 0 32 0 0 0 20 0 0 0 0 0 0 0 5 150 128 327 6 0 171 31 0 35 0 42 0
4422456 0 0 8 0 27 0 0 26 5 0 0 0 0 0 3 4 2 0 0 0 0 0 0 0 0 0 0 0 15 0 0 0 0 0 1 24 0 0 0 7 0 0 0 0 0 0 0 3 57 68 187 1 0 124 15 0 30 0 27 0
3801267 23 0 39 0 3 0 0 0 0 9 1 0 0 0 4 0 259 0 0 0 0 0 0 0 0 0 0 0 0 38 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 49 0 0 3 0 0 0 0 3 0 7 0 0 0
128382 0 0 0 7 0 0 48 61 0 0 0 211 0 0 1 0 15 0 60 0 0 0 0 0 25 0 0 0 0 0 0 0 0 2 154 0 0 0 13 158 0 0 0 34 0 0 0 42 0 0 0 0 0 0 26 0 13 0 0 26
4318284 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 117 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
753638 0 0 3 1 10 0 0 5 12 0 0 0 0 0 5 4 0 0 32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 102 5 0 2 29 0 0 0 0 0 0 0 0 0 0 0 0 20 5 0 30
496787 0 15 12 0 0 0 21 186 0 5 3 4 207 270 109 26 64 0 48 64 8 0 24 16 12 2 4 1476 0 0 0 10 0 3 10 34 0 40 13 25 10 48 6 32 28 1 0 11 0 84 5 6 17 8 0 0 11 0 17 38
377874 0 8 5 0 22 0 45 10 17 0 1 34 21 1 466 10 212 12 81 2 10 0 51 15 12 29 57 2 15 9 0 7 10 5 14 1 112 2 23 8 427 9 45 48 50 1 27 6 1 210 14 36 30 74 223 0 0 0 47 12
115186 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 123 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
934488 263 35 111 0 356 53 209 114 24 117 9 61 84 38 187 321 42 239 131 49 3 2 184 7 248 74 60 26 93 0 0 24 10 126 95 180 362 15 1 17 195 27 33 332 27 5 142 0 36 84 8 14 0 45 65 0 13 25 233 50
356733 2789 454 1125 8 4284 467 1880 1354 198 1075 118 763 1145 512 2111 3761 352 1315 1026 382 39 12 2028 172 1949 555 646 288 771 4 0 283 234 950 533 1188 2359 90 4 79 1191 143 230 2219 136 42 1016 11 322 992 54 196 98 545 770 72 252 188 2926 632
868615 0 0 3 0 0 5 0 0 0 0 1 890 0 1 0 4 0 0 0 4919 3162 3358 153 4517 23 36 47 3036 0 2943 398 5997 6582 0 1 3 1757 5183 464 323 2 0 1 26 5021 0 0 5157 0 3 0 1 6434 312 0 7095 0 0 8 1380
517754 6 2 41 0 12 0 0 0 102 0 0 0 0 2 3 6 170 0 0 0 0 0 0 0 0 0 39 0 0 0 0 0 0 0 0 25 0 0 0 4 0 1 2 0 0 0 24 0 119 57 0 0 0 0 0 0 13 0 0 0
561636 0 0 2 0 0 0 0 0 38 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 3 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 5 0 0 0 0 410 0 0 2 0 0 0
965500 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 1 2 0 2 0 11 1 0 0 5 0 0 148 21 0 9 0 0 0 0 99 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0
871442 0 109 696 0 52 0 0 4 342 14 17 14 137 11 19 177 251 21 3 452 23 186 33 1077 117 531 389 238 49 700 5022 684 16 426 26 82 56 17 3661 28 17 32 13 16 28 2 87 0 982 180 1 200 7 394 31 20 391 0 208 4
4439603 0 0 18 0 3 0 0 2 451 0 0 2 2 0 1 12 4 0 3 1 2 0 0 0 0 0 17 0 0 128 0 0 0 0 1 0 0 0 0 0 9 0 0 7 9 0 2 0 211 0 57 14 0 4 1 0 0 0 0 0
2901965 0 3 12 0 0 2 0 1 143 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 4 0 0 77 0 0 7 0 1 6 11 0 0 0 12 12 0 1 3 0 0 0 22 3 4 2 0 2 4 0 0 5 10 0
579608 0 5 23 0 170 0 0 10 40 0 0 2 0 0 4 52 87 0 0 0 11 0 0 0 0 0 0 0 0 30 1 0 0 0 1 0 6 4 1 4 0 18 6 0 4 0 125 0 0 82 0 0 10 31 28 0 65 0 3 0
516966 0 3 15 0 0 2 0 1 161 0 0 1 0 0 1 1 0 0 0 1 0 0 0 0 0 2 9 0 0 111 0 0 0 0 1 9 17 0 0 5 16 16 0 1 4 0 0 0 53 5 5 0 0 6 10 0 0 18 10 0
967427 0 0 13 0 1 0 0 0 11 0 2 1 4 0 7 4 16 0 0 8 0 1 0 1 0 1 18 1 0 0 0 0 0 0 4 0 0 0 0 17 0 5 2 3 0 0 58 0 1 49 4 0 3 3 20 0 122 0 42 0
172804 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 220 0 0 10 0 0 0 0 0 0
557570 0 0 10 0 90 0 0 0 24 0 0 0 0 0 0 18 1 6 0 0 11 0 0 0 0 1 4 0 10 786 0 0 0 0 1 3 0 0 0 2 2 5 3 0 0 0 7 0 0 0 0 0 0 1 0 16 0 0 84 0
291149 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 72 78 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0
301921 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 216 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0
1995182 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 333 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1076668 0 0 0 0 21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 232 0 0 0 0 0 0 7 0 4 0
885628 0 0 7 0 0 0 17 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 26 2 0 0 0 0 0 0 0 0 0 699 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0
882921 6 0 0 5990 3 381 185 0 0 1258 805 0 126 482 0 0 820 1037 1 33 6 13 229 159 437 3151 3 2 0 0 0 0 0 0 0 2485 162 0 215 0 70 0 3796 0 0 2850 2 0 678 5 2909 2747 7 2 125 3 1903 3915 0 2
590982 0 2 0 0 1 0 0 62 1 52 0 0 9 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 0 0 0 24 0 0 0 0 0 3 21 4 0 10 82 0 0 0 0 0 0 7 0 0 0
4331684 0 4321 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
196428 0 0 0 0 0 0 0 0 162 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0
366195 0 0 0 0 0 0 0 0 103 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
328105 0 0 0 0 2 0 0 0 82 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 1 0 0 0 9 0 5 1 0 9 8 0 7 0 0 0
227000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 128 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
347382 0 0 0 0 0 0 0 0 100 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 9 0 0 0 0 0
386088 17 7 36 1 33 9 0 95 94 66 0 90 0 3 38 15 9 18 39 0 2 0 7 0 27 2 0 0 0 0 0 3 0 12 21 7 0 8 1 1 0 4 2 43 7 0 0 2 7 27 1 2 0 61 36 0 0 5 1 28
403853 46 248 117 3 11 18 182 605 389 324 5 35 43 76 18 52 64 206 209 12 0 0 6 1 6 2 1 2 10 0 0 10 3 99 60 80 28 133 16 103 1 14 21 5 5 4 49 14 39 33 14 3 30 45 84 40 150 5 38 24
4294749 0 0 3 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 3 0 0 0 0 1 6 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 3 0 87 0 0 0 0 23 0 11 0 0 2 0 1 0
903426 0 0 4 0 23 0 0 0 6 0 0 1 0 0 3 5 180 0 0 4 4 0 0 2 0 1 3642 39 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 5 0 0 0 2 0 0 0 0 5 0
866280 0 0 5 0 1 0 0 0 125 0 1 1 0 0 7 4 1 0 0 4 0 1 0 0 0 0 23 0 0 26 0 0 0 0 3 9 0 0 0 2 0 0 0 4 0 0 11 0 35 0 0 0 0 0 0 0 2 0 0 0
4466006 6 0 2 0 4 0 0 5 103 0 0 0 7 0 0 2 4 3 0 14 5 0 20 0 0 1 63 0 0 9 0 0 0 3 0 15 0 0 0 8 0 2 10 0 1 0 105 0 0 0 3 17 13 12 8 0 26 0 1 0
1017181 0 0 1 0 1 0 0 0 37 0 0 0 0 0 0 0 5 0 0 0 0 1 0 0 0 0 77 1 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 27 0 32 0 0 0 0 0 1 0 15 0 0 0
930926 0 0 0 616 608 857 0 0 0 977 78 0 3132 1081 202 0 1913 3 0 31 0 1 320 825 810 1549 2 1 3908 0 0 0 0 5 0 1238 72 720 133 0 1804 392 2253 0 0 4279 0 2 0 0 120 2657 0 0 0 0 0 1209 0 0
974249 0 0 0 0 1 0 0 0 0 0 1 0 0 0 4 1 1 0 1 0 0 0 17 0 1 5 5 2 0 0 0 7 0 0 4 0 0 0 0 0 41 5 0 4 11 0 22 0 20 38 2 5 27 3 186 0 2 0 0 2
4296075 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 295 0 0 0 0 0
128390 0 0 0 3 0 5 0 0 0 0 0 3 0 0 0 0 0 0 11 0 0 0 0 0 0 0 0 29 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 4 4 0 0 0 409 0 0 0 0 0 8 0 0 0 0 0
552026 0 0 0 0 0 0 0 0 11 0 0 0 0 4 1 0 8 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 8 0 0 0 0 20 0 166 0 2 0 0 1 47 0 0 0 0 0
28246 0 0 0 0 0 0 0 0 0 0 0 0 209 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 0 425 0 0 0 0 11 0 0 0 0 0 0 0 0
544480 0 0 0 0 0 0 0 26 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 2 1 0 20 0 134 0 0 0 0 0 19 0 0 0 2 0
912997 0 0 1 1 0 0 0 0 0 0 0 0 0 0 2 0 128 0 0 0 0 3 0 0 0 0 0 24 0 0 25 0 0 0 0 0 0 0 4 0 0 133 0 0 0 0 371 0 0 30 30 0 0 0 0 0 0 0 3 0
495096 0 0 0 0 0 0 0 0 4 0 0 9 0 1 0 2 0 6 0 0 899 0 0 1 1 0 2 61 0 0 0 0 0 0 0 0 0 0 0 0 1 6 0 0 4 0 0 0 8 0 0 0 0 6 17 0 0 0 0 0
446403 17 47 80 0 0 2 0 0 9 0 0 0 0 0 1 1 0 0 52 0 5 0 0 0 0 5 1 1173 0 0 0 10 0 0 0 109 0 8 0 0 22 28 0 0 1 78 0 0 0 0 0 76 27 0 0 0 0 0 45 0
1047041 6 0 215 0 0 83 889 0 2 94 25 1 0 160 16 214 4 36 42 70 292 0 188 0 9 0 1 0 0 0 0 0 0 0 0 1 0 273 0 247 0 3 2 29 0 0 0 6 0 0 28 8 47 12 8 0 13 23 0 90
505626 0 0 0 24 8 0 0 98 0 0 0 0 4 0 130 0 0 161 21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 0 104 0 0 1 124 0 0 0 0 0 0 30 0 0 0 0 0 0 0 0 8
858026 6 0 0 0 0 0 0 17 0 0 0 0 1 1 1 24 0 0 29 0 0 0 3 0 3 0 1 0 0 0 0 0 0 6 91 0 6 0 0 7 32 1 1 0 0 0 0 0 90 0 0 0 0 2 1 0 2 0 0 0
467198 23 13 5 0 11 25 0 190 0 0 3 2 31 4 13 268 4 0 271 1 0 0 37 0 29 1 14 2 0 0 0 0 0 18 702 0 117 0 1 79 225 4 4 2 14 0 7 0 1461 25 1 0 0 7 18 0 7 38 2 2
861807 29 137 168 10 62 0 124 1203 46 1096 17 871 58 440 1239 286 197 0 625 75 8 5 802 140 812 3 86 385 2171 0 0 66 16 2938 106 15 1500 585 217 133 136 44 11 798 117 0 1682 31 1104 545 63 9 162 67 1082 7 17 3 28 144
377613 2714 680 2322 403 2 4116 902 1180 1210 16 28 398 1789 4296 1237 982 1069 2180 1289 754 2161 4 1863 4 663 1283 1674 275 24 9 1 0 0 6 2224 3 11 25 1876 1955 2296 4351 566 5 5 2 11 1482 9 1992 700 939 229 2590 546 0 555 1642 1964 1896 
# Calculate the relative abundance for each sample
ps.ng.tax_abund_rel <- transform_sample_counts(ps.ng.tax_abund, function(x) x/sum(x))
datamat_ = as.data.frame(otu_table(ps.ng.tax_abund))
# datamat <- datamat_[c('8','9','10','12','13','14',
# '21','22','23','24','25','26','27','28',
# '33','34','35','36','37','38','39','51', '47','48','49','50','52','53','55')]
#--take all columns--
datamat <- datamat_
hr <- hclust(as.dist(1 - cor(t(datamat), method = "pearson")), method = "complete")
hc <- hclust(as.dist(1 - cor(datamat, method = "spearman")), method = "complete")
mycl = cutree(hr, h = max(hr$height)/1.08)
mycol = c("YELLOW", "DARKBLUE", "DARKORANGE", "DARKMAGENTA", "DARKCYAN", "DARKRED",
    "MAROON", "DARKGREEN", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN", "LIGHTGREEN",
    "BLUE", "ORANGE", "CYAN", "RED", "GREEN")
mycol = mycol[as.vector(mycl)]
sampleCols <- rep("GREY", ncol(datamat))
# Colors: '#a6cee3', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99', '#e31a1c',
# '#cab2d6', '#6a3d9a' Control (19)
sampleCols[colnames(datamat) == "kg001"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg002"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg003"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg004"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg005"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg007"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg009"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg015"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg016"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg019"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg020"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg021"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg022"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg023"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg025"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg026"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg027"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg028"] <- "#1f78b4"
sampleCols[colnames(datamat) == "kg029"] <- "#1f78b4"
# GPA (21)
sampleCols[colnames(datamat) == "micro1"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro3"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro4"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro6"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro7"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro12"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro13"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro16"] <- "#33a02c"
sampleCols[colnames(datamat) == "micro17"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw001"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw004"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw005"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw006"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw007"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw009"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw010"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw013"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw014"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw015"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw017"] <- "#33a02c"
sampleCols[colnames(datamat) == "mw021"] <- "#33a02c"
# RA (20)
sampleCols[colnames(datamat) == "ra002"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra003"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra004"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra005"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra006"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra007"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra008"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra009"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra010"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra013"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra014"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra015"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra017"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra018"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra019"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra020"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra022"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra023"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra024"] <- "#e31a1c"
sampleCols[colnames(datamat) == "ra025"] <- "#e31a1c"

library(RColorBrewer)
custom_palette <- colorRampPalette(brewer.pal(9, "Blues"))
heatmap_colors <- custom_palette(100)
png("figures/heatmap.png", width = 1200, height = 2400)
heatmap.2(as.matrix(datamat), Rowv = as.dendrogram(hr), Colv = NA, dendrogram = "row",
    scale = "row", trace = "none", col = heatmap_colors, cexRow = 1.2, cexCol = 1.5,
    RowSideColors = mycol, ColSideColors = sampleCols, srtCol = 15, labRow = row.names(datamat),
    key = TRUE, margins = c(10, 15), lhei = c(0.7, 15), lwid = c(1, 8))
dev.off()
png 
  2 
knitr::include_graphics("./figures/heatmap.png")
Heatmap

Heatmap

Taxonomic summary

Bar plots in phylum level

library(ggplot2)
geom.text.size = 6
theme.size = 8  #(14/5) * geom.text.size
ps.ng.tax_most = phyloseq::filter_taxa(ps.ng.tax_rel, function(x) mean(x) > 0.001,
    TRUE)
ps.ng.tax_most_ = transform_sample_counts(ps.ng.tax_most, function(x) x/sum(x))
plot_bar(ps.ng.tax_most_, fill = "Phylum") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

ps.ng.tax_most_copied <- data.table::copy(ps.ng.tax_most_)

Bar plots in class level

plot_bar(ps.ng.tax_most_copied, fill = "Class") + geom_bar(aes(), stat = "identity",
    position = "stack") + scale_fill_manual(values = c("darkblue", "darkgoldenrod1",
    "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen",
    "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4",
    "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue",
    "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1",
    "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Bar plots in order level

plot_bar(ps.ng.tax_most_copied, fill = "Order") + geom_bar(aes(), stat = "identity",
    position = "stack") + scale_fill_manual(values = c("darkblue", "darkgoldenrod1",
    "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen",
    "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4",
    "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue",
    "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1",
    "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 4))

Bar plots in family level

plot_bar(ps.ng.tax_most_copied, fill = "Family") + geom_bar(aes(), stat = "identity",
    position = "stack") + scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF",
    "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080",
    "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000",
    "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF",
    "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9",
    "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017",
    "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117",
    "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117",
    "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = 5, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 8))

Alpha diversity

Plot Chao1 richness estimator, Observed OTUs, Shannon index, and Phylogenetic diversity. Regroup together samples from the same group.

hmp.div_qiime <- read.csv("adiv_even.txt", sep = "\t")
colnames(hmp.div_qiime) <- c("sam_name", "chao1", "observed_otus", "shannon", "PD_whole_tree")
row.names(hmp.div_qiime) <- hmp.div_qiime$sam_name
div.df <- merge(hmp.div_qiime, hmp.meta, by = "sam_name")
div.df2 <- div.df[, c("SampleType", "chao1", "shannon", "observed_otus", "PD_whole_tree")]
colnames(div.df2) <- c("Group", "Chao-1", "Shannon", "OTU", "Phylogenetic Diversity")
options(max.print = 999999)
knitr::kable(div.df2) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Group Chao-1 Shannon OTU Phylogenetic Diversity
control 38.00000 2.5126324 31 4.65148
control 82.11111 2.6519845 67 7.86756
control 199.66667 4.1349045 160 11.50759
control 42.00000 1.2463443 40 5.86806
control 188.68182 3.0985879 155 12.05980
control 61.42857 2.4456359 42 3.94668
control 39.50000 3.4289966 36 4.13313
control 88.66667 3.8060304 80 8.20430
control 142.60000 5.5176171 139 12.43166
control 41.60000 3.6188130 41 5.36125
control 95.00000 0.9732663 72 7.70470
control 109.23529 3.6748240 106 9.15602
control 62.75000 2.6264060 49 4.90319
control 69.27273 2.1585520 52 4.21946
control 127.58824 3.4188787 102 9.61382
control 148.00000 2.9950552 129 9.87727
control 132.60000 4.1923767 119 10.65810
control 73.11111 3.3864724 58 6.59343
control 120.40000 4.2420857 105 9.31995
GPA 129.76923 2.2335361 112 9.96287
GPA 147.20000 3.7174971 143 10.92196
GPA 156.23077 2.4347508 125 10.36394
GPA 225.96875 2.7675173 204 12.85912
GPA 147.75000 2.8307342 109 9.46643
GPA 183.33333 2.5857569 150 10.94254
GPA 106.50000 1.6027140 81 7.92421
GPA 153.37500 3.7046056 132 8.88776
GPA 101.61538 2.1616398 87 7.56155
GPA 32.20000 1.9577693 28 4.72075
GPA 62.00000 3.5884207 57 6.42758
GPA 44.00000 1.9468255 35 4.55259
GPA 37.62500 1.2476340 32 4.69088
GPA 25.66667 0.8671918 24 3.64861
GPA 70.12500 2.9601222 66 6.98732
GPA 103.30000 3.5152901 88 8.74152
GPA 138.00000 3.5230938 111 9.42402
GPA 39.66667 2.9541263 35 3.37007
GPA 43.00000 1.8556788 40 4.73671
GPA 84.00000 2.5440048 73 6.83341
GPA 103.68750 3.9434930 93 8.57517
RA 121.00000 3.3449784 115 9.61511
RA 191.00000 2.9762825 184 13.02167
RA 107.00000 2.1913073 97 8.83295
RA 173.43750 3.7432315 170 13.24291
RA 152.06667 2.6193442 125 11.40126
RA 49.42857 1.3489865 40 5.83268
RA 113.66667 4.3555040 91 8.86575
RA 48.00000 1.6248900 44 5.91250
RA 93.60000 4.2709559 88 6.78461
RA 84.00000 4.5113078 80 7.68684
RA 111.68750 3.5108714 106 9.23421
RA 115.00000 2.5577960 100 9.89423
RA 41.50000 1.1329113 38 4.23374
RA 177.06667 4.5805528 168 13.79084
RA 148.15385 4.7489461 135 11.36418
RA 29.66667 0.4577077 25 4.35241
RA 142.15385 4.6629967 126 11.02192
RA 39.33333 2.0711502 37 5.38240
RA 111.57143 3.1565146 106 9.39170
RA 68.33333 3.5058156 55 6.08455 
# https://uc-r.github.io/t_test We can perform the test with t.test and
# transform our data and we can also perform the nonparametric test with the
# wilcox.test function.
stat.test.Shannon <- compare_means(Shannon ~ Group, data = div.df2, method = "t.test")
knitr::kable(stat.test.Shannon) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
.y. group1 group2 p p.adj p.format p.signif method
Shannon control GPA 0.0857343 0.26 0.086 ns T-test
Shannon control RA 0.8014361 0.80 0.801 ns T-test
Shannon GPA RA 0.1979828 0.40 0.198 ns T-test 
div_df_melt <- reshape2::melt(div.df2)
# head(div_df_melt) https://plot.ly/r/box-plots/#horizontal-boxplot
# http://www.sthda.com/english/wiki/print.php?id=177
# https://rpkgs.datanovia.com/ggpubr/reference/as_ggplot.html
# http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/82-ggplot2-easy-way-to-change-graphical-parameters/
# https://plot.ly/r/box-plots/#horizontal-boxplot library('gridExtra')
# par(mfrow=c(4,1))
p <- ggboxplot(div_df_melt, x = "Group", y = "value", facet.by = "variable", scales = "free",
    width = 0.5, fill = "gray", legend = "right")
# ggpar(p, xlab = FALSE, ylab = FALSE)

lev <- levels(factor(div_df_melt$Group))  # get the variables
# Using my_stat_compare_means (see
# https://stackoverflow.com/questions/47839988/indicating-significance-with-ggplot2-in-a-boxplot-with-multiple-groups)
L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i])  #%>% filter(p.signif != 'ns')
my_stat_compare_means <- function(mapping = NULL, data = NULL, method = NULL, paired = FALSE,
    method.args = list(), ref.group = NULL, comparisons = NULL, hide.ns = FALSE,
    label.sep = ", ", label = NULL, label.x.npc = "left", label.y.npc = "top", label.x = NULL,
    label.y = NULL, tip.length = 0.03, symnum.args = list(), geom = "text", position = "identity",
    na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, ...) {
    if (!is.null(comparisons)) {
        method.info <- ggpubr:::.method_info(method)
        method <- method.info$method
        method.args <- ggpubr:::.add_item(method.args, paired = paired)
        if (method == "wilcox.test")
            method.args$exact <- FALSE
        pms <- list(...)
        size <- ifelse(is.null(pms$size), 0.3, pms$size)
        color <- ifelse(is.null(pms$color), "black", pms$color)
        map_signif_level <- FALSE
        if (is.null(label))
            label <- "p.format"
        if (ggpubr:::.is_p.signif_in_mapping(mapping) | (label %in% "p.signif")) {
            if (ggpubr:::.is_empty(symnum.args)) {
                map_signif_level <- c(`****` = 1e-04, `***` = 0.001, `**` = 0.01,
                  `*` = 0.05, ns = 1)
            } else {
                map_signif_level <- symnum.args
            }
            if (hide.ns)
                names(map_signif_level)[5] <- " "
        }
        step_increase <- ifelse(is.null(label.y), 0.12, 0)
        ggsignif::geom_signif(comparisons = comparisons, y_position = label.y, test = method,
            test.args = method.args, step_increase = step_increase, size = size,
            color = color, map_signif_level = map_signif_level, tip_length = tip.length,
            data = data)
    } else {
        mapping <- ggpubr:::.update_mapping(mapping, label)
        layer(stat = StatCompareMeans, data = data, mapping = mapping, geom = geom,
            position = position, show.legend = show.legend, inherit.aes = inherit.aes,
            params = list(label.x.npc = label.x.npc, label.y.npc = label.y.npc, label.x = label.x,
                label.y = label.y, label.sep = label.sep, method = method, method.args = method.args,
                paired = paired, ref.group = ref.group, symnum.args = symnum.args,
                hide.ns = hide.ns, na.rm = na.rm, ...))
    }
}

# FITTING_3: #> Using lev ('control', 'GPA', 'RA') for definition of the
# following comparisons. comparisons = L.pairs, symnum.args <- list(cutpoints =
# c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c('****', '***', '**', '*')),
p3 <- p + stat_compare_means(method = "t.test", comparisons = list(c("GPA", "control"),
    c("RA", "control"), c("GPA", "RA")), label = "p.signif", symnum.args <- list(cutpoints = c(0,
    1e-04, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")))
# stat_pvalue_manual print(p2)
# https://stackoverflow.com/questions/20500706/saving-multiple-ggplots-from-ls-into-one-and-separate-files-in-r
# FITTING_4: mkdir figures
ggsave("./figures/alpha_diversity_SampleType.png", device = "png", height = 10, width = 12)
ggsave("./figures/alpha_diversity_SampleType.svg", device = "svg", height = 10, width = 12)

Selected alpha diversity

knitr::include_graphics("./figures/alpha_diversity_SampleType.png")
Alpha diversity

Alpha diversity 

# FITTING_5: Preparing the file selected_alpha_diversities.txt
# selected_alpha_diversities<-read.csv('selected_alpha_diversities.txt',sep='\t')
# knitr::kable(selected_alpha_diversities) %>% kable_styling(bootstrap_options
# = c('striped', 'hover', 'condensed', 'responsive'))

Beta diversity (weighted_unifrac)

# The tests of significance were performed using a two-sided Student's
# two-sample t-test.  Alternative hypothesis: Group 1 mean != Group 2 mean The
# nonparametric p-values were calculated using 999 Monte Carlo permutations.
# The nonparametric p-values contain the correct number of significant digits.
# Entries marked with 'N/A' could not be calculated because at least one of the
# groups of distances was empty, both groups each contained only a single
# distance, or the test could not be performed (e.g. no variance in groups with
# the same mean).  PREPARING_FILES: cp
# bdiv_even1300_SampleType/weighted_unifrac_boxplots/SampleType_Stats.txt
# bdiv_even1300_SampleType/weighted_unifrac_boxplots/SampleType_Stats_.txt and
# delete all headers, add # before 'Group 1'.  cp
# bdiv_even1300_SampleType/unweighted_unifrac_boxplots/SampleType_Stats.txt
# bdiv_even1300_SampleType/unweighted_unifrac_boxplots/SampleType_Stats_.txt
# and delete all headers, add # before 'Group 1'.
# file:///home/jhuang/DATA/Data_Marius_16S/core_diversity_e42369/bdiv_even42369_Group/unweighted_unifrac_boxplots/Group_Stats.txt
weighted_unifrac_stats <- read.csv(file = "bdiv_even1300_SampleType/weighted_unifrac_boxplots/SampleType_Stats_.txt",
    sep = "\t")
knitr::kable(weighted_unifrac_stats) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
X.Group.1 Group.2 t.statistic Parametric.p.value Parametric.p.value..Bonferroni.corrected. Nonparametric.p.value Nonparametric.p.value..Bonferroni.corrected.
All within SampleType All between SampleType -0.9448843 0.3448471 1.0000000 0.336 1.000
All within SampleType control vs. control 0.0742240 0.9408522 1.0000000 0.945 1.000
All within SampleType GPA vs. GPA -1.2065252 0.2279813 1.0000000 0.244 1.000
All within SampleType RA vs. RA 1.3214246 0.1867580 1.0000000 0.178 1.000
All within SampleType control vs. GPA -2.4760610 0.0134542 0.3767165 0.007 0.196
All within SampleType control vs. RA 0.6098859 0.5420834 1.0000000 0.542 1.000
All within SampleType GPA vs. RA -0.1972929 0.8436389 1.0000000 0.834 1.000
All between SampleType control vs. control 0.6629703 0.5074612 1.0000000 0.521 1.000
All between SampleType GPA vs. GPA -0.7058209 0.4804164 1.0000000 0.472 1.000
All between SampleType RA vs. RA 2.0469735 0.0408488 1.0000000 0.031 0.868
All between SampleType control vs. GPA -2.0543775 0.0401017 1.0000000 0.040 1.000
All between SampleType control vs. RA 1.5251598 0.1274198 1.0000000 0.135 1.000
All between SampleType GPA vs. RA 0.6349431 0.5255554 1.0000000 0.535 1.000
control vs. control GPA vs. GPA -0.9373706 0.3491649 1.0000000 0.347 1.000
control vs. control RA vs. RA 0.9865003 0.3245521 1.0000000 0.325 1.000
control vs. control control vs. GPA -1.7467733 0.0812172 1.0000000 0.069 1.000
control vs. control control vs. RA 0.3592265 0.7195638 1.0000000 0.740 1.000
control vs. control GPA vs. RA -0.2102543 0.8335419 1.0000000 0.824 1.000
GPA vs. GPA RA vs. RA 2.1631960 0.0311208 0.8713822 0.039 1.000
GPA vs. GPA control vs. GPA -0.7225328 0.4702453 1.0000000 0.468 1.000
GPA vs. GPA control vs. RA 1.6587818 0.0976930 1.0000000 0.103 1.000
GPA vs. GPA GPA vs. RA 1.0373044 0.2999933 1.0000000 0.304 1.000
RA vs. RA control vs. GPA -3.1483258 0.0017255 0.0483127 0.002 0.056
RA vs. RA control vs. RA -0.8667216 0.3864605 1.0000000 0.389 1.000
RA vs. RA GPA vs. RA -1.5020995 0.1335906 1.0000000 0.130 1.000
control vs. GPA control vs. RA 2.8913391 0.0039431 0.1104061 0.006 0.168
control vs. GPA GPA vs. RA 2.1828773 0.0293288 0.8212076 0.033 0.924
control vs. RA GPA vs. RA -0.7821717 0.4343457 1.0000000 0.418 1.000 
# BUG:How to generate the plot?
# knitr::include_graphics('./figures/PCoA_weighted_unifrac.png')

Beta diversity (unweighted_unifrac)

# file:///home/jhuang/DATA/Data_Marius_16S/core_diversity_e42369/bdiv_even42369_Group/unweighted_unifrac_boxplots/Group_Stats.txt
unweighted_unifrac_stats <- read.csv(file = "bdiv_even1300_SampleType/unweighted_unifrac_boxplots/SampleType_Stats_.txt",
    sep = "\t")
knitr::kable(unweighted_unifrac_stats) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
X.Group.1 Group.2 t.statistic Parametric.p.value Parametric.p.value..Bonferroni.corrected. Nonparametric.p.value Nonparametric.p.value..Bonferroni.corrected.
All within SampleType All between SampleType -1.9887655 0.0468810 1.0000000 0.059 1.000
All within SampleType control vs. control 0.8077441 0.4194975 1.0000000 0.387 1.000
All within SampleType GPA vs. GPA -0.5117013 0.6090051 1.0000000 0.583 1.000
All within SampleType RA vs. RA -0.2380432 0.8119118 1.0000000 0.824 1.000
All within SampleType control vs. GPA -0.8835497 0.3771587 1.0000000 0.388 1.000
All within SampleType control vs. RA -0.5284031 0.5973431 1.0000000 0.609 1.000
All within SampleType GPA vs. RA -3.2013911 0.0014113 0.0395154 0.004 0.112
All between SampleType control vs. control 2.0875347 0.0370245 1.0000000 0.039 1.000
All between SampleType GPA vs. GPA 0.8135408 0.4160457 1.0000000 0.426 1.000
All between SampleType RA vs. RA 1.0470685 0.2952504 1.0000000 0.295 1.000
All between SampleType control vs. GPA 0.7234669 0.4694991 1.0000000 0.474 1.000
All between SampleType control vs. RA 1.1250900 0.2607221 1.0000000 0.287 1.000
All between SampleType GPA vs. RA -1.7868299 0.0741523 1.0000000 0.093 1.000
control vs. control GPA vs. GPA -1.0747519 0.2831696 1.0000000 0.279 1.000
control vs. control RA vs. RA -0.8469234 0.3976021 1.0000000 0.379 1.000
control vs. control control vs. GPA -1.3439725 0.1794938 1.0000000 0.193 1.000
control vs. control control vs. RA -1.1275912 0.2599852 1.0000000 0.270 1.000
control vs. control GPA vs. RA -3.0258939 0.0025872 0.0724408 0.003 0.084
GPA vs. GPA RA vs. RA 0.2193061 0.8265240 1.0000000 0.823 1.000
GPA vs. GPA control vs. GPA -0.2089473 0.8345595 1.0000000 0.831 1.000
GPA vs. GPA control vs. RA 0.0685205 0.9453946 1.0000000 0.949 1.000
GPA vs. GPA GPA vs. RA -1.9986350 0.0460785 1.0000000 0.046 1.000
RA vs. RA control vs. GPA -0.4336156 0.6647269 1.0000000 0.649 1.000
RA vs. RA control vs. RA -0.1734713 0.8623428 1.0000000 0.858 1.000
RA vs. RA GPA vs. RA -2.1755872 0.0299706 0.8391767 0.029 0.812
control vs. GPA control vs. RA 0.3251898 0.7451250 1.0000000 0.757 1.000
control vs. GPA GPA vs. RA -2.0376445 0.0419072 1.0000000 0.050 1.000
control vs. RA GPA vs. RA -2.4146038 0.0159763 0.4473367 0.022 0.616 
# BUG:How to generate the plot?
# knitr::include_graphics('./figures/PCoA_unweighted_unifrac.png')

Differential abundance analysis

Differential abundance analysis aims to find the differences in the abundance of each taxa between two groups of samples, assigning a significance value to each comparison. ## GPA vs control

library("DESeq2")
# ALTERNATIVE using ps.ng.tax_most_copied: ps.ng.tax (40594) vs.
# ps.ng.tax_most_copied (166)
ps.ng.tax_sel <- ps.ng.tax
# FITTING5: correct the id of the group members, see FITTING6
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("micro1", "micro3", "micro4",
    "micro6", "micro7", "micro12", "micro13", "micro16", "micro17", "mw001", "mw004",
    "mw005", "mw006", "mw007", "mw009", "mw010", "mw013", "mw014", "mw015", "mw017",
    "mw021", "kg001", "kg002", "kg003", "kg004", "kg005", "kg007", "kg009", "kg015",
    "kg016", "kg019", "kg020", "kg021", "kg022", "kg023", "kg025", "kg026", "kg027",
    "kg028", "kg029")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~SampleType)
diagdds$SampleType <- relevel(diagdds$SampleType, "control")

## Filter out rows where all values are zero --> from 631 to 551 keep <-
## rowSums(counts(diagdds) > 0) > 0 # diagdds <- diagdds[keep,]

# The poscounts method uses a modified geometric mean that excludes zero values
# from the calculation, which should prevent the error. counts(diagdds) <-
# as.matrix(round(counts(diagdds) + 1))
diagdds <- estimateSizeFactors(diagdds, type = "poscounts")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)
[1] "Intercept"                 "SampleType_GPA_vs_control"
res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
baseMean log2FoldChange lfcSE stat pvalue padj Domain Phylum Class Order Family Genus Species
861881 21.95936 -7.810001 2.1649848 -3.607416 0.0003093 0.0162834 k__Bacteria p__Proteobacteria c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Moraxella s__
451449 31.29902 -3.722008 1.1916255 -3.123471 0.0017873 0.0277033 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
868615 1406.00583 4.925705 1.4243705 3.458163 0.0005439 0.0162834 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Planococcaceae g__ s__
871442 427.03319 2.990001 0.8906742 3.357009 0.0007879 0.0162834 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__ 
# rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied))
write.xlsx(sigtab, file = "sigtab_GPA_vs_control.xlsx")
# subv %in% v returns a vector TRUE FALSE is.element(subv, v) returns a vector
# TRUE FALSE
library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))

# Error in checkForExperimentalReplicates(object, modelMatrix) : The design
# matrix has the same number of samples and coefficients to fit, so estimation
# of dispersion is not possible. Treating samples as replicates was deprecated
# in v1.20 and no longer supported since v1.22.

RA vs control

ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("ra002", "ra003", "ra004", "ra005",
    "ra006", "ra007", "ra008", "ra009", "ra010", "ra013", "ra014", "ra015", "ra017",
    "ra018", "ra019", "ra020", "ra022", "ra023", "ra024", "ra025", "kg001", "kg002",
    "kg003", "kg004", "kg005", "kg007", "kg009", "kg015", "kg016", "kg019", "kg020",
    "kg021", "kg022", "kg023", "kg025", "kg026", "kg027", "kg028", "kg029")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~SampleType)
diagdds$SampleType <- relevel(diagdds$SampleType, "control")
diagdds <- estimateSizeFactors(diagdds, type = "poscounts")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)
[1] "Intercept"                "SampleType_RA_vs_control"
res = results(diagdds, cooksCutoff = FALSE)
alpha = 2
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
baseMean log2FoldChange lfcSE stat pvalue padj Domain Phylum Class Order Family Genus Species
80113 4.702974 3.3268103 1.7757464 1.8734714 0.0610033 0.2707022 k__Bacteria p__Proteobacteria c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Rhizobium s__leguminosarum
521073 14.699905 2.7238523 1.6247136 1.6765123 0.0936379 0.3221164 k__Bacteria p__Proteobacteria c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Agrobacterium s__
341460 7.080004 -1.0820521 1.2030820 -0.8994001 0.3684396 0.7474060 k__Bacteria p__Proteobacteria c__Gammaproteobacteria o__Pasteurellales f__Pasteurellaceae g__Haemophilus NA
345362 9.416377 4.0951487 1.3182900 3.1064096 0.0018937 0.0268911 k__Bacteria p__Proteobacteria c__Gammaproteobacteria o__Enterobacteriales f__Enterobacteriaceae g__ s__
861881 24.012929 -23.6929822 2.2009854 -10.7647156 0.0000000 0.0000000 k__Bacteria p__Proteobacteria c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Moraxella s__
574102 4.311373 -0.6302778 1.1274145 -0.5590471 0.5761296 0.8703234 k__Bacteria p__Proteobacteria c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Enhydrobacter s__
68621 12.772930 1.1863358 0.7992093 1.4843869 0.1377063 0.4250934 k__Bacteria p__Proteobacteria c__Betaproteobacteria o__Burkholderiales f__Comamonadaceae g__Delftia s__
900973 53.442548 -0.1916407 1.1930155 -0.1606356 0.8723804 0.9582657 k__Bacteria p__Proteobacteria c__Betaproteobacteria o__Neisseriales f__Neisseriaceae g__ s__
933546 22.943268 1.1692183 2.0029770 0.5837403 0.5593951 0.8634141 k__Bacteria p__Proteobacteria c__Betaproteobacteria o__Neisseriales f__Neisseriaceae g__ s__
3393186 6.286072 1.2147295 1.9514679 0.6224696 0.5336331 0.8555449 k__Bacteria p__Proteobacteria c__Betaproteobacteria o__Neisseriales f__Neisseriaceae g__ s__
386273 22.584250 0.8961318 1.3965390 0.6416805 0.5210807 0.8555449 k__Bacteria p__Proteobacteria c__Epsilonproteobacteria o__Campylobacterales f__Campylobacteraceae g__Campylobacter s__
760967 10.949050 -0.2064463 1.9463925 -0.1060661 0.9155299 0.9582657 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Prevotellaceae g__Prevotella s__
1077373 2.811412 -1.5816399 1.9130362 -0.8267694 0.4083678 0.7491079 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Prevotellaceae g__Prevotella s__
4451252 2.657968 -0.0376489 1.6155966 -0.0233034 0.9814083 0.9814083 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Prevotellaceae g__Prevotella s__melaninogenica
4420570 15.009561 -0.1414098 1.1853273 -0.1193002 0.9050375 0.9582657 k__Bacteria p__Cyanobacteria c__Chloroplast o__Streptophyta f__ g__ s__
503315 237.998264 -1.0951498 0.8287542 -1.3214410 0.1863544 0.5292464 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Finegoldia s__
504674 66.366126 -0.3939525 0.9587696 -0.4108938 0.6811504 0.9078506 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
937813 12.784780 -0.2092928 1.0331575 -0.2025759 0.8394665 0.9582657 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
451449 40.995465 -2.5921855 1.2012800 -2.1578529 0.0309393 0.1830573 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
2751958 138.100169 -1.3061125 0.9373319 -1.3934366 0.1634878 0.4836515 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
403258 4.302868 -5.7373231 2.3118547 -2.4816971 0.0130758 0.0928384 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
1029165 18.730338 0.2182649 1.5510799 0.1407180 0.8880927 0.9582657 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Anaerococcus s__
494906 395.795222 -1.6668181 1.0023961 -1.6628339 0.0963457 0.3221164 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__[Tissierellaceae] g__Peptoniphilus s__
4475758 3.751514 1.0033482 1.4552409 0.6894723 0.4905261 0.8292227 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Veillonellaceae g__Veillonella s__dispar
1143045 8.394818 0.6740782 1.3374643 0.5039972 0.6142634 0.8871863 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Veillonellaceae g__Veillonella s__dispar
4310208 8.657849 2.5279829 1.4061298 1.7978303 0.0722039 0.3015575 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Veillonellaceae g__Veillonella s__dispar
4422456 6.654802 0.7648789 1.4258777 0.5364267 0.5916637 0.8751692 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Veillonellaceae g__Veillonella s__dispar
128382 8.078932 -0.2779723 1.8484550 -0.1503809 0.8804641 0.9582657 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Veillonellaceae g__Dialister s__
496787 50.793948 -2.6176979 1.0271446 -2.5485194 0.0108181 0.0898058 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Staphylococcaceae g__Staphylococcus s__
377874 29.041455 0.2123153 0.8390605 0.2530393 0.8002378 0.9582657 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Staphylococcaceae g__Staphylococcus s__
934488 105.114037 -1.3782006 0.7188910 -1.9171203 0.0552227 0.2706982 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Staphylococcaceae g__Staphylococcus NA
356733 1091.743683 -1.7297271 0.6263231 -2.7617167 0.0057498 0.0680397 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Staphylococcaceae g__Staphylococcus s__
868615 145.387441 7.6557230 1.5920538 4.8087087 0.0000015 0.0000539 k__Bacteria p__Firmicutes c__Bacilli o__Bacillales f__Planococcaceae g__ s__
517754 3.827454 -1.8093946 1.3934602 -1.2984903 0.1941189 0.5300940 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__
871442 90.636562 -0.3449414 0.9006785 -0.3829795 0.7017350 0.9078506 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__
4439603 2.817320 -0.4835726 1.1616468 -0.4162820 0.6772037 0.9078506 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__
579608 10.540264 -0.4849768 1.2731809 -0.3809174 0.7032645 0.9078506 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__
516966 4.791375 0.7490395 1.2291043 0.6094190 0.5422468 0.8555449 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__
882921 682.372085 1.0305857 1.3566966 0.7596287 0.4474766 0.7748984 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Aerococcaceae g__Alloiococcus s__
386088 21.595572 -2.0621066 0.8657189 -2.3819587 0.0172208 0.1111526 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Propionibacteriaceae g__Propionibacterium s__
403853 68.887221 -2.3676750 0.5369151 -4.4097754 0.0000103 0.0002449 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Propionibacteriaceae g__Propionibacterium s__acnes
4466006 4.435592 0.9148156 1.1287567 0.8104631 0.4176741 0.7491079 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Micrococcaceae g__Rothia s__dentocariosa
930926 829.404786 0.4306626 1.7974047 0.2396025 0.8106384 0.9582657 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__
446403 14.326158 0.7925423 1.6204476 0.4890885 0.6247791 0.8871863 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium NA
1047041 56.210085 -2.8614629 1.1127296 -2.5715708 0.0101238 0.0898058 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__
505626 8.730665 -2.5427326 2.0353998 -1.2492546 0.2115720 0.5452725 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__
467198 23.941453 -0.9323028 1.0018691 -0.9305635 0.3520794 0.7371856 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__
861807 355.023652 -1.5716651 0.8263841 -1.9018578 0.0571898 0.2706982 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__
377613 1556.584022 -0.8018036 0.8047079 -0.9963909 0.3190603 0.7079150 k__Bacteria p__Actinobacteria c__Actinobacteria o__Actinomycetales f__Corynebacteriaceae g__Corynebacterium s__ 
library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))

GPA vs RA

ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("micro1", "micro3", "micro4",
    "micro6", "micro7", "micro12", "micro13", "micro16", "micro17", "mw001", "mw004",
    "mw005", "mw006", "mw007", "mw009", "mw010", "mw013", "mw014", "mw015", "mw017",
    "mw021", "ra002", "ra003", "ra004", "ra005", "ra006", "ra007", "ra008", "ra009",
    "ra010", "ra013", "ra014", "ra015", "ra017", "ra018", "ra019", "ra020", "ra022",
    "ra023", "ra024", "ra025")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~SampleType)
diagdds$SampleType <- relevel(diagdds$SampleType, "RA")
diagdds <- estimateSizeFactors(diagdds, type = "poscounts")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)
[1] "Intercept"            "SampleType_GPA_vs_RA"
res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
baseMean log2FoldChange lfcSE stat pvalue padj Domain Phylum Class Order Family Genus Species
871442 424.8328 3.219397 0.8031873 4.008276 6.12e-05 0.0337623 k__Bacteria p__Firmicutes c__Bacilli o__Lactobacillales f__Streptococcaceae g__Streptococcus s__ 
library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))