Modelling soil bacterial traits across environments with the trait sequence database ampliconTraits and the R package MicEnvMod - Applying the workflow to real world data
classication of traits
as described in https://github.com/jdonhauser/ampliconTraits using usearch v11.0.667 e.g.
usearch -sinaps ASVs.fasta -db genome_size10.fasta -attr genome_size10 -tabbedout genomeclassification.txt -strand pluscalculation of sequence identity with the top hit
e.g for genome size
usearch -usearch_global dna-sequences.fasta -db genome_size10.fasta -strand plus -id 0.5 \
-maxaccepts 8 -maxrejects 128 -top_hit_only \
-userout genome_size10_topHit.txt -userfields query+target+id libraries
library(vegan)
library(grid)
library(reshape2)
library(egg)
library(RColorBrewer)
library(viridis)
library(ggforce)
library(sp)
library(missMDA)
library(cluster)
library(usdm)
library(MASS)
library(corrplot)
library(ecospat)
library(randomForest)
library(fitdistrplus)
library(Hmisc)
library(rgdal)
library(raster)
library(dismo)
library(maptools)
library(terra)
library(viridis)
library(MicEnvMod)
library(stars)
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] stars_0.6-4 sf_1.0-12 abind_1.4-5
## [4] MicEnvMod_0.0.0.9000 maptools_1.1-4 dismo_1.3-9
## [7] raster_3.5-29 rgdal_1.6-4 Hmisc_4.7-1
## [10] Formula_1.2-4 fitdistrplus_1.1-8 survival_3.4-0
## [13] randomForest_4.7-1.1 ecospat_3.5.1 corrplot_0.92
## [16] MASS_7.3-57 usdm_2.1-6 terra_1.7-23
## [19] cluster_2.1.4 missMDA_1.18 sp_2.0-0
## [22] ggforce_0.4.1 viridis_0.6.4 viridisLite_0.4.2
## [25] RColorBrewer_1.1-3 egg_0.4.5 ggplot2_3.4.4
## [28] gridExtra_2.3 reshape2_1.4.4 vegan_2.6-2
## [31] lattice_0.20-45 permute_0.9-7
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 plyr_1.8.7 splines_4.1.3
## [4] digest_0.6.29 foreach_1.5.2 htmltools_0.5.5
## [7] earth_5.3.2 fansi_1.0.3 magrittr_2.0.3
## [10] checkmate_2.1.0 doParallel_1.0.17 ks_1.14.0
## [13] jpeg_0.1-9 colorspace_2.1-0 ggrepel_0.9.1
## [16] xfun_0.32 dplyr_1.0.9 jsonlite_1.8.0
## [19] iterators_1.0.14 ape_5.6-2 glue_1.6.2
## [22] polyclip_1.10-4 PresenceAbsence_1.1.11 gtable_0.3.4
## [25] emmeans_1.8.7 scales_1.2.1 mvtnorm_1.1-3
## [28] DBI_1.1.3 Rcpp_1.0.10 plotrix_3.8-2
## [31] xtable_1.8-4 htmlTable_2.4.1 units_0.8-0
## [34] flashClust_1.01-2 foreign_0.8-82 proxy_0.4-27
## [37] mclust_6.0.0 DT_0.28 htmlwidgets_1.6.2
## [40] nabor_0.5.0 mice_3.15.0 pkgconfig_2.0.3
## [43] reshape_0.8.9 farver_2.1.1 nnet_7.3-17
## [46] multcompView_0.1-9 sass_0.4.2 deldir_1.0-6
## [49] utf8_1.2.2 tidyselect_1.1.2 rlang_1.1.1
## [52] munsell_0.5.0 TeachingDemos_2.12 tools_4.1.3
## [55] cachem_1.0.6 xgboost_1.6.0.1 cli_3.3.0
## [58] generics_0.1.3 ade4_1.7-19 broom_1.0.1
## [61] evaluate_0.16 stringr_1.4.1 fastmap_1.1.0
## [64] yaml_2.3.5 maxnet_0.1.4 knitr_1.40
## [67] purrr_0.3.4 nlme_3.1-159 pracma_2.4.2
## [70] leaps_3.1 biomod2_4.2-4 compiler_4.1.3
## [73] rstudioapi_0.14 png_0.1-7 e1071_1.7-11
## [76] tibble_3.2.1 tweenr_2.0.2 bslib_0.4.0
## [79] stringi_1.7.6 plotmo_3.6.2 poibin_1.5
## [82] Matrix_1.4-1 classInt_0.4-7 gbm_2.1.8.1
## [85] vctrs_0.6.1 pillar_1.9.0 lifecycle_1.0.4
## [88] jquerylib_0.1.4 estimability_1.4.1 data.table_1.14.2
## [91] R6_2.5.1 latticeExtra_0.6-30 KernSmooth_2.23-20
## [94] codetools_0.2-18 gtools_3.9.3 assertthat_0.2.1
## [97] withr_2.5.2 mgcv_1.8-40 parallel_4.1.3
## [100] rpart_4.1.16 tidyr_1.2.0 coda_0.19-4
## [103] class_7.3-20 rmarkdown_2.16 mda_0.5-3
## [106] pROC_1.18.0 scatterplot3d_0.3-44 base64enc_0.1-3
## [109] FactoMineR_2.8 interp_1.1-3setup
plot.theme1 <- theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "white",
colour = "black",
size = 0.5, linetype = "solid"),
panel.border= element_rect(fill=NA,size = 0.5, linetype = 'solid',colour = "black"),
axis.text.x = element_text(size=13),axis.text.y = element_text(size=13),legend.text = element_text(size=13),
axis.title = element_text(size=14),
legend.title = element_text(color = "black", size = 14),
strip.text.x = element_text(size=14),
strip.background = element_rect(colour="black", fill="white")
)
# colors for landcover
covLeg <- read.csv("./igbp_legend.csv")
# create color vector with the colors indicated in the legend table
colLC <- c()
for (i in 1:nrow(covLeg)){
colLC[i] <- rgb(covLeg[i,'r']/360,covLeg[i,'g']/360,covLeg[i,'b']/360)
}
names(colLC) <- covLeg$Name
### Function for aggregating asvs *********************************************************************** ####
#countab: counttable,taxa are rows
#taxo=taxonomy table
#col2matchcount vector of rownames or column in counttab to which order of taxonomy table should be matched
#col2matchtax:vector of rownames or column in taxo to match
#Taxlevel: taxonomic level on which should be aggregated
#Samp: sample data
#fac: factor in Sample of which the mean should be calculated
#Summarize: should rare taxa be summarized and represented as others
#sumlevel: abundance threshold below which taxa are summarkued as others
# long: transform to long format, default TRUE
AbuBarTable <- function(countab,taxo,col2matchcount,col2matchtax,Taxlevel,Samp,fac,Summarize=F,sumlevel,long=T){
#replace NA with unclassified
taxo <- as.data.frame(apply(taxo,2,function(x){
sapply(x,function(y){ifelse(is.na(y),"unclassified",y)})
}))
tax <- taxo[match(col2matchcount,col2matchtax),]
taxabu <- aggregate(countab,list(tax[,Taxlevel]),sum)
rownames(taxabu)<- taxabu[,1]
taxabu <- taxabu[,-1]
taxabu.rel <- apply(taxabu,2,function(x)x/sum(x))
taxabu.rel <- as.data.frame(t(taxabu.rel))
taxabu.rel <- taxabu.rel[rownames(Samp),]
for(i in 1:ncol(Samp)){
Samp[,i] <- as.character(Samp[,i])
}
taxabu.rel.mean <- aggregate(taxabu.rel,list(Samp[,fac]),mean)
colnames(taxabu.rel.mean)[1] <- fac
rownames(taxabu.rel.mean) <- taxabu.rel.mean[,1]
a=c()
for (i in colnames(taxabu.rel.mean)){
a[i] <- is.numeric(taxabu.rel.mean[,i])
}
taxabu.rel.mean <- taxabu.rel.mean[,a]
if(Summarize==T){
num=taxabu.rel.mean
num.l=as.list(as.data.frame(t(num)))
num.l=lapply(num.l,function(x){
names(x)=colnames(num)
Others=sum(x[which(x<sumlevel)])
x=x[-which(x<sumlevel)]
names(Others)="Others"
x=c(x,Others)
})
taxabu.rel.mean <- as.data.frame(do.call(rbind, lapply(num.l, "[",unique(as.vector(unlist(sapply(num.l,names)))))))
colnames(taxabu.rel.mean) <- unique(as.vector(unlist(sapply(num.l,names))))
taxabu.rel.mean <- apply(taxabu.rel.mean,2,function(x){
sapply(x,function(y){ifelse(is.na(y),0,y)})
})
}
a <- as.data.frame(Samp[!duplicated(Samp[,fac]),])
if(all(rownames(taxabu.rel.mean)==a[,fac])){
taxabu.rel.mean <- cbind(taxabu.rel.mean,a)
}else{
taxabu.rel.mean <- taxabu.rel.mean[match(a[,fac],rownames(taxabu.rel.mean)),]
taxabu.rel.mean <- cbind(taxabu.rel.mean,a)
}
require(reshape2)
if (long){
taxabu.rel.mean.long <- melt(taxabu.rel.mean)
a=unique(as.character(taxabu.rel.mean.long$variable))
if("Others"%in%taxabu.rel.mean.long$variable){
taxabu.rel.mean.long$variable <- factor(taxabu.rel.mean.long$variable,levels=c(a[-(which(a=="Others"))],"Others"))
}
taxabu.rel.mean.long[,fac] <- factor(taxabu.rel.mean.long[,fac],levels=unique(taxabu.rel.mean.long[,fac]))
return(taxabu.rel.mean.long)
}else{
return(taxabu.rel.mean)
}
}Import files
load("./traitAnnotationProperties.RData" )Properties of trait annotations across the dataset
Calculate community weighted trait means for continuous traits for bootstrap > 70 and sequence identity with tophit > 80
# based on mean of interval
# d1_lo5 has only one value for this dataset, omit
# continuous traits
cont <- c("d1_lo10","d1_lo20","d1_lo30","d1_up10","d1_up20","d1_up30","d1_up5","d2_lo10","d2_lo20","d2_lo30",
"d2_lo5","d2_up10","d2_up20","d2_up30","d2_up5","doubling_h10","doubling_h20","doubling_h30","doubling_h40","doubling_h50",
"doubling_h5","genome_size10","genome_size20","genome_size5","optimum_ph10","optimum_ph20",
"optimum_ph5","optimum_tmp10","optimum_tmp20","rRNA16S_genes_exact_int","rRNA16S_genes10",
"rRNA16S_genes5")
traits_WA <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 70 or identity is < 80 or NA
a$Value <- ifelse(a$Bootstrap < 70 | a$topHit_ID < 80 | is.na(a$topHit_ID), 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
traits_WA2 <- traits_WA
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df <- do.call(cbind, traits_WA2)
colnames(traits_WA_df) <- names(traits_WA2)Correlation matrices for different intervals for the same trait
# d1lo
ecospat.cor.plot(traits_WA_df[,grep('d1_lo',colnames(traits_WA_df))])
# d1_up
ecospat.cor.plot(traits_WA_df[,grep('d1_up',colnames(traits_WA_df))][,c(4,1,2,3)])
# d2_lo
ecospat.cor.plot(traits_WA_df[,grep('d2_lo',colnames(traits_WA_df))][,c(4,1,2,3)])
# d2_up
ecospat.cor.plot(traits_WA_df[,grep('d2_up',colnames(traits_WA_df))][,c(4,1,2,3)])
#doubling_h
ecospat.cor.plot(traits_WA_df[,grep('doubling_h',colnames(traits_WA_df))][,c(6,1,2,3,4,5)])
#genome_size
ecospat.cor.plot(traits_WA_df[,grep('genome_size',colnames(traits_WA_df))][,c(3,1,2)])
#optimum_ph
ecospat.cor.plot(traits_WA_df[,grep('optimum_ph',colnames(traits_WA_df))][,c(3,1,2)])
#optimum_tmp
ecospat.cor.plot(traits_WA_df[,grep('optimum_tmp',colnames(traits_WA_df))])
#rRNA16S_genes
ecospat.cor.plot(traits_WA_df[,grep('rRNA16S_genes',colnames(traits_WA_df))][,c(3,1,2)])
Correlation matrices of abundance weighted average for different bootstrap cutoffs
with and without cutoff for sequence identity ### bootstrap > 70 considered classified, no seqID cutoff
traits_WA2 <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 70
a$Value <- ifelse(a$Bootstrap < 70, 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA2[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df_70 <- do.call(cbind, traits_WA2)
colnames(traits_WA_df_70) <- names(traits_WA2)bootstrap > 60 considered classified, seqID > 80
traits_WA2 <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 60 or identity is < 80 or NA
a$Value <- ifelse(a$Bootstrap < 60 | a$topHit_ID < 80 | is.na(a$topHit_ID), 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA2[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df_60_80 <- do.call(cbind, traits_WA2)
colnames(traits_WA_df_60_80) <- names(traits_WA2)bootstrap > 60 considered classified, no seqID cutoff
traits_WA2 <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 60
a$Value <- ifelse(a$Bootstrap < 60, 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA2[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df_60 <- do.call(cbind, traits_WA2)
colnames(traits_WA_df_60) <- names(traits_WA2)bootstrap > 50 considered classified, seqID > 80
traits_WA2 <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 50 or identity is < 80 or NA
a$Value <- ifelse(a$Bootstrap < 50 | a$topHit_ID < 80 | is.na(a$topHit_ID), 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA2[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df_50_80 <- do.call(cbind, traits_WA2)
colnames(traits_WA_df_50_80) <- names(traits_WA2)bootstrap > 50 considered classified, no seqID cutoff
traits_WA2 <- list()
for (i in cont){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 50
a$Value <- ifelse(a$Bootstrap < 50, 'unclassified', a$Value)
# extract lower and upper boundary of interval and calculate mean
for (j in unique(a$Value)){
if (j=='unclassified'){
a[a$Value==j,'lo'] <- NA
a[a$Value==j,'up'] <- NA
}else{
a[a$Value==j,'lo'] <- gsub('-.*','', j)
a[a$Value==j,'up'] <- gsub('.*-','', j)
}
}
a$lo <- as.numeric(a$lo)
a$up <- as.numeric(a$up)
# calculate mean
a$mean <- (a$lo+a$up)/2
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "mean", sam, "Sample.ID")
# calculated weighted average
for (j in unique(abu$Sample.ID)){
b <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','variable']))
# abundance relative to sum of classified
c <- as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])) /
sum(as.numeric(as.character(abu[abu$Sample.ID==j & abu$variable!='unclassified','value'])))
abu[abu$Sample.ID==j, paste0(i,'_weightedMean')] <- sum(b*c)
rm(b,c)
}
# dereplicate, so each sample occurs only once
abu <- abu[!duplicated(abu$Sample.ID),]
traits_WA2[[i]] <- abu
}## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables
## Using Sample.ID, Gradient, Site, GxS as id variables# keep only last column( weighted average)
for (i in names(traits_WA2)){
traits_WA2[[i]] <- as.data.frame(traits_WA2[[i]][,ncol(traits_WA2[[i]])])
}
traits_WA_df_50 <- do.call(cbind, traits_WA2)
colnames(traits_WA_df_50) <- names(traits_WA2)Correlation matrices for different bootstrap cutoffs
### add cutoff to column names and merge tables
colnames(traits_WA_df) <- paste0(colnames(traits_WA_df),'_7080')
colnames(traits_WA_df_70) <- paste0(colnames(traits_WA_df_70),'_70')
colnames(traits_WA_df_60_80) <- paste0(colnames(traits_WA_df_60_80),'_6080')
colnames(traits_WA_df_60) <- paste0(colnames(traits_WA_df_60),'_60')
colnames(traits_WA_df_50_80) <- paste0(colnames(traits_WA_df_50_80),'_5080')
colnames(traits_WA_df_50) <- paste0(colnames(traits_WA_df_50),'_50')
# merge tables
traits_WA_df_all <- cbind(traits_WA_df,traits_WA_df_70,traits_WA_df_60_80,traits_WA_df_60,traits_WA_df_50_80,
traits_WA_df_50)
# for one of the intervals where stable results were obtained
# d1lo20
ecospat.cor.plot(traits_WA_df_all[,grep('d1_lo20',colnames(traits_WA_df_all))])
# d1lo20 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('d1_lo20.*80$',colnames(traits_WA_df_all))])
# d1up30
ecospat.cor.plot(traits_WA_df_all[,grep('d1_up30',colnames(traits_WA_df_all))])
# d1up30 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('d1_up30.*80$',colnames(traits_WA_df_all))])
# d2lo30
ecospat.cor.plot(traits_WA_df_all[,grep('d2_lo30',colnames(traits_WA_df_all))])
# d2lo30 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('d2_lo30.*80$',colnames(traits_WA_df_all))])
# d2up30
ecospat.cor.plot(traits_WA_df_all[,grep('d2_up30',colnames(traits_WA_df_all))])
# d2up30 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('d1_up30.*80$',colnames(traits_WA_df_all))])
# doubling_h30
ecospat.cor.plot(traits_WA_df_all[,grep('doubling_h30',colnames(traits_WA_df_all))])
# doubling_h30 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('doubling_h30.*80$',colnames(traits_WA_df_all))])
# genome_size10
ecospat.cor.plot(traits_WA_df_all[,grep('genome_size10',colnames(traits_WA_df_all))])
# genome_size10 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('genome_size10.*80$',colnames(traits_WA_df_all))])
# optimum_ph20
ecospat.cor.plot(traits_WA_df_all[,grep('optimum_ph20',colnames(traits_WA_df_all))])
# optimum_ph20 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('optimum_ph20.*80$',colnames(traits_WA_df_all))])
# optimum_tmp20
ecospat.cor.plot(traits_WA_df_all[,grep('optimum_tmp20',colnames(traits_WA_df_all))])
# optimum_tmp20 - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('optimum_tmp20.*80$',colnames(traits_WA_df_all))])
# rRNA16S_genes exact
ecospat.cor.plot(traits_WA_df_all[,grep('rRNA16S_genes_exact_int',colnames(traits_WA_df_all))])
# rRNA16S_genes exact - only with sequence identity cutoff
ecospat.cor.plot(traits_WA_df_all[,grep('rRNA16S_genes_exact_int.*80$',colnames(traits_WA_df_all))])
Fraction of unclassified samples
# aggregate abundances of trait levels
sam <- sam[,c(1:2,4)] # subset to non-numeric
traits_abu <- list()
for (i in names(traits)){
a <- traits[[i]]
rownames(a) <- a[,1]
# replace values with unclassified if bootstrap is < 70 or identity is < 80 or NA
a$Value <- ifelse(a$Bootstrap < 70 | a$topHit_ID < 80 | is.na(a$topHit_ID), 'unclassified', a$Value)
abu <- AbuBarTable(asv, a, rownames(asv), rownames(a), "Value", sam, "Sample.ID")
traits_abu[[i]] <- abu
}## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables
## Using Sample.ID, Gradient, GxS as id variables# add name of trait
for (i in names(traits_abu)){
traits_abu[[i]]$trait <- i
}
traits_abu_df <- do.call(rbind, traits_abu)
# subset to unclassified
traits_abu_df <- traits_abu_df[traits_abu_df$variable=='unclassified',]
# subset with continuous traits
traits_abu_df_cont <- traits_abu_df[grep('[0-9]$|_int',traits_abu_df$trait),]
# add column with interval and trait without number
traits_abu_df_cont$trait2 <- gsub('[0-9]*$','',traits_abu_df_cont$trait)
traits_abu_df_cont$intervals <- gsub('.*_','',traits_abu_df_cont$trait)
traits_abu_df_cont$intervals <- gsub('[a-z]*','',traits_abu_df_cont$intervals)
traits_abu_df_cont[traits_abu_df_cont$trait=="rRNA16S_genes_exact_int",'intervals'] <- 'exact'
traits_abu_df_cont[traits_abu_df_cont$trait2=="rRNA16S_genes_exact_int",'trait2'] <- 'rRNA16S_genes'
# order intervals
traits_abu_df_cont$intervals <- factor(traits_abu_df_cont$intervals, levels = c("5", "10", "20", "30", "40", "50", "exact"))
# boxplot for continuous variables as a function of interval
ggplot(traits_abu_df_cont, aes(x = intervals, y = value)) + geom_boxplot(outlier.size = 0.3) +
facet_grid(~trait2, scales = 'free', space = 'free') + plot.theme1 +
geom_jitter(size = 0.3, alpha = 0.3, color = 'blue') +
labs(y = '% unclassified', x = 'number of intervals') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
# subset with categorigal traits
traits_abu_df_cat <- traits_abu_df[-(grep('[0-9]$|exact',traits_abu_df$trait)),]
# boxplot for categorical variables
ggplot(traits_abu_df_cat, aes(x = trait, y = value)) + geom_boxplot(outlier.size = 0.3) +
plot.theme1 +
geom_jitter(size = 0.3, alpha = 0.3, color = 'blue') +
labs(y = '% unclassified', x = 'Trait') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
bootstrap values as a function of number of intervals
traits_df <- do.call(rbind,traits)
# subset with continuous traits
traits_df_cont <- traits_df[grep('[0-9]$|_int',traits_df$Trait),]
# add column with interval and trait without number
traits_df_cont$trait2 <- gsub('[0-9]*$','',traits_df_cont$Trait)
traits_df_cont$intervals <- gsub('.*_','',traits_df_cont$Trait)
traits_df_cont$intervals <- gsub('[a-z]*','',traits_df_cont$intervals)
traits_df_cont[traits_df_cont$Trait=="rRNA16S_genes_exact_int",'intervals'] <- 'exact'
traits_df_cont[traits_df_cont$trait2=="rRNA16S_genes_exact_int",'trait2'] <- 'rRNA16S_genes'
# order intervals
traits_df_cont$intervals <- factor(traits_df_cont$intervals, levels = c("5", "10", "20", "30", "40", "50", "exact"))
# boxplot for continuous variables as a function of interval
ggplot(traits_df_cont, aes(x = intervals, y = Bootstrap)) + geom_boxplot(outlier.size = 0.3) +
facet_grid(~trait2, scales = 'free', space = 'free') + plot.theme1 +
# geom_jitter(size = 0.3, alpha = 0.1, color = 'blue') +
labs(y = 'Bootstrap value', x = 'number of intervals') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
# for categorical traits
traits_df_cat <- traits_df[-(grep('[0-9]$|exact',traits_df$Trait)),]
ggplot(traits_df_cat, aes(x = Trait, y = Bootstrap)) + geom_boxplot(outlier.size = 0.3) +
plot.theme1 +
# geom_jitter(size = 0.3, alpha = 0.3, color = 'blue') +
labs(y = '% unclassified', x = 'Trait') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
sequence identity with tophit properties
distribution across dataset
#retain only one interval for each continuous trait (seq ID with tophit should be the same for different intervals)
traits_df2 <- traits_df[traits_df$Trait%in%c("cell_shape","d1_lo10", "d1_up10","d2_lo10","d2_up10", "doubling_h10", "genome_size10",
"gram_stain", "metabolism", "motility", "optimum_ph10","optimum_tmp10",
"range_salinity", "range_tmp", "rRNA16S_genes_exact", "sporulation"),]
ggplot(traits_df2, aes(x = Trait, y = topHit_ID)) + geom_boxplot(outlier.size = 0.3) +
plot.theme1 +
# geom_jitter(size = 0.3, alpha = 0.1, color = 'blue') +
labs(y = '% sequence identity top hit', x = 'number of intervals') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
# for each sample and trait: average of ASVs that occur in this sample (not weighted by abundance)
meanseqID <- data.frame()
for (i in unique(traits_df2$Trait)){
a <- traits_df2[traits_df2$Trait==i,]
for(j in colnames(asv)){
asvList <- rownames(asv[asv[,j]!=0,])
meanseqID[j,i] <- mean(a[a$ASV%in%asvList,'topHit_ID'],na.rm = T)
}
}
rm(asvList,a)
meanseqID <- melt(meanseqID)## No id variables; using all as measure variablesggplot(meanseqID, aes(x=variable,y=value)) + geom_boxplot(outlier.size = 0.3) +
plot.theme1 +
geom_jitter(size = 0.3, alpha = 0.1, color = 'blue') +
labs(y = 'mean % sequence identity \n top hit', x = 'number of intervals') +
theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1))
### bootstrap values as function of sequence identity for each trait
ggplot(traits_df, aes(x = topHit_ID, y = Bootstrap)) + geom_point(alpha = 0.2, size = 0.8) +
geom_smooth(formula = y ~ s(x)) +
plot.theme1 +
facet_wrap(~Trait)## `geom_smooth()` using method = 'gam'