Supplementary materials for Latvia WW publication

Wastewater microbiome in Latvia Edgars Liepa 14/08/2025 This is a code for analysis of Kraken2 results of wastewater microbiome data from short read metagenomic sequencing. Lets start with loading necessary libraries. ## Loading required package: permute ## ## Attaching package: ’dplyr’ ## The following objects are masked from ’package:stats’: ## ## filter, lag ## The following objects are masked from ’package:base’: ## ## intersect, setdiff, setequal, union ## Registered S3 method overwritten by ’gplots’: ## method from ## reorder.factor DescTools ## ## Attaching package: ’microbiomeMarker’ ## The following object is masked from ’package:phyloseq’: ## ## plot_heatmap Load metadata METADATA_PATH ="../sampleMetadata.csv" metadata <- read.csv(METADATA_PATH, header = TRUE, colClasses = c(Date = "character", Dairy_farming="character", Meat_production="character", Metal_processing="character", Washrooms="character")) 1 Bacterial abundance Read summarized Kraken2 report biome file that was generated with kraken-biom tool from Kraken2 otputs reading sequencing rads merged_metagenomes <- import_biom(BIOMFILE) # Clean Tax names merged_metagenomes@[email protected] <- substring(merged_metagenomes@[email protected], 4) # Set collumn names colnames(merged_metagenomes@[email protected])<- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") # Fix Sample names sample_names(merged_metagenomes) <- gsub("_k2_report$","",sample_names(merged_metagenomes)) # Add metadata sample_data(merged_metagenomes) <- sample_data(metadata) Set Sample names to descriptive City names # Create a named vector for mapping name_mapping <- setNames(metadata$Sample, rownames(metadata)) # Use the mapping to rename the columns sample_names(merged_metagenomes) <- name_mapping[sample_names(merged_metagenomes)] filtering Subset bacteria only, since we are only interested in those. merged_metagenomes <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria") Read counts in samples used in analysis plot_bar(merged_metagenomes) 2 0e+00 5e+06 1e+07 Cesis.1 Cesis.2 Cesis.3 Dobele.1 Dobele.2 Dobele.3 Jelgava.5 Jelgava.6 Jelgava.7 JurmalaS9.1 JurmalaS9.2 JurmalaS9.3 KuldigaT.1 KuldigaT.2 KuldigaT.3 Liepaja.1 Liepaja.2 Liepaja.3 Madona.1 Madona.2 Madona.3 Salaspils.1 Salaspils.2 Salaspils.3 Saldus.1 Saldus.2 Saldus.3 Sigulda.1 Sigulda.2 Sigulda.3 Smiltene.1 Smiltene.2 Smiltene.3 Talsi.1 Talsi.2 Talsi.3 Tukums.1 Tukums.2 Tukums.3 Valmiera.1 Valmiera.2 Valmiera.3 Ventspils.1 Ventspils.2 Ventspils.3 Sample Abundance Histogram of sample size distribution. We can see that one sample is cleare outlier hist(log10(sample_sums(merged_metagenomes)), breaks=50, main="Sample size distribution", xlab="Sample size (log10)",ylab="Frequency", col="#007c80") 3 Sample size distribution Sample size (log10) Frequency 6.4 6.6 6.8 7.0 02468 Show Total sequence and OTU counts. Total sequence count before filtering sum(sample_sums(merged_metagenomes)) ## [1] 376704805 Count of uniqe taxas reconstructed before filtering dim(otu_table(merged_metagenomes)) ## [1] 6837 45 Remove Outliers Talsi 3, was compromised by an unexpected industrial event - a dairy company wastewater discharge that occurred during the sampling period. This contamination is clearly visible in the data as an abnormal spike in Lactobacillus helveticus abundance. Here it can be seen that low abundance of reads in Salaspils show increased relative diversity of most abundant Species. Remove samples based on sample size distribution and outlier event on Salspils.2 and Talsi.3 4 ps_good =subset_samples(merged_metagenomes, sample_data(merged_metagenomes)$Sample != "Salaspils.2") ps_good =subset_samples(ps_good, sample_data(ps_good)$Sample != "Talsi.3") hist(log10(sample_sums(ps_good)), breaks=50, main="Sample size distribution", xlab="Sample size (log10)",ylab="Frequency", col="#007c80") Sample size distribution Sample size (log10) Frequency 6.7 6.8 6.9 7.0 0 1 2 3 4 5 boxplot(sample_sums(ps_good), main="Sequencing depth across samples remove singletons", xlab="",ylab="Number of reads", col="#a6093d") text(y=boxplot.stats(sample_sums(ps_good))$stats, labels=boxplot.stats(sample_sums(ps_good))$stats, x=1.25) 5 5.0e+06 8.0e+06 1.1e+07 Sequencing depth across samples remove singletons Number of reads 5249574 7641844.5 8611705 9329820 10807040 Difference between samples with largest and smallest sample deapth. The max difference in sequencing depth is max(sample_sums(ps_good))/min(sample_sums(ps_good)) ## [1] 2.379365 ### # Different options for filters ### # Filter Singletons ps_bact_filtered =filter_taxa(ps_good, function(x) sum(x) >1,TRUE) # Remove taxa not seen more than 3 times in at least 20% (9) of the samples. This protects against an OTU with small mean & trivially large C.V. #ps_bact_filtered = filter_taxa(ps_bact, function(x) sum(x > 3) > (0.2*length(x)), TRUE) # Filter OTUs to include only those with at least 7 reads in 5% (2) of the samples #ps_bact_filtered = filter_taxa(merged_metagenomes, function(x) sum(x > 7) > (0.05 * length(x)), TRUE) Filter low abundance tax Difference between samples with largest and smallest sample deapth after filtering singletons. 6 max(sample_sums(ps_bact_filtered))/min(sample_sums(ps_bact_filtered)) ## [1] 2.379365 Phyloseq ARG scaffolds after filtering Print the dimensions of the filtered data dim(data.frame(otu_table(ps_bact_filtered))) ## [1] 6798 43 Print the total sequence count after filtering Filtered sequence count sum(data.frame(otu_table(ps_bact_filtered))) ## [1] 361540544 Calculate and print the percentage of sequences dropped from the original dataset original_seq_count =sum(data.frame(otu_table(ps_good))) filtered_seq_count =sum(data.frame(otu_table(ps_bact_filtered))) seq_dropped_percentage =((original_seq_count -filtered_seq_count) /original_seq_count) *100 Sequence % dropped from the dataset: seq_dropped_percentage ## [1] 9.4042e-06 options(repr.plot.width=4,repr.plot.height=5) boxplot(sample_sums(ps_bact_filtered), main="Sequencing depth across samples remove singletons", xlab="", ylab="Number of reads",col="#a6093d") text(y=boxplot.stats(sample_sums(ps_bact_filtered))$stats, labels=boxplot.stats(sample_sums(ps_bact_filtered))$stats, x=1.25) 7 5.0e+06 8.0e+06 1.1e+07 Sequencing depth across samples remove singletons Number of reads 5249574 7641843 8611703 9329819 10807037 Sequence Count per Sample sample_sums(ps_bact_filtered) ## JurmalaS9.1 JurmalaS9.2 Jelgava.5 Jelgava.6 Jelgava.7 KuldigaT.1 ## 7719165 8220414 6271949 5404918 8138633 8841971 ## JurmalaS9.3 KuldigaT.2 KuldigaT.3 Talsi.1 Madona.1 Saldus.1 ## 9219143 8405969 8057134 11922250 6998348 9398233 ## Dobele.1 Tukums.1 Liepaja.2 Ventspils.1 Valmiera.1 Salaspils.1 ## 9378991 8080723 8909024 9582486 5936425 8611703 ## Dobele.2 Tukums.2 Talsi.2 Madona.2 Saldus.2 Valmiera.2 ## 10428432 9877085 10626485 10807037 7658892 9225767 ## Liepaja.3 Smiltene.2 Cesis.2 Smiltene.3 Cesis.3 Ventspils.2 ## 7624794 9020475 9530824 7416145 8258439 7932781 ## Sigulda.1 Sigulda.2 Salaspils.3 Smiltene.1 Cesis.1 Dobele.3 ## 5010685 5249574 8902442 6369985 9187228 9263098 ## Tukums.3 Saldus.3 Madona.3 Liepaja.1 Ventspils.3 Valmiera.3 ## 9295140 9862912 9109765 9364498 8432351 7169475 ## Sigulda.3 ## 6818756 Normalization Normalization by subsampling (rarefaction) 8 tab <- otu_table(ps_bact_filtered) class(tab) <- "matrix" tab <- t(tab) rare <- rarecurve(tab, step=1000,lwd=2,ylab="OTU",label=F) 0.0e+00 2.0e+06 4.0e+06 6.0e+06 8.0e+06 1.0e+07 1.2e+07 0 1000 3000 5000 Sample Size OTU # rare rarefy samples if needed. Not used currently # ps_rare = rarefy_even_depth(ps_bact_filtered, sample.size=5010685, replace=FALSE, rngseed=123, verbose=FALSE) Normalization by cumulative sum scaling (CSS) ps_CSS =microbiomeMarker::normalize(ps_bact_filtered, method="CSS") ## Default value being used. Tax Composition 9 species_abundance <- data.frame( Species = c( "Arcobacter","Aeromonas","Bacteroides", "Acinetobacter","Lactobacillus","Lactococcus", "Pseudoarcobacter","Pseudomonas" , "Cloacibacterium","Raoultella" ) ) # Calculate total counts and relative abundance species_abundance$Total_Counts <- sapply(species_abundance$Species, function(x) sum(total_counts[x])) species_abundance Genus ## Species Total_Counts ## 1 Arcobacter 78833549 ## 2 Aeromonas 19543577 ## 3 Bacteroides 22564275 ## 4 Acinetobacter 21427891 ## 5 Lactobacillus 2829250 ## 6 Lactococcus 5879146 ## 7 Pseudoarcobacter 8563781 ## 8 Pseudomonas 12293570 ## 9 Cloacibacterium 5483351 ## 10 Raoultella 1432469 Calculate mean relative abundance and variance genus_transformed =transform_sample_counts(ps_bact_filtered_genus, function(x) x /sum(x) ) abu_table <- as.data.frame(otu_table(genus_transformed)) tax_table <- as.data.frame(tax_table(genus_transformed)) species_names <- tax_table$Genus rownames(abu_table) <- tax_table$Genus variance <- apply(abu_table, 1, var) mean_percentage <- apply(abu_table, 1, mean) *100 std_dev <- sqrt(apply(abu_table, 1, var)) *100 # Select genera of interest genus_abundance <- data.frame( Genus = c("Arcobacter","Aeromonas","Bacteroides", "Acinetobacter","Citrobacter","Escherichia", "Enterobacter","Lactobacillus","Lactococcus", "Pseudoarcobacter","Pseudomonas","Cloacibacterium", "Raoultella","Klebsiella") ) 16 genus_abundance$mean_percentage <- sapply(genus_abundance$Genus, function(genus) unname(mean_percentage[genus])) genus_abundance$std_dev <- sapply(genus_abundance$Genus, function(genus) unname(std_dev[genus])) # Round relative abundance to 2 decimal places genus_abundance$mean_percentage <- round(genus_abundance$mean_percentage, 4) genus_abundance$std_dev <- round(genus_abundance$std_dev, 4) genus_abundance ## Genus mean_percentage std_dev ## 1 Arcobacter 22.2456 6.7168 ## 2 Aeromonas 5.5156 1.8185 ## 3 Bacteroides 6.4812 2.5672 ## 4 Acinetobacter 6.0865 2.9751 ## 5 Citrobacter 0.5873 0.7816 ## 6 Escherichia 0.2865 0.0938 ## 7 Enterobacter 0.2849 0.2759 ## 8 Lactobacillus 0.6896 1.6491 ## 9 Lactococcus 1.5683 2.5533 ## 10 Pseudoarcobacter 2.4541 1.3143 ## 11 Pseudomonas 3.5543 1.3135 ## 12 Cloacibacterium 1.6522 3.6093 ## 13 Raoultella 0.4085 0.8938 ## 14 Klebsiella 0.6962 0.9297 # Format output with both mean percentage and standard deviation formatted_output <- sapply(genus_abundance["Genus"], function(genus) { sprintf("%s (%.3f %% +/- %.3f %%)", genus, mean_percentage[genus], std_dev[genus]) }) # Display formatted output formatted_output ## Genus ## [1,] "Arcobacter (22.246 % +/- 6.717 %)" ## [2,] "Aeromonas (5.516 % +/- 1.818 %)" ## [3,] "Bacteroides (6.481 % +/- 2.567 %)" ## [4,] "Acinetobacter (6.086 % +/- 2.975 %)" ## [5,] "Citrobacter (0.587 % +/- 0.782 %)" ## [6,] "Escherichia (0.287 % +/- 0.094 %)" ## [7,] "Enterobacter (0.285 % +/- 0.276 %)" ## [8,] "Lactobacillus (0.690 % +/- 1.649 %)" ## [9,] "Lactococcus (1.568 % +/- 2.553 %)" ## [10,] "Pseudoarcobacter (2.454 % +/- 1.314 %)" 17 ## [11,] "Pseudomonas (3.554 % +/- 1.313 %)" ## [12,] "Cloacibacterium (1.652 % +/- 3.609 %)" ## [13,] "Raoultella (0.408 % +/- 0.894 %)" ## [14,] "Klebsiella (0.696 % +/- 0.930 %)" Now we look at bacteria that look like dominating in some cities or are outliers. sample_df <- data.frame( value = as.numeric((abu_table["Streptococcus",]))*100, city = colnames(abu_table) ) ggplot(sample_df, aes(x=as.factor(city), y=value))+geom_bar(stat = "identity")+ theme(axis.text.x = element_text(angle = 90,hjust = 1))+ geom_hline(yintercept = mean_percentage["Streptococcus"], color = "red",linewidth = 1)+ geom_text(aes(x=1,y=mean_percentage["Streptococcus"]), label = paste("Mean:",round(mean_percentage["Streptococcus"], 4)), color = "red",hjust = 0,vjust = -1.5)+ labs(y="Relative Abundance",x="City",title = "abundance ") Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926Mean: 0.6926 0 2 4 6 Cesis.1 Cesis.2 Cesis.3 Dobele.1 Dobele.2 Dobele.3 Jelgava.5 Jelgava.6 Jelgava.7 JurmalaS9.1 JurmalaS9.2 JurmalaS9.3 KuldigaT.1 KuldigaT.2 KuldigaT.3 Liepaja.1 Liepaja.2 Liepaja.3 Madona.1 Madona.2 Madona.3 Salaspils.1 Salaspils.3 Saldus.1 Saldus.2 Saldus.3 Sigulda.1 Sigulda.2 Sigulda.3 Smiltene.1 Smiltene.2 Smiltene.3 Talsi.1 Talsi.2 Tukums.1 Tukums.2 Tukums.3 Valmiera.1 Valmiera.2 Valmiera.3 Ventspils.1 Ventspils.2 Ventspils.3 City Relative Abundance abundance print(abu_table["Streptococcus",]*100) ## JurmalaS9.1 JurmalaS9.2 Jelgava.5 Jelgava.6 Jelgava.7 KuldigaT.1 ## Streptococcus 0.6940389 0.6086488 1.02555 0.8803811 1.382882 0.4111996 ## JurmalaS9.3 KuldigaT.2 KuldigaT.3 Talsi.1 Madona.1 Saldus.1 18 ## Streptococcus 0.6779063 0.4038788 0.4248247 0.4992554 0.8955967 0.3817327 ## Dobele.1 Tukums.1 Liepaja.2 Ventspils.1 Valmiera.1 Salaspils.1 ## Streptococcus 0.4812225 0.4909717 0.611809 0.4842722 0.3553665 1.155518 ## Dobele.2 Tukums.2 Talsi.2 Madona.2 Saldus.2 Valmiera.2 ## Streptococcus 0.493598 0.5508311 0.6301463 5.776454 0.4405335 0.3325698 ## Liepaja.3 Smiltene.2 Cesis.2 Smiltene.3 Cesis.3 Ventspils.2 ## Streptococcus 0.6108063 0.5805507 0.1715233 0.5953808 0.2891122 0.3269763 ## Sigulda.1 Sigulda.2 Salaspils.3 Smiltene.1 Cesis.1 Dobele.3 ## Streptococcus 0.27933 0.2979804 0.8799479 0.3089534 0.3485771 0.4380089 ## Tukums.3 Saldus.3 Madona.3 Liepaja.1 Ventspils.3 Valmiera.3 ## Streptococcus 1.156354 0.3466632 1.533399 0.5051496 0.2845456 0.4698584 ## Sigulda.3 ## Streptococcus 0.2682042 std_dev["Streptococcus"] ## Streptococcus ## 0.8503037 mean_percentage["Streptococcus"] ## Streptococcus ## 0.69257 abu_table["Streptococcus","Madona.2"]*100 ## [1] 5.776454 Valmiera had the highest levels of Cloacibacterium (Valmiera mean 13.71 % +/- 5.54 %), comprising more than a quarter of the genera. This was consistent across all sample repeats and indicated a clear difference from other cities. sample_df <- data.frame( value = as.numeric((abu_table["Cloacibacterium",]))*100, city = colnames(abu_table) ) ggplot(sample_df, aes(x=as.factor(city), y=value))+geom_bar(stat = "identity")+ geom_hline(yintercept = mean_percentage["Cloacibacterium"], color = "red",linewidth = 1)+ geom_text(aes(x=1,y=mean_percentage["Cloacibacterium"]), label = paste("Mean:",round(mean_percentage["Cloacibacterium"], 4)), color = "red",hjust = 0,vjust = -1)+ theme(axis.text.x = element_text(angle = 90,hjust = 1))+ labs(y="Relative Abundance",x="City",title = "Cloacibacterium abundance") 19 Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522Mean: 1.6522 0 5 10 15 20 Cesis.1 Cesis.2 Cesis.3 Dobele.1 Dobele.2 Dobele.3 Jelgava.5 Jelgava.6 Jelgava.7 JurmalaS9.1 JurmalaS9.2 JurmalaS9.3 KuldigaT.1 KuldigaT.2 KuldigaT.3 Liepaja.1 Liepaja.2 Liepaja.3 Madona.1 Madona.2 Madona.3 Salaspils.1 Salaspils.3 Saldus.1 Saldus.2 Saldus.3 Sigulda.1 Sigulda.2 Sigulda.3 Smiltene.1 Smiltene.2 Smiltene.3 Talsi.1 Talsi.2 Tukums.1 Tukums.2 Tukums.3 Valmiera.1 Valmiera.2 Valmiera.3 Ventspils.1 Ventspils.2 Ventspils.3 City Relative Abundance Cloacibacterium abundance print(abu_table["Cloacibacterium",]*100) ## JurmalaS9.1 JurmalaS9.2 Jelgava.5 Jelgava.6 Jelgava.7 ## Cloacibacterium 0.889483 0.7499702 0.7619668 0.5672876 0.8830507 ## KuldigaT.1 JurmalaS9.3 KuldigaT.2 KuldigaT.3 Talsi.1 Madona.1 ## Cloacibacterium 0.7173336 0.6375597 0.7414016 0.8119117 0.2701291 0.720799 ## Saldus.1 Dobele.1 Tukums.1 Liepaja.2 Ventspils.1 Valmiera.1 ## Cloacibacterium 0.3936393 0.7282066 0.2570669 0.6180769 0.7537409 8.193565 ## Salaspils.1 Dobele.2 Tukums.2 Talsi.2 Madona.2 Saldus.2 ## Cloacibacterium 1.222326 0.5104464 0.2655611 0.2336022 1.660026 0.3767129 ## Valmiera.2 Liepaja.3 Smiltene.2 Cesis.2 Smiltene.3 Cesis.3 ## Cloacibacterium 19.26936 0.5892093 2.712127 0.3117746 1.680225 0.2504079 ## Ventspils.2 Sigulda.1 Sigulda.2 Salaspils.3 Smiltene.1 Cesis.1 ## Cloacibacterium 0.6823112 0.4036743 0.4750814 0.9444778 0.7875611 0.278382 ## Dobele.3 Tukums.3 Saldus.3 Madona.3 Liepaja.1 Ventspils.3 ## Cloacibacterium 0.4015482 0.1549319 0.2689403 3.601326 0.563031 0.7197879 ## Valmiera.3 Sigulda.3 ## Cloacibacterium 13.65774 0.3303974 # Check outlier cities print( paste( "Cloacibacterium mean realative abundance:", round(mean_percentage["Cloacibacterium"],2), "% +/-", 20 round(std_dev["Cloacibacterium"],2),"%" ) ) ## [1] "Cloacibacterium mean realative abundance: 1.65 % +/- 3.61 %" print( paste( "Madona.3 relative abundance ", round(abu_table["Cloacibacterium","Madona.3"]*100,2), "%" ) ) ## [1] "Madona.3 relative abundance 3.6 %" print( paste( "Valmiera mean realative abundance:", round(apply(abu_table["Cloacibacterium",c("Valmiera.1","Valmiera.2","Valmiera.3")], 1, mean) *100,2), "% +/-", round(sqrt(apply(abu_table["Cloacibacterium",c("Valmiera.1","Valmiera.2","Valmiera.3")], 1, var)) *100,2), "%" ) ) ## [1] "Valmiera mean realative abundance: 13.71 % +/- 5.54 %" Species Show most combined read count for top 15 bacteria. abu_table <- as.data.frame(otu_table(ps_bact_filtered_species)) tax_table <- as.data.frame(tax_table(ps_bact_filtered_species)) rownames(abu_table) <- tax_table$Genus_species top_organisms <- sort(rowSums(abu_table), decreasing = TRUE)%>% head(n=15) top_organisms ## Arcobacter_cryaerophilus Pseudoarcobacter_acticola ## 49169875 8926393 ## Arcobacter_suis Aeromonas_media ## 7349119 6414859 ## Cloacibacterium_normanense Lactococcus_raffinolactis ## 5491624 5334047 ## Bacteroides_vulgatus Faecalibacterium_prausnitzii ## 4973217 4851055 ## Lactobacillus_helveticus Arcobacter_skirrowii ## 4625859 4205955 ## Acinetobacter_johnsonii Arcobacter_venerupis ## 4136710 3670322 21 ## Tolumonas_auensis Arcobacter_cibarius ## 3504160 3051260 ## Aeromonas_veronii ## 2776537 Calculate the mean relative abundance and standard deviation of the top 15 species genus_transformed =transform_sample_counts(ps_bact_filtered_species, function(x) x /sum(x) ) abu_table <- as.data.frame(otu_table(genus_transformed)) tax_table <- as.data.frame(tax_table(genus_transformed)) species_names <- tax_table$Genus_species rownames(abu_table) <- tax_table$Genus_species mean_percentage <- apply(abu_table, 1, mean) *100 std_dev <- sqrt(apply(abu_table, 1, var)) *100 top_species <- sort(mean_percentage, decreasing = TRUE)%>% head(n=15)%>% as.data.frame %>% rownames # Format output with both mean percentage and standard deviation formatted_output <- sapply(top_species, function(species) { sprintf("%s (%.2f%% +/- %.2f%%)", species, unname(mean_percentage[species]), unname(std_dev[species]) ) }) # Display formatted output formatted_output ## Arcobacter_cryaerophilus ## "Arcobacter_cryaerophilus (15.36% +/- 9.60%)" ## Pseudoarcobacter_acticola ## "Pseudoarcobacter_acticola (2.74% +/- 1.51%)" ## Arcobacter_suis ## "Arcobacter_suis (2.19% +/- 1.24%)" ## Aeromonas_media ## "Aeromonas_media (1.87% +/- 0.90%)" ## Cloacibacterium_normanense ## "Cloacibacterium_normanense (1.71% +/- 3.80%)" ## Bacteroides_vulgatus ## "Bacteroides_vulgatus (1.47% +/- 0.62%)" ## Faecalibacterium_prausnitzii ## "Faecalibacterium_prausnitzii (1.45% +/- 0.59%)" ## Lactococcus_raffinolactis ## "Lactococcus_raffinolactis (1.41% +/- 2.33%)" ## Arcobacter_skirrowii ## "Arcobacter_skirrowii (1.31% +/- 0.77%)" 22 ## Acinetobacter_johnsonii ## "Acinetobacter_johnsonii (1.21% +/- 0.53%)" ## Arcobacter_venerupis ## "Arcobacter_venerupis (1.10% +/- 0.62%)" ## Tolumonas_auensis ## "Tolumonas_auensis (1.05% +/- 0.80%)" ## Arcobacter_cibarius ## "Arcobacter_cibarius (0.95% +/- 0.45%)" ## Lactobacillus_helveticus ## "Lactobacillus_helveticus (0.94% +/- 4.87%)" ## Acidovorax_carolinensis ## "Acidovorax_carolinensis (0.83% +/- 0.60%)" Bacteria with clinical relevance Calculate the mean relative abundance and standard deviation for ESKAPE pathogens: The ESKAPE pathogens — Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa and Enterobacter spp. — are identified as critical multidrugresistant bacteria for which effective therapies are needed. species_abundance <- data.frame( Species = c( "Klebsiella_pneumoniae","Klebsiella_huaxiensis", "Acinetobacter_baumannii","Acinetobacter_johnsonii", "Pseudomonas_aeruginosa","Pseudomonas_alcaligenes", "Escherichia_coli","Citrobacter_freundii", "Citrobacter_portucalensis","Citrobacter_braakii", "Aeromonas_caviae","Aeromonas_dhakensis", "Aeromonas_veronii","Enterobacter_cloacae", "Enterococcus_faecium","Staphylococcus_aureus" ) ) # Calculate total counts and relative abundance species_abundance$Total_Counts <- sapply(species_abundance$Species, function(x) sum(total_counts[x])) species_abundance$mean_percentage <- sapply(species_abundance$Species, function(species) unname(mean_percentage[species])) species_abundance$std_dev <- sapply(species_abundance$Species, function(species) unname(std_dev[species])) # Round relative abundance to 4 decimal places species_abundance$mean_percentage <- round(species_abundance$mean_percentage, 4) species_abundance$std_dev <- round(species_abundance$std_dev, 4) species_abundance ## Species Total_Counts mean_percentage std_dev ## 1 Klebsiella_pneumoniae NA 0.1547 0.1785 23 ## 2 Klebsiella_huaxiensis NA 0.0699 0.0648 ## 3 Acinetobacter_baumannii NA 0.1501 0.0860 ## 4 Acinetobacter_johnsonii NA 1.2131 0.5254 ## 5 Pseudomonas_aeruginosa NA 0.1938 0.0728 ## 6 Pseudomonas_alcaligenes NA 0.2745 0.1566 ## 7 Escherichia_coli NA 0.2681 0.1007 ## 8 Citrobacter_freundii NA 0.1472 0.1398 ## 9 Citrobacter_portucalensis NA 0.1080 0.2804 ## 10 Citrobacter_braakii NA 0.0194 0.0104 ## 11 Aeromonas_caviae NA 0.2995 0.1583 ## 12 Aeromonas_dhakensis NA 0.0120 0.0048 ## 13 Aeromonas_veronii NA 0.8209 0.4321 ## 14 Enterobacter_cloacae NA 0.1000 0.2252 ## 15 Enterococcus_faecium NA 0.0524 0.0239 ## 16 Staphylococcus_aureus NA 0.0063 0.0017 sapply(species_abundance["Species"], function(species) { sprintf("%s (%.3f %% +/- %.3f %%)", species, unname(mean_percentage[species]), unname(std_dev[species]) ) }) ## Species ## [1,] "Klebsiella_pneumoniae (0.155 % +/- 0.178 %)" ## [2,] "Klebsiella_huaxiensis (0.070 % +/- 0.065 %)" ## [3,] "Acinetobacter_baumannii (0.150 % +/- 0.086 %)" ## [4,] "Acinetobacter_johnsonii (1.213 % +/- 0.525 %)" ## [5,] "Pseudomonas_aeruginosa (0.194 % +/- 0.073 %)" ## [6,] "Pseudomonas_alcaligenes (0.275 % +/- 0.157 %)" ## [7,] "Escherichia_coli (0.268 % +/- 0.101 %)" ## [8,] "Citrobacter_freundii (0.147 % +/- 0.140 %)" ## [9,] "Citrobacter_portucalensis (0.108 % +/- 0.280 %)" ## [10,] "Citrobacter_braakii (0.019 % +/- 0.010 %)" ## [11,] "Aeromonas_caviae (0.299 % +/- 0.158 %)" ## [12,] "Aeromonas_dhakensis (0.012 % +/- 0.005 %)" ## [13,] "Aeromonas_veronii (0.821 % +/- 0.432 %)" ## [14,] "Enterobacter_cloacae (0.100 % +/- 0.225 %)" ## [15,] "Enterococcus_faecium (0.052 % +/- 0.024 %)" ## [16,] "Staphylococcus_aureus (0.006 % +/- 0.002 %)" Arcobacter_cryaerophilus sample_df <- data.frame( value = as.numeric((abu_table["Arcobacter_cryaerophilus",]))*100, city = colnames(abu_table) ) ggplot(sample_df, aes(x=as.factor(city), y=value))+geom_bar(stat = "identity")+ geom_hline(yintercept = mean_percentage["Arcobacter_cryaerophilus"], color = "red",linewidth = 1)+ geom_text(aes(x=1,y=mean_percentage["Arcobacter_cryaerophilus"]), label = paste("Mean:",round(mean_percentage["Arcobacter_cryaerophilus"], 4)), 24 color = "red",hjust = 0,vjust = -1)+ theme(axis.text.x = element_text(angle = 90,hjust = 1))+ labs(y="Relative Abundance",x="City",title = "Arcobacter_cryaerophilus abundance") Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605Mean: 15.3605 0 20 40 60 Cesis.1 Cesis.2 Cesis.3 Dobele.1 Dobele.2 Dobele.3 Jelgava.5 Jelgava.6 Jelgava.7 JurmalaS9.1 JurmalaS9.2 JurmalaS9.3 KuldigaT.1 KuldigaT.2 KuldigaT.3 Liepaja.1 Liepaja.2 Liepaja.3 Madona.1 Madona.2 Madona.3 Salaspils.1 Salaspils.2 Salaspils.3 Saldus.1 Saldus.2 Saldus.3 Sigulda.1 Sigulda.2 Sigulda.3 Smiltene.1 Smiltene.2 Smiltene.3 Talsi.1 Talsi.2 Talsi.3 Tukums.1 Tukums.2 Tukums.3 Valmiera.1 Valmiera.2 Valmiera.3 Ventspils.1 Ventspils.2 Ventspils.3 City Relative Abundance Arcobacter_cryaerophilus abundance print(abu_table["Arcobacter_cryaerophilus",]*100) ## JurmalaS9.1 JurmalaS9.2 Jelgava.5 Jelgava.6 Jelgava.7 ## Arcobacter_cryaerophilus 12.99523 13.07077 7.254617 6.439419 9.276442 ## KuldigaT.1 JurmalaS9.3 KuldigaT.2 KuldigaT.3 Talsi.1 ## Arcobacter_cryaerophilus 18.01341 13.05277 16.0905 16.11237 16.82096 ## Madona.1 Saldus.1 Dobele.1 Tukums.1 Liepaja.2 ## Arcobacter_cryaerophilus 18.43214 24.88376 23.14072 13.85596 15.79005 ## Ventspils.1 Valmiera.1 Salaspils.1 Dobele.2 Tukums.2 ## Arcobacter_cryaerophilus 10.91912 7.674535 17.40022 21.40384 16.53625 ## Talsi.2 Madona.2 Saldus.2 Valmiera.2 Liepaja.3 ## Arcobacter_cryaerophilus 16.41362 14.30856 17.95213 8.645385 13.96083 ## Smiltene.2 Cesis.2 Smiltene.3 Cesis.3 Ventspils.2 ## Arcobacter_cryaerophilus 10.09143 8.250459 11.55965 8.904996 6.097155 ## Salaspils.2 Sigulda.1 Sigulda.2 Salaspils.3 Talsi.3 ## Arcobacter_cryaerophilus 66.83835 9.52775 10.9779 21.19461 5.465491 ## Smiltene.1 Cesis.1 Dobele.3 Tukums.3 Saldus.3 ## Arcobacter_cryaerophilus 11.04986 10.89016 29.29645 14.39981 25.8735 ## Madona.3 Liepaja.1 Ventspils.3 Valmiera.3 Sigulda.3 ## Arcobacter_cryaerophilus 17.68552 20.7957 9.392108 10.60929 11.87948 25 min(alpha_indexes_bact$InvSimpson) ## [1] 13.63599 max(alpha_indexes_bact$InvSimpson) ## [1] 169.039 min(alpha_indexes_bact$Shannon) ## [1] 5.084606 max(alpha_indexes_bact$Shannon) ## [1] 6.890168 Calculate Dunes Test to pinpoint which specific means between municipalities are most significant. library(dunn.test) dunn_results <- dunn.test(alpha_indexes_bact$InvSimpson, sample_data(ps_bact_filtered)$City, method="bh") ## Kruskal-Wallis rank sum test ## ## data: x and group ## Kruskal-Wallis chi-squared = 34.6913, df = 14, p-value = 0 ## ## ## Comparison of x by group ## (Benjamini-Hochberg) ## Col Mean-| ## Row Mean | Cesis Dobele Jelgava Jurmala Kuldiga Liep¯ aja ## ---------+------------------------------------------------------------------ ## Dobele | 3.023690 ## | 0.0262 ## | ## Jelgava | -0.650256 -3.673946 ## | 0.3515 0.0125* ## | ## Jurmala | 0.747794 -2.275896 1.398050 ## | 0.3409 0.0750 0.1850 ## | ## Kuldiga | 1.755691 -1.267999 2.405947 1.007896 ## | 0.1298 0.2150 0.0651 0.2698 ## | ## Liep¯ aja | 1.625640 -1.398050 2.275896 0.877845 -0.130051 ## | 0.1517 0.1811 0.0800 0.3023 0.4614 ## | ## Madona | 1.820717 -1.202973 2.470973 1.072922 0.065025 0.195076 32 ## | 0.1243 0.2268 0.0643 0.2564 0.4833 0.4672 ## | ## Salaspil | 2.021083 -0.683387 2.602690 1.352235 0.450745 0.567066 ## | 0.0947 0.3507 0.0486 0.1928 0.3981 0.3699 ## | ## Saldus | 2.763588 -0.260102 3.413844 2.015793 1.007896 1.137948 ## | 0.0333 0.4535 0.0168* 0.0920 0.2743 0.2392 ## | ## Sigulda | -0.130051 -3.153742 0.520204 -0.877845 -1.885742 -1.755691 ## | 0.4660 0.0282 0.3768 0.3069 0.1154 0.1385 ## | ## Smiltene | 0.552717 -2.470973 1.202973 -0.195076 -1.202973 -1.072922 ## | 0.3672 0.0590 0.2312 0.4721 0.2226 0.2609 ## | ## Talsi | 2.195565 -0.508905 2.777172 1.526717 0.625227 0.741548 ## | 0.0738 0.3773 0.0411 0.1624 0.3580 0.3389 ## | ## Tukums | 1.463076 -1.560614 2.113332 0.715281 -0.292615 -0.162564 ## | 0.1712 0.1557 0.0825 0.3459 0.4491 0.4763 ## | ## Valmiera | 1.593127 -1.430563 2.243383 0.845332 -0.162564 -0.032512 ## | 0.1535 0.1780 0.0768 0.3118 0.4713 0.4917 ## | ## Ventspil | 0.000000 -3.023690 0.650256 -0.747794 -1.755691 -1.625640 ## | 0.5000 0.0328 0.3561 0.3459 0.1340 0.1560 ## Col Mean-| ## Row Mean | Madona Salaspil Saldus Sigulda Smiltene Talsi ## ---------+------------------------------------------------------------------ ## Salaspil | 0.392584 ## | 0.4144 ## | ## Saldus | 0.942871 0.450745 ## | 0.2928 0.3936 ## | ## Sigulda | -1.950768 -2.137404 -2.893639 ## | 0.1032 0.0814 0.0333 ## | ## Smiltene | -1.267999 -1.526717 -2.210870 0.682768 ## | 0.2108 0.1585 0.0789 0.3463 ## | ## Talsi | 0.567066 0.159279 -0.276263 2.311886 1.701199 ## | 0.3745 0.4679 0.4514 0.0779 0.1373 ## | ## Tukums | -0.357640 -0.712468 -1.300512 1.593127 0.910358 -0.886950 ## | 0.4251 0.3425 0.2072 0.1496 0.3022 0.3077 ## | ## Valmiera | -0.227589 -0.596146 -1.170460 1.723178 1.040409 -0.770628 ## | 0.4629 0.3662 0.2308 0.1350 0.2653 0.3404 ## | ## Ventspil | -1.820717 -2.021083 -2.763588 0.130051 -0.552717 -2.195565 ## | 0.1287 0.0988 0.0375 0.4754 0.3716 0.0777 ## Col Mean-| ## Row Mean | Tukums Valmiera ## ---------+---------------------- 33 ## Valmiera | 0.130051 ## | 0.4707 ## | ## Ventspil | -1.463076 -1.593127 ## | 0.1751 0.1577 ## ## alpha = 0.05 ## Reject Ho if p <= alpha/2 dunn_results$P.adjusted[dunn_results$P.adjusted <0.05] ## [1] 0.02621971 0.01253870 0.04856015 0.03334894 0.01681395 0.02820851 ## [7] 0.03332044 0.04112560 0.03277463 0.03751756 dunn_results$chi2 ## [1] 34.69133 dunn_results$Z[dunn_results$P.adjusted <0.05] ## [1] 3.023691 -3.673947 2.602690 2.763588 3.413844 -3.153742 -2.893640 ## [8] 2.777172 -3.023691 -2.763588 dunn_results$comparisons[dunn_results$P.adjusted <0.05] ## [1] "Cesis - Dobele" "Dobele - Jelgava" "Jelgava - Salaspils" ## [4] "Cesis - Saldus" "Jelgava - Saldus" "Dobele - Sigulda" ## [7] "Saldus - Sigulda" "Jelgava - Talsi" "Dobele - Ventspils" ## [10] "Saldus - Ventspils" dunn_results <- dunn.test(alpha_indexes_bact$Shannon, sample_data(ps_bact_filtered)$City, method="bh") ## Kruskal-Wallis rank sum test ## ## data: x and group ## Kruskal-Wallis chi-squared = 29.1469, df = 14, p-value = 0.01 ## ## ## Comparison of x by group ## (Benjamini-Hochberg) ## Col Mean-| ## Row Mean | Cesis Dobele Jelgava Jurmala Kuldiga Liep¯ aja ## ---------+------------------------------------------------------------------ ## Dobele | 1.625640 ## | 0.1883 ## | ## Jelgava | -1.625640 -3.251280 ## | 0.1707 0.0603 34 ## | ## Jurmala | -0.260102 -1.885742 1.365537 ## | 0.4439 0.1557 0.2151 ## | ## Kuldiga | 1.007896 -0.617743 2.633537 1.267999 ## | 0.2939 0.3708 0.0634 0.2240 ## | ## Liep¯ aja | 0.032512 -1.593127 1.658153 0.292615 -0.975384 ## | 0.4917 0.1716 0.1892 0.4393 0.2981 ## | ## Madona | 0.000000 -1.625640 1.625640 0.260102 -1.007896 -0.032512 ## | 0.5000 0.1950 0.1762 0.4487 0.2993 0.4965 ## | ## Salaspil | 0.785169 -0.668847 2.239185 1.017811 -0.116321 0.756088 ## | 0.3388 0.3573 0.1015 0.3002 0.4812 0.3421 ## | ## Saldus | 1.333024 -0.292615 2.958665 1.593127 0.325128 1.300512 ## | 0.2129 0.4441 0.0324 0.1768 0.4395 0.2161 ## | ## Sigulda | -1.560614 -3.186254 0.065025 -1.300512 -2.568511 -1.593127 ## | 0.1683 0.0378 0.4880 0.2208 0.0670 0.1667 ## | ## Smiltene | -0.715281 -2.340921 0.910358 -0.455179 -1.723178 -0.747794 ## | 0.3412 0.1010 0.3121 0.4008 0.1937 0.3409 ## | ## Talsi | 1.701199 0.247182 3.155216 1.933842 0.799709 1.672119 ## | 0.1867 0.4447 0.0281 0.1468 0.3424 0.1908 ## | ## Tukums | 0.520204 -1.105435 2.145845 0.780307 -0.487692 0.487692 ## | 0.3908 0.2664 0.1116 0.3360 0.3911 0.3958 ## | ## Valmiera | -0.617743 -2.243383 1.007896 -0.357640 -1.625640 -0.650256 ## | 0.3613 0.1088 0.2888 0.4299 0.1820 0.3609 ## | ## Ventspil | -0.195076 -1.820717 1.430563 0.065025 -1.202973 -0.227589 ## | 0.4529 0.1638 0.2054 0.4929 0.2357 0.4484 ## Col Mean-| ## Row Mean | Madona Salaspil Saldus Sigulda Smiltene Talsi ## ---------+------------------------------------------------------------------ ## Salaspil | 0.785169 ## | 0.3439 ## | ## Saldus | 1.333024 0.407124 ## | 0.2178 0.4127 ## | ## Sigulda | -1.560614 -2.181025 -2.893639 ## | 0.1730 0.1094 0.0333 ## | ## Smiltene | -0.715281 -1.424936 -2.048306 0.845332 ## | 0.3459 0.2024 0.1252 0.3370 ## | ## Talsi | 1.701199 0.836217 0.508905 3.097055 2.340967 ## | 0.1945 0.3359 0.3911 0.0257 0.1122 ## | 35 ## Tukums | 0.520204 -0.319883 -0.812820 2.080819 1.235486 -1.235914 ## | 0.3957 0.4369 0.3415 0.1229 0.2275 0.2320 ## | ## Valmiera | -0.617743 -1.337695 -1.950768 0.942871 0.097538 -2.253726 ## | 0.3660 0.2210 0.1490 0.3025 0.4842 0.1156 ## | ## Ventspil | -0.195076 -0.959651 -1.528101 1.365537 0.520204 -1.875681 ## | 0.4575 0.3001 0.1748 0.2204 0.4007 0.1517 ## Col Mean-| ## Row Mean | Tukums Valmiera ## ---------+---------------------- ## Valmiera | -1.137948 ## | 0.2576 ## | ## Ventspil | -0.715281 0.422666 ## | 0.3508 0.4106 ## ## alpha = 0.05 ## Reject Ho if p <= alpha/2 dunn_results$P.adjusted[dunn_results$P.adjusted <0.05] ## [1] 0.03244233 0.03783349 0.03332044 0.02806633 0.02565321 dunn_results$chi2 ## [1] 29.14693 dunn_results$Z[dunn_results$P.adjusted <0.05] ## [1] 2.958665 -3.186255 -2.893640 3.155216 3.097056 dunn_results$comparisons[dunn_results$P.adjusted <0.05] ## [1] "Jelgava - Saldus" "Dobele - Sigulda" "Saldus - Sigulda" "Jelgava - Talsi" ## [5] "Sigulda - Talsi" Beta Diversity Now we create NMDS Bray–Curtis distance ordination plot, to discover how similar are cities between each other. theme_set(theme_bw()) #https://joey711.github.io/phyloseq/plot_ordination-examples.html GP.ord <- ordinate(ps_CSS, "NMDS","bray") 36 ## Square root transformation ## Wisconsin double standardization ## Run 0 stress 0.1070927 ## Run 1 stress 0.1353133 ## Run 2 stress 0.138403 ## Run 3 stress 0.1064777 ## ... New best solution ## ... Procrustes: rmse 0.008399121 max resid 0.04866991 ## Run 4 stress 0.1346305 ## Run 5 stress 0.1259668 ## Run 6 stress 0.1037676 ## ... New best solution ## ... Procrustes: rmse 0.05279716 max resid 0.264843 ## Run 7 stress 0.1466088 ## Run 8 stress 0.147242 ## Run 9 stress 0.09897481 ## ... New best solution ## ... Procrustes: rmse 0.01983423 max resid 0.06576992 ## Run 10 stress 0.1064749 ## Run 11 stress 0.1286362 ## Run 12 stress 0.1475917 ## Run 13 stress 0.1394729 ## Run 14 stress 0.1328386 ## Run 15 stress 0.1372535 ## Run 16 stress 0.1064749 ## Run 17 stress 0.1070927 ## Run 18 stress 0.0993338 ## ... Procrustes: rmse 0.007192529 max resid 0.02752483 ## Run 19 stress 0.1525577 ## Run 20 stress 0.1097565 ## *** Best solution was not repeated -- monoMDS stopping criteria: ## 14: stress ratio > sratmax ## 6: scale factor of the gradient < sfgrmin GP.ord ## ## Call: ## metaMDS(comm = veganifyOTU(physeq), distance = distance) ## ## global Multidimensional Scaling using monoMDS ## ## Data: wisconsin(sqrt(veganifyOTU(physeq))) ## Distance: bray ## ## Dimensions: 2 ## Stress: 0.09897481 ## Stress type 1, weak ties ## Best solution was not repeated after 20 tries ## The best solution was from try 9 (random start) ## Scaling: centring, PC rotation, halfchange scaling ## Species: expanded scores based on ’wisconsin(sqrt(veganifyOTU(physeq)))’ 37 p1 =plot_ordination(ps_CSS, GP.ord, type="Sample",color="City")+geom_point(size=6)+ geom_polygon(aes(fill=City), alpha = 1/2)+ theme( legend.position = "right", axis.text.x = element_text(angle = 90,hjust = 1,size = 14), axis.text.y = element_text(size = 14), axis.title = element_text(size = 14), legend.text = element_text(size = 19), legend.title = element_text(size = 18) ) p1 −0.10 −0.05 0.00 0.05 0.10 0.15 −0.1 0.0 0.1 0.2 NMDS1 NMDS2 City Cesis Dobele Jelgava Jurmala Kuldiga Liep.ja Madona Salaspils Saldus Sigulda Smiltene Talsi Tukums Valmiera Ventspils Now calculate Bray–Curtis distance again to use for Analysis of similarities (ANOSIM) and permutational multivariate analysis of variance using distance matrices (PERMANOVA). tax_bray <- phyloseq::distance(ps_CSS, method = "bray") anosim(tax_bray, sample_data(ps_CSS)$City, distance = "bray",permutations = 999) ## ## Call: ## anosim(x = tax_bray, grouping = sample_data(ps_CSS)$City, permutations = 999, distance = "bray") ## Dissimilarity: bray ## 38 ## ANOSIM statistic R: 0.8328 ## Significance: 0.001 ## ## Permutation: free ## Number of permutations: 999 adonis2_rez<-adonis2(tax_bray ~sample_data(ps_CSS)$City) print(adonis2_rez) ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = tax_bray ~ sample_data(ps_CSS)$City) ## Df SumOfSqs R2 F Pr(>F) ## Model 14 1.89889 0.8323 9.9258 0.001 *** ## Residual 28 0.38262 0.1677 ## Total 42 2.28151 1.0000 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Lets calculate beta dispersion to verify similarity analysis results between municipalities. beta <- betadisper(tax_bray, sample_data(ps_CSS)$City) ## Warning in betadisper(tax_bray, sample_data(ps_CSS)$City): some squared ## distances are negative and changed to zero print(anova(beta)) ## Analysis of Variance Table ## ## Response: Distances ## Df Sum Sq Mean Sq F value Pr(>F) ## Groups 14 0.025106 0.0017933 0.6778 0.7759 ## Residuals 28 0.074083 0.0026458 plot(beta) 39 Cesis Dobele Jelgava Jurmala Kuldiga Liep.ja Madona Salaspils Saldus Sigulda Smiltene Talsi Tukums Valmiera Ventspils beta method = "bray" PCoA 1 PCoA 2 −0.4 −0.2 0.0 0.2 −0.1 0.0 0.1 0.2 Hospital One of the main research questions were how does the hospital influence local waste-water microbiome and resistance gene pool? First we compare differences between sample alpha diversity measures in samples with and without regional Hospital facilities. wilcox.test(alpha_indexes_bact$Shannon~sample_data(ps_bact_filtered)$Regional_Hospital, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_bact$Shannon by sample_data(ps_bact_filtered)$Regional_Hospital ## W = 198, p-value = 0.4746 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_bact$InvSimpson~sample_data(ps_bact_filtered)$Regional_Hospital, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## 40 ## data: alpha_indexes_bact$InvSimpson by sample_data(ps_bact_filtered)$Regional_Hospital ## W = 198, p-value = 0.4746 ## alternative hypothesis: true location shift is not equal to 0 Here differences in municipalities with regional hospital between Shannon and InvSimpson indices were not found. kruskal.test(alpha_indexes_bact$Shannon ~sample_data(ps_bact_filtered)$Hospital_Type) Hospital Types ## ## Kruskal-Wallis rank sum test ## ## data: alpha_indexes_bact$Shannon by sample_data(ps_bact_filtered)$Hospital_Type ## Kruskal-Wallis chi-squared = 20.715, df = 7, p-value = 0.004215 Display Shannon and InvSimpson values in Plot plot_richness(ps_bact_filtered_species, x="Hospital_Type", measures=c("Shannon","InvSimpson"), nrow = 2,sortby = "Shannon")+ geom_boxplot() +geom_point(size=1,alpha=0.5)+labs( title = "Species Hospital type diversity") InvSimpson Shannon 0 Branch Other 2 3 Red_cross 4 Specilised 2 3 4 5 6 7 0 50 100 Hospital_Type Alpha Diversity Measure Species Hospital type diversity 41 ## ## data: alpha_indexes_ww_impact$Shannon by sample_data(ps_tax_ww_impact)$Dairy_farming ## W = 215, p-value = 0.7267 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Shannon~sample_data(ps_tax_ww_impact)$Meat_production, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Shannon by sample_data(ps_tax_ww_impact)$Meat_production ## W = 121, p-value = 0.5706 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Shannon~sample_data(ps_tax_ww_impact)$Metal_processing, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Shannon by sample_data(ps_tax_ww_impact)$Metal_processing ## W = 193, p-value = 0.1021 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Shannon~sample_data(ps_tax_ww_impact)$Washrooms, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Shannon by sample_data(ps_tax_ww_impact)$Washrooms ## W = 207, p-value = 0.005998 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Simpson~sample_data(ps_tax_ww_impact)$Dairy_farming, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Simpson by sample_data(ps_tax_ww_impact)$Dairy_farming ## W = 259, p-value = 0.4916 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Simpson~sample_data(ps_tax_ww_impact)$Meat_production, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test 48 ## ## data: alpha_indexes_ww_impact$Simpson by sample_data(ps_tax_ww_impact)$Meat_production ## W = 69, p-value = 0.02571 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Simpson~sample_data(ps_tax_ww_impact)$Metal_processing, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Simpson by sample_data(ps_tax_ww_impact)$Metal_processing ## W = 223, p-value = 0.008078 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(alpha_indexes_ww_impact$Simpson~sample_data(ps_tax_ww_impact)$Washrooms, p.adjust.method = "BH") ## ## Wilcoxon rank sum exact test ## ## data: alpha_indexes_ww_impact$Simpson by sample_data(ps_tax_ww_impact)$Washrooms ## W = 188, p-value = 0.04119 ## alternative hypothesis: true location shift is not equal to 0 Beta diversity Industrial_wastewater_impact ps_tax_ww_impact =subset_samples(ps_CSS, sample_data(ps_CSS)$Industrial_wastewater_impact != "None") tax_bray_ww <- phyloseq::distance(ps_tax_ww_impact, method = "bray") adonis2_rez<-adonis2(tax_bray_ww ~sample_data(ps_tax_ww_impact)$Industrial_wastewater_impact) print(adonis2_rez) ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = tax_bray_ww ~ sample_data(ps_tax_ww_impact)$Industrial_wastewater_impact) ## Df SumOfSqs R2 F Pr(>F) ## Model 3 0.36558 0.16024 2.4805 0.005 ** ## Residual 39 1.91593 0.83976 ## Total 42 2.28151 1.00000 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 From Food 49 ps_tax_ww_impact =subset_samples(ps_CSS, sample_data(ps_CSS)$Industrial_wastewater_impact_from_food != "None") tax_bray_ww <- phyloseq::distance(ps_tax_ww_impact, method = "bray") adonis2_rez<-adonis2(tax_bray_ww ~sample_data(ps_tax_ww_impact)$Industrial_wastewater_impact_from_food) print(adonis2_rez) ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = tax_bray_ww ~ sample_data(ps_tax_ww_impact)$Industrial_wastewater_impact_from_food) ## Df SumOfSqs R2 F Pr(>F) ## Model 3 0.3193 0.16851 2.0941 0.023 * ## Residual 31 1.5756 0.83149 ## Total 34 1.8949 1.00000 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 PERMANOVA test for industrial variables. set.seed(7) # take rest of the variables that we are looking at metadata_col =colnames(sample_data(ps_CSS))[9:13] for (param in metadata_col){ # Using Bray-Curtis distance by default if (param == "norm_factor" || param == "Sample" || param == "Date"){ next } print(param) formula_str <- paste("tax_bray ~", param) adonis2_rez<-adonis2(as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) # View results print(adonis2_rez) } ## [1] "Industrial_wastewater_impact_from_food" ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) ## Df SumOfSqs R2 F Pr(>F) ## Model 4 0.42769 0.18746 2.1917 0.003 ** ## Residual 38 1.85382 0.81254 ## Total 42 2.28151 1.00000 50 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 ## [1] "Dairy_farming" ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) ## Df SumOfSqs R2 F Pr(>F) ## Model 1 0.14634 0.06414 2.8101 0.021 * ## Residual 41 2.13516 0.93586 ## Total 42 2.28151 1.00000 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 ## [1] "Meat_production" ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) ## Df SumOfSqs R2 F Pr(>F) ## Model 1 0.07877 0.03452 1.4661 0.17 ## Residual 41 2.20274 0.96548 ## Total 42 2.28151 1.00000 ## [1] "Metal_processing" ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) ## Df SumOfSqs R2 F Pr(>F) ## Model 1 0.07088 0.03107 1.3147 0.251 ## Residual 41 2.21062 0.96893 ## Total 42 2.28151 1.00000 ## [1] "Washrooms" ## Permutation test for adonis under reduced model ## Permutation: free ## Number of permutations: 999 ## ## adonis2(formula = as.formula(formula_str), data = data.frame(sample_data(ps_CSS))) ## Df SumOfSqs R2 F Pr(>F) ## Model 1 0.16717 0.07327 3.2417 0.012 * ## Residual 41 2.11433 0.92673 ## Total 42 2.28151 1.00000 ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Core Tax Show how many OTUs are in 70% of the samples with abundance larger than and more than 0.01% relative abundance 51 dim(as.data.frame(ps_bact_filtered@tax_table [microbiome::core_members(microbiomeMarker::normalize(ps_bact_filtered, method="TSS"), detection = 1/10000 ,prevalence = 70/100)])) ## [1] 658 7 Core genes in Hospitals # Filter for hospital samples physeq_hospital <- subset_samples(ps_bact_filtered, Regional_Hospital == "yes") # Filter for non-hospital samples physeq_nonhospital <- subset_samples(ps_bact_filtered, Regional_Hospital == "no") ### We employ extra normalisation to transform absolute values to relative with TSS scalling. core_hospital <- microbiome::core(microbiomeMarker::normalize(physeq_hospital, method="TSS"), detection=1/10000,prevalence=70/100) # Core in non-hospital samples core_nonhospital <- microbiome::core(microbiomeMarker::normalize(physeq_nonhospital, method="TSS"), detection=1/10000,prevalence=70/100) # Extract core feature names hospital_core_genes <- taxa_names(core_hospital) nonhospital_core_genes <- taxa_names(core_nonhospital) # Compare (intersect, union, setdiff) shared_cores <- intersect(hospital_core_genes, nonhospital_core_genes) unique_hospital <- setdiff(hospital_core_genes, nonhospital_core_genes) unique_nonhospital <- setdiff(nonhospital_core_genes, hospital_core_genes) Print core taxa unique to hospital samples physeq_nonhospital@tax_table[unique_hospital][,c("Genus","Species")] ## Taxonomy Table: [73 taxa by 2 taxonomic ranks]: ## Genus Species ## 631220 "Desulfovibrio" "sp. G11" ## 44742 "Desulfovibrio" "fairfieldensis" ## 345370 "Desulfovibrio" "carbinoliphilus" ## 2293 "Desulfobacter" "postgatei" ## 253237 "Pseudomonas" "sp. phDV1" ## 1028989 "Pseudomonas" "sp. StFLB209" 52 ## 2320270 "Pseudomonas" "sp. DG56-2" ## 1981174 "Pseudomonas" "sp. M30-35" ## 47878 "Pseudomonas" "azotoformans" ## 237610 "Pseudomonas" "psychrotolerans" ## 101564 "Pseudomonas" "alcaliphila" ## 122355 "Pseudomonas" "psychrophila" ## 353 "Azotobacter" "chroococcum" ## 2560029 "Aeromonas" "sp. 2692-1" ## 196024 "Aeromonas" "dhakensis" ## 347534 "Zobellella" "denitrificans" ## 83771 "Succinivibrio" "dextrinosolvens" ## 57706 "Citrobacter" "braakii" ## 561 "Escherichia" "" ## 61646 "Lelliottia" "amnigena" ## 204037 "Dickeya" "" ## 453783 "Lysobacter" "soli" ## 1605891 "Lysobacter" "maris" ## 56458 "Xanthomonas" "sacchari" ## 1176533 "Luteimonas" "sp. Gr-4" ## 2508168 "Luteimonas" "sp. YGD11-2" ## 666 "Vibrio" "cholerae" ## 729 "Haemophilus" "parainfluenzae" ## 87883 "Burkholderia" "multivorans" ## 105560 "Methylibium" "petroleiphilum" ## 2494234 "Sutterella" "megalosphaeroides" ## 1499392 "" "" ## 1867715 "Bosea" "sp. Tri-49" ## 82115 "" "" ## 1842534 "Sinorhizobium" "sp. RAC02" ## 1545044 "Paracoccus" "sanguinis" ## 135740 "Paracoccus" "kondratievae" ## 2065379 "Paracoccus" "jeotgali" ## 34003 "Paracoccus" "aminophilus" ## 1850250 "Rhodobacter" "sp. LPB0142" ## 33050 "Sphingopyxis" "macrogoltabida" ## 165695 "Sphingobium" "" ## 293 "Brevundimonas" "diminuta" ## 930 "Acidithiobacillus" "thiooxidans" ## 1898203 "" "Lachnospiraceae bacterium" ## 1531 "Enterocloster" "clostridioformis" ## 29370 "Lacrimispora" "sphenoides" ## 2696063 "Anaerocolumna" "sp. CBA3638" ## 1161942 "Ruminococcus" "champanellensis" ## 1731 "Peptoclostridium" "acidaminophilum" ## 1042156 "Clostridium" "sp. SY8519" ## 51515 "Dehalobacterium" "formicoaceticum" ## 1297617 "Intestinimonas" "butyriciproducens" ## 33952 "Acetobacterium" "woodii" ## 1624 "Lactobacillus" "salivarius" ## 1354 "Enterococcus" "hirae" ## 33945 "Enterococcus" "avium" ## 29466 "Veillonella" "parvula" ## 33025 "Phascolarctobacterium" "faecium" ## 626940 "Phascolarctobacterium" "succinatutens" 53 ## 187327 "Acidaminococcus" "intestini" ## 1716 "Corynebacterium" "" ## 1725 "Corynebacterium" "xerosis" ## 53461 "Nakamurella" "multipartita" ## 1433126 "Mucinivorans" "hirudinis" ## 2133944 "Alloprevotella" "sp. E39" ## 2518971 "Duncaniella" "dubosii" ## 2530390 "Muribaculum" "sp. TLL-A4" ## 1796646 "Muribaculum" "intestinale" ## 84566 "" "" ## 330214 "Nitrospira" "defluvii" ## 1796921 "Ereboglobus" "luteus" ## 856 "Fusobacterium" "varium" 54