#library(dada2)
library(phyloseq)
library(permute)
library(lattice)
library(vegan)
library(ggplot2)
library(tidyr)
library(dplyr)
##
## 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
library(pheatmap)
library(GUniFrac)
## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httr
#library(metagMisc)
#library(raster)
library(pals)
library(RColorBrewer)
library(ragg)
library(ggpubr)
library("phyloseq")
library(readr)
#Set theme
theme_set(theme_bw() + theme(
plot.title = element_text(size=20, color="black"),
axis.text.x = element_text(size=15, color="black"),
axis.text.y = element_text(size=15, color="black"),
axis.title.x = element_text(size=15),
axis.title.y = element_text(size=15),
legend.text = element_text(size=12),
legend.title = element_text(size=15),
# legend.position = "bottom",
# legend.key=element_blank(),
# legend.key.size = unit(0.5, "cm"),
# legend.spacing.x = unit(0.1, "cm"),
# legend.spacing.y = unit(0.1, "cm"),
panel.background = element_blank(),
#panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.background = element_blank()))
#Create phyloseq object
#read in files
asvtab <- read.table("16S_dada2_ASVs_04.23.2025.count_table", header=T, row.names=1,
check.names=F, sep="\t")
metadata_full <- read.csv("fin_read_track_16S_04162025.csv", row.names=1)
metadata_full$treatment <- as.factor(metadata_full$treatment)
taxtab <- as.matrix(read.csv("16S_dada2_ASV_taxonomy_04.23.2025.csv", header = T, row.names=1))
asvtab <- otu_table(asvtab, taxa_are_rows = FALSE)
taxtab <- tax_table(taxtab)
metadata <- sample_data(metadata_full)
#combine into phyloseq object
pseqtest <- phyloseq(asvtab, taxtab, metadata)
#extract easy-to-use sample data
samdftest <- data.frame(sample_data(pseqtest))
cured_asvs<- pseqtest@otu_table
summary(rowSums(cured_asvs))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13139 27540 33702 32554 37856 49380
sum(rowSums(cured_asvs))
## [1] 1302144
ncol(cured_asvs)
## [1] 12461
##remove singletons, doubletons, and ASVs with less than 5 reads across all samples.
cured_asvs <- cured_asvs[,colSums(cured_asvs)>5]
summary(rowSums(cured_asvs))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13029 27331 33495 32341 37665 49117
#rowSums(cured_asvs)
pseqtest <- phyloseq(cured_asvs, taxtab, metadata)
#extract easy-to-use sample data
samdftest <- data.frame(sample_data(pseqtest))
sum(rowSums(cured_asvs))
## [1] 1293658
ncol(cured_asvs)
## [1] 9528
#NMDS of 16s, rarefied data rarefaction
summary(rowSums(cured_asvs))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13029 27331 33495 32341 37665 49117
sum(rowSums(cured_asvs))
## [1] 1293658
rowSums(cured_asvs)
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4 InCGS1-5 InCGS1-6 InCGS1-7
## 24472 34097 23968 30479 19431 33452 33634 25889
## InCGS1-8 InCGS1-9 InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4 InCGS2-5
## 27511 30799 43593 34118 45095 29467 25473 37620
## InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9 OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3
## 29426 22248 35118 14442 33537 25080 33214 26792
## OuCGS1-4 OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9 OuCGS2-1 OuCGS2-10
## 38159 31681 33266 35123 38570 27788 37799 33882
## OuCGS2-2 OuCGS2-3 OuCGS2-4 OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## 13029 43110 43611 35740 39800 49117 38285 35743
### NMDS no convergent solutions at the lower read depth.
set.seed(2024)
str(cured_asvs)
## Formal class 'otu_table' [package "phyloseq"] with 2 slots
## ..@ .Data : int [1:40, 1:9528] 177 213 109 307 89 139 251 148 103 137 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:40] "InCGS1-1" "InCGS1-10" "InCGS1-2" "InCGS1-3" ...
## .. .. ..$ : chr [1:9528] "ASV_1" "ASV_2" "ASV_3" "ASV_4" ...
## ..@ taxa_are_rows: logi FALSE
## ..$ dim : int [1:2] 40 9528
## ..$ dimnames:List of 2
## .. ..$ : chr [1:40] "InCGS1-1" "InCGS1-10" "InCGS1-2" "InCGS1-3" ...
## .. ..$ : chr [1:9528] "ASV_1" "ASV_2" "ASV_3" "ASV_4" ...
# ensure it's a data.frame
cured_asvs <- as.data.frame(cured_asvs)
# drop any non-numeric columns
cured_asvs <- cured_asvs[, sapply(cured_asvs, is.numeric)]
rowSums(cured_asvs)
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4 InCGS1-5 InCGS1-6 InCGS1-7
## 24472 34097 23968 30479 19431 33452 33634 25889
## InCGS1-8 InCGS1-9 InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4 InCGS2-5
## 27511 30799 43593 34118 45095 29467 25473 37620
## InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9 OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3
## 29426 22248 35118 14442 33537 25080 33214 26792
## OuCGS1-4 OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9 OuCGS2-1 OuCGS2-10
## 38159 31681 33266 35123 38570 27788 37799 33882
## OuCGS2-2 OuCGS2-3 OuCGS2-4 OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## 13029 43110 43611 35740 39800 49117 38285 35743
set.seed(2024)
rarefied <- rrarefy(cured_asvs, sample = 13029)
original_reads <- rowSums(cured_asvs)
rarefied_reads <- rowSums(rarefied)
discarded_reads <- original_reads - rarefied_reads
class(rarefied)
## [1] "matrix" "array"
str(rarefied)
## int [1:40, 1:9528] 88 79 54 141 54 57 103 73 47 55 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:40] "InCGS1-1" "InCGS1-10" "InCGS1-2" "InCGS1-3" ...
## ..$ : chr [1:9528] "ASV_1" "ASV_2" "ASV_3" "ASV_4" ...
# If rarefied is a matrix of ASV counts
# Convert to data frame
rarasv <- as.data.frame(rarefied)
# Remove ASVs (columns) that have zero counts across all samples
rarasv <- rarasv[, colSums(rarasv) > 0]
rarasv_1<-as.data.frame(t(rarasv)) #convert rarified asv count table to data frame and transpose it (ASV ID as rows, samples as columns)
rarasv_1 <- tibble::rownames_to_column(rarasv_1, "ASV_ID") #labels the first row of the first column (i.e., the column containing ASV IDs from row 2 onward)
write_tsv(rarasv_1, file="cogongrass_16s.tsv")
rarasv_2<- read.table("cogongrass_16s.tsv")
colnames(rarasv_2)<-rarasv_2[1,]
rarasv_3<-rarasv_2[-1,]
write_tsv(rarasv_3, file="cogongrass_16s.tsv")
# Summarize total read counts per sample after rarefaction
summary(rowSums(rarasv))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13029 13029 13029 13029 13029 13029
## look at read depths for the different substrates:
asvtab <- otu_table(rarasv, taxa_are_rows = FALSE)
taxtab <- tax_table(taxtab)
#combine into phyloseq object
pseqtest <- phyloseq(asvtab, taxtab, metadata)
#extract easy-to-use sample data
samdftest <- data.frame(sample_data(pseqtest))
##Beta diversity
#generate distance matrices
braytest <- phyloseq::distance(pseqtest, method = "bray")
#jactest <- phyloseq::distance(pseqtest, method="jaccard", binary=TRUE)
#perform NMDS
#Bray-Curtis
braymds <- metaMDS(braytest, k=2, trymax=500, wascores = T)
## Run 0 stress 0.1187049
## Run 1 stress 0.1191334
## ... Procrustes: rmse 0.01578087 max resid 0.0813796
## Run 2 stress 0.1190433
## ... Procrustes: rmse 0.0115116 max resid 0.06206553
## Run 3 stress 0.1701139
## Run 4 stress 0.1676428
## Run 5 stress 0.1196892
## Run 6 stress 0.1187053
## ... Procrustes: rmse 0.0009975231 max resid 0.004440977
## ... Similar to previous best
## Run 7 stress 0.1191338
## ... Procrustes: rmse 0.01597329 max resid 0.08191243
## Run 8 stress 0.119134
## ... Procrustes: rmse 0.01602776 max resid 0.08206076
## Run 9 stress 0.1187051
## ... Procrustes: rmse 0.0009498957 max resid 0.004229137
## ... Similar to previous best
## Run 10 stress 0.1191338
## ... Procrustes: rmse 0.01597972 max resid 0.08193144
## Run 11 stress 0.1187045
## ... New best solution
## ... Procrustes: rmse 0.000593883 max resid 0.002643882
## ... Similar to previous best
## Run 12 stress 0.1187054
## ... Procrustes: rmse 0.0004074497 max resid 0.001808442
## ... Similar to previous best
## Run 13 stress 0.1191335
## ... Procrustes: rmse 0.01623265 max resid 0.08266346
## Run 14 stress 0.1187051
## ... Procrustes: rmse 0.0006566119 max resid 0.002922768
## ... Similar to previous best
## Run 15 stress 0.1191336
## ... Procrustes: rmse 0.01629797 max resid 0.08283759
## Run 16 stress 0.1191335
## ... Procrustes: rmse 0.01601484 max resid 0.08206216
## Run 17 stress 0.1191337
## ... Procrustes: rmse 0.01634079 max resid 0.08295259
## Run 18 stress 0.1187054
## ... Procrustes: rmse 0.0004700561 max resid 0.002092734
## ... Similar to previous best
## Run 19 stress 0.1191334
## ... Procrustes: rmse 0.01603785 max resid 0.08213296
## Run 20 stress 0.1187052
## ... Procrustes: rmse 0.0002656048 max resid 0.00116298
## ... Similar to previous best
## *** Best solution repeated 5 times
## Warning in metaMDS(braytest, k = 2, trymax = 500, wascores = T): stress is
## (nearly) zero: you may have insufficient data
brayscores <- scores(braymds)
#Jaccard
#jacmds <- metaMDS(jactest, k=3, trymax=500, wascores = T)
#jacscores <- scores(jacmds)
braymds
##
## Call:
## metaMDS(comm = braytest, k = 2, trymax = 500, wascores = T)
##
## global Multidimensional Scaling using monoMDS
##
## Data: braytest
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1187045
## Stress type 1, weak ties
## Best solution was repeated 5 times in 20 tries
## The best solution was from try 11 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: scores missing
stressplot(braymds)
factorfit(brayscores, samdftest$treatment, permutations = 999)
## Centroids:
## NMDS1 NMDS2
## PInCG 0.0342 -0.0141
## POuCG -0.0342 0.0141
##
## Goodness of fit:
## r2 Pr(>r)
## P 0.0032 0.899
## Permutation: free
## Number of permutations: 999
myColors <- c("#9ACD32","#CD853F","#6495ED")
nmdsp3 <- plot_ordination(pseqtest, braymds, type="samples", color="treatment")+
geom_point(size=2)+scale_colour_manual(values=myColors)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the phyloseq package.
## Please report the issue at <https://github.com/joey711/phyloseq/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
nmdsp3
### PERMANOVA for treatment effects
adonis2(braytest ~ treatment, data = samdftest, permutations = 99999, by="terms")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 99999
##
## adonis2(formula = braytest ~ treatment, data = samdftest, permutations = 99999, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## treatment 1 0.3936 0.04563 1.8169 0.01278 *
## Residual 38 8.2323 0.95437
## Total 39 8.6259 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### testing for differences in dispersion between treatments
bds<- betadisper(braytest, samdftest$treatment)
anova(bds)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.014354 0.0143535 2.0467 0.1607
## Residuals 38 0.266488 0.0070128
no significant differences in dispersion between treatments.
### dispersion test for treatment
bds<- betadisper(braytest, samdftest$treatment)
permutest(bds)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.014354 0.0143535 2.0467 999 0.16
## Residuals 38 0.266488 0.0070128
plot(bds)
boxplot(bds)
mod.HSD<- TukeyHSD(bds)
mod.HSD
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## OuCG-InCG 0.03788607 -0.01572349 0.09149563 0.1607039
plot(mod.HSD)
#calculate alpha diversity
alpha <- estimate_richness(pseqtest, split=TRUE, measures=NULL)
row.names(alpha)<- row.names(samdftest)
alpha <- merge(samdftest, alpha, by= 'row.names', all=TRUE)
row.names(alpha) <- alpha$Row.names
pseq_richplot<- plot_richness(pseqtest, x = "treatment", color = "treatment", measures = c("Observed","Shannon"))+ geom_boxplot()
pseq_richplot
#Generalized linear models for alpha diversity
### checkinig for normality of response variables (required for gaussian)
# Shapiro-Wilk test for normality
shapiro_test <- alpha %>%
group_by(treatment) %>%
summarize(p_value = shapiro.test(Observed)$p.value)
shapiro_test
## # A tibble: 2 × 2
## treatment p_value
## <fct> <dbl>
## 1 InCG 0.975
## 2 OuCG 0.536
# Shapiro-Wilk test for normality
shapiro_test <- alpha %>%
group_by(treatment) %>%
summarize(p_value = shapiro.test(Shannon)$p.value)
shapiro_test
## # A tibble: 2 × 2
## treatment p_value
## <fct> <dbl>
## 1 InCG 0.261
## 2 OuCG 0.769
# Anova Alpha Diversity
alphaan<- aov(InvSimpson ~ treatment, data = alpha)
summary(alphaan)
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 5278 5278 0.268 0.608
## Residuals 38 748193 19689
alphawil<- wilcox.test(Observed ~ treatment, data = alpha)
alphawil
##
## Wilcoxon rank sum exact test
##
## data: Observed by treatment
## W = 225, p-value = 0.5117
## alternative hypothesis: true location shift is not equal to 0
#TAXONOMIC COMPOSITION:
#PHYLUM ~ class ~ family Taxonomy table
#make taxonomy object by phylum
phyla_counts_tab <- otu_table(tax_glom(pseqtest, taxrank="Phylum"), taxa_are_rows = FALSE)
phyla_counts_tab <- t(phyla_counts_tab)
#make vector of phyla names to set as row names
phyla_tax_vec <- as.vector(tax_table(tax_glom(pseqtest, taxrank="Phylum"))[,2])
rownames(phyla_counts_tab) <- as.vector(phyla_tax_vec)
phyla_counts_tab
## OTU Table: [32 taxa and 40 samples]
## taxa are rows
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4
## Pseudomonadota 3890 5279 3062 3627 3257
## Actinomycetota 4147 4207 5093 3993 4713
## Verrucomicrobiota 287 509 290 419 462
## Acidobacteriota 1809 1230 1230 2189 1577
## Bacteroidota 46 48 198 88 83
## Chloroflexota 1864 1029 2261 1492 1944
## RCP2-54 222 116 91 267 148
## Bacillota 242 243 327 279 333
## Candidatus Eremiobacterota 49 142 33 66 77
## Myxococcota 228 93 268 231 169
## Bdellovibrionota 52 40 26 63 22
## Gemmatimonadota 44 15 34 60 56
## Chlamydiota 20 26 11 24 37
## Thermodesulfobacteriota 18 7 29 42 37
## NB1-j 11 2 3 11 6
## Methylomirabilota 15 11 17 73 42
## Cyanobacteriota 5 11 9 5 8
## Planctomycetota 22 7 14 36 27
## Dependentiae 14 10 10 10 6
## Elusimicrobiota 0 1 2 0 3
## MBNT15 4 0 0 3 0
## Entotheonellaeota 0 0 4 13 0
## FCPU426 4 0 0 0 4
## GAL15 0 0 0 1 1
## Candidatus Kryptonia 0 0 5 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 0 0
## Patescibacteria 0 0 2 0 0
## Thermoproteota 0 0 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 0 0 2 0
## Deinococcota 0 0 0 0 0
## InCGS1-5 InCGS1-6 InCGS1-7 InCGS1-8 InCGS1-9
## Pseudomonadota 4407 4552 3739 4056 3796
## Actinomycetota 4861 3778 4663 2750 4102
## Verrucomicrobiota 278 490 427 567 422
## Acidobacteriota 1759 1982 1451 2733 2090
## Bacteroidota 17 40 166 29 55
## Chloroflexota 883 1263 1833 1831 1530
## RCP2-54 222 270 177 430 267
## Bacillota 155 109 126 70 305
## Candidatus Eremiobacterota 70 63 48 62 62
## Myxococcota 140 149 163 114 123
## Bdellovibrionota 19 79 44 23 32
## Gemmatimonadota 27 64 46 60 71
## Chlamydiota 51 43 10 19 33
## Thermodesulfobacteriota 18 25 35 123 16
## NB1-j 0 1 3 8 9
## Methylomirabilota 14 7 11 45 22
## Cyanobacteriota 12 15 15 3 17
## Planctomycetota 32 26 28 42 29
## Dependentiae 35 49 15 17 23
## Elusimicrobiota 8 5 3 0 1
## MBNT15 6 5 2 4 0
## Entotheonellaeota 0 0 0 5 3
## FCPU426 1 0 0 0 1
## GAL15 0 2 0 1 0
## Candidatus Kryptonia 2 2 2 8 4
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 0 0
## Patescibacteria 0 0 3 0 0
## Thermoproteota 3 0 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 0 0 0 3
## Deinococcota 0 0 0 0 0
## InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4
## Pseudomonadota 4538 4222 4389 3029 3047
## Actinomycetota 4039 4497 3658 4881 4607
## Verrucomicrobiota 218 277 558 385 783
## Acidobacteriota 1758 1565 2004 1493 1504
## Bacteroidota 39 95 20 144 41
## Chloroflexota 1424 1145 1426 1916 2121
## RCP2-54 232 164 289 95 171
## Bacillota 233 498 104 624 283
## Candidatus Eremiobacterota 56 143 116 89 159
## Myxococcota 181 187 222 217 121
## Bdellovibrionota 47 72 49 45 41
## Gemmatimonadota 29 28 24 23 37
## Chlamydiota 48 27 22 14 4
## Thermodesulfobacteriota 9 0 7 7 4
## NB1-j 2 0 1 0 4
## Methylomirabilota 3 1 5 0 9
## Cyanobacteriota 2 28 20 24 21
## Planctomycetota 43 28 50 18 8
## Dependentiae 87 24 45 7 35
## Elusimicrobiota 3 0 1 0 0
## MBNT15 0 1 0 0 0
## Entotheonellaeota 0 0 0 0 0
## FCPU426 0 3 1 3 0
## GAL15 0 0 0 0 4
## Candidatus Kryptonia 0 0 0 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 1 0 4 0 0
## Patescibacteria 0 3 0 0 13
## Thermoproteota 4 1 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 0 0 0 0
## Deinococcota 0 0 0 0 0
## InCGS2-5 InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9
## Pseudomonadota 4088 4671 5209 4119 5259
## Actinomycetota 3780 3356 4733 3676 4870
## Verrucomicrobiota 298 361 161 601 117
## Acidobacteriota 1586 1933 926 1759 1310
## Bacteroidota 65 74 270 79 81
## Chloroflexota 1801 1143 1101 1922 691
## RCP2-54 121 300 96 216 108
## Bacillota 648 661 284 112 88
## Candidatus Eremiobacterota 219 80 84 127 124
## Myxococcota 173 174 92 174 212
## Bdellovibrionota 68 89 18 95 48
## Gemmatimonadota 31 9 9 19 29
## Chlamydiota 5 20 10 15 5
## Thermodesulfobacteriota 1 5 0 14 0
## NB1-j 0 0 0 4 1
## Methylomirabilota 1 1 0 3 0
## Cyanobacteriota 22 32 4 20 14
## Planctomycetota 72 38 21 18 58
## Dependentiae 15 37 8 17 13
## Elusimicrobiota 0 6 0 10 0
## MBNT15 0 0 1 0 0
## Entotheonellaeota 0 0 0 0 0
## FCPU426 5 5 2 0 1
## GAL15 0 0 0 2 0
## Candidatus Kryptonia 0 12 0 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 0 0
## Patescibacteria 14 0 0 8 0
## Thermoproteota 0 1 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 2 0 0 0
## Deinococcota 0 0 0 0 0
## OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3 OuCGS1-4
## Pseudomonadota 4581 4187 4536 3604 4402
## Actinomycetota 4094 5117 4006 6147 4366
## Verrucomicrobiota 365 276 456 96 348
## Acidobacteriota 1881 970 1870 776 1720
## Bacteroidota 119 361 71 325 111
## Chloroflexota 944 1245 945 1235 1083
## RCP2-54 224 59 202 23 152
## Bacillota 259 238 432 279 234
## Candidatus Eremiobacterota 30 31 48 151 91
## Myxococcota 241 229 188 230 264
## Bdellovibrionota 54 68 90 38 87
## Gemmatimonadota 36 44 24 54 13
## Chlamydiota 33 0 37 1 19
## Thermodesulfobacteriota 48 36 11 1 12
## NB1-j 2 25 0 1 4
## Methylomirabilota 18 32 0 1 14
## Cyanobacteriota 15 29 23 7 19
## Planctomycetota 26 21 42 40 26
## Dependentiae 42 3 28 4 33
## Elusimicrobiota 8 1 6 0 6
## MBNT15 2 0 2 0 3
## Entotheonellaeota 0 14 0 0 0
## FCPU426 0 0 0 0 0
## GAL15 0 0 2 0 0
## Candidatus Kryptonia 0 0 0 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 0 0
## Patescibacteria 0 0 0 3 0
## Thermoproteota 0 0 0 0 0
## Fibrobacterota 0 3 0 0 0
## Candidatus Kapabacteria 0 0 0 0 0
## Deinococcota 0 0 0 0 0
## OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9
## Pseudomonadota 4357 4730 4184 4986 3618
## Actinomycetota 4259 4040 3470 3098 5469
## Verrucomicrobiota 328 403 175 396 268
## Acidobacteriota 1427 1616 867 2058 727
## Bacteroidota 301 114 1659 72 402
## Chloroflexota 1127 1145 898 1154 1073
## RCP2-54 140 194 72 273 82
## Bacillota 451 211 1179 194 265
## Candidatus Eremiobacterota 86 93 39 38 0
## Myxococcota 257 225 293 237 362
## Bdellovibrionota 44 74 52 84 60
## Gemmatimonadota 53 36 28 62 198
## Chlamydiota 15 23 6 96 12
## Thermodesulfobacteriota 41 22 24 37 93
## NB1-j 2 0 1 3 128
## Methylomirabilota 31 10 15 49 59
## Cyanobacteriota 23 25 5 16 10
## Planctomycetota 36 33 35 48 54
## Dependentiae 8 20 11 58 4
## Elusimicrobiota 6 0 0 9 8
## MBNT15 6 0 4 0 15
## Entotheonellaeota 6 3 2 14 45
## FCPU426 1 0 0 0 7
## GAL15 0 2 0 1 0
## Candidatus Kryptonia 0 0 0 0 0
## Latescibacterota 0 0 0 8 29
## Nitrospirota 4 0 1 0 6
## Patescibacteria 0 0 0 0 0
## Thermoproteota 0 0 0 1 0
## Fibrobacterota 0 0 0 0 18
## Candidatus Kapabacteria 0 0 0 3 0
## Deinococcota 0 0 4 0 0
## OuCGS2-1 OuCGS2-10 OuCGS2-2 OuCGS2-3 OuCGS2-4
## Pseudomonadota 4325 4596 4717 4943 4831
## Actinomycetota 4872 4090 4497 3941 4101
## Verrucomicrobiota 262 371 160 696 458
## Acidobacteriota 1441 1854 1367 1901 1606
## Bacteroidota 89 49 137 31 35
## Chloroflexota 1401 1297 1266 535 1201
## RCP2-54 109 248 96 374 161
## Bacillota 44 67 56 47 139
## Candidatus Eremiobacterota 162 138 149 28 134
## Myxococcota 127 133 324 249 171
## Bdellovibrionota 96 53 84 39 50
## Gemmatimonadota 12 25 29 12 20
## Chlamydiota 8 8 6 46 15
## Thermodesulfobacteriota 0 17 1 5 2
## NB1-j 0 1 3 8 2
## Methylomirabilota 0 0 0 4 5
## Cyanobacteriota 16 21 20 9 21
## Planctomycetota 46 51 95 61 30
## Dependentiae 9 6 8 93 28
## Elusimicrobiota 1 0 0 0 10
## MBNT15 0 0 0 0 0
## Entotheonellaeota 0 0 0 0 0
## FCPU426 0 1 0 0 0
## GAL15 0 0 0 0 0
## Candidatus Kryptonia 0 0 0 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 0 0
## Patescibacteria 3 2 9 2 0
## Thermoproteota 0 0 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 0 2 2 0
## Deinococcota 0 0 0 0 0
## OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## Pseudomonadota 5253 4189 5754 4372 5479
## Actinomycetota 4201 3482 2713 5793 3971
## Verrucomicrobiota 301 623 292 244 279
## Acidobacteriota 1960 1978 2644 1517 1488
## Bacteroidota 19 15 15 44 55
## Chloroflexota 703 1716 591 454 1058
## RCP2-54 227 314 373 106 151
## Bacillota 47 160 156 71 121
## Candidatus Eremiobacterota 43 125 15 128 105
## Myxococcota 143 208 197 100 157
## Bdellovibrionota 21 52 38 81 71
## Gemmatimonadota 11 26 29 12 12
## Chlamydiota 40 13 21 6 17
## Thermodesulfobacteriota 0 3 24 0 0
## NB1-j 0 3 27 0 0
## Methylomirabilota 0 11 14 0 1
## Cyanobacteriota 15 8 19 18 31
## Planctomycetota 29 33 34 38 17
## Dependentiae 12 45 37 36 9
## Elusimicrobiota 0 0 19 4 0
## MBNT15 4 0 1 0 0
## Entotheonellaeota 0 0 0 0 0
## FCPU426 0 0 4 2 0
## GAL15 0 0 0 0 0
## Candidatus Kryptonia 0 0 0 0 0
## Latescibacterota 0 0 0 0 0
## Nitrospirota 0 0 0 1 0
## Patescibacteria 0 3 1 0 2
## Thermoproteota 0 0 0 0 0
## Fibrobacterota 0 0 0 0 0
## Candidatus Kapabacteria 0 0 0 1 0
## Deinococcota 0 0 0 0 0
write.csv(phyla_counts_tab, 'phyla_counts_rar_05.12.2025.csv')
asv_counts <- pseqtest@otu_table
#determine the number of unclassified seqs at the phylum level
unclassified_tax_counts <- colSums(t(asv_counts)) - colSums(phyla_counts_tab)
#Add a row of "unclassified" to the phylum count table
phyla_and_unidentified_counts_tab <- rbind(phyla_counts_tab, "Unclassified_Phylum"=unclassified_tax_counts)
#split Pseudomonadota into classes
#remove Pseudomonadota for the purpose of separating it into classes later
temp_major_taxa_counts_tab <-
phyla_and_unidentified_counts_tab[!row.names(phyla_and_unidentified_counts_tab)
%in% "Pseudomonadota",]
#make count table broken down by class
class_counts_tab <- otu_table(tax_glom(t(pseqtest), taxrank="Class"))
#make table containing phylum and class data
class_tax_phy_tab <- tax_table(tax_glom(t(pseqtest), taxrank="Class"))
phy_tmp_vec <- class_tax_phy_tab[,2]
class_tmp_vec <- class_tax_phy_tab[,3]
rows_tmp <- row.names(class_tax_phy_tab)
class_tax_tab <- data.frame("Phylum"=phy_tmp_vec, "Class"=class_tmp_vec, row.names = rows_tmp)
#make vector of just the Pseudomonadota classes
pseudomonadata_classes_vec <- as.vector(class_tax_tab[class_tax_tab$Phylum == "Pseudomonadota", "Class"])
row.names(class_counts_tab) <- as.vector(class_tax_tab$Class)
#make table of Pseudomonadota classes
pseudomonadata_class_counts_tab <- class_counts_tab[row.names(class_counts_tab) %in% pseudomonadata_classes_vec, ]
#make table of pseudomonadata not identified to the class level
pseudomonadata_no_class_annotated_counts <-
phyla_and_unidentified_counts_tab[row.names(phyla_and_unidentified_counts_tab) %in% "Pseudomonadata",]-
colSums(pseudomonadata_class_counts_tab)
#Now combine the tables
major_taxa_counts_tab <- rbind(temp_major_taxa_counts_tab, pseudomonadata_class_counts_tab,
"Unclassified_Pseudomonadata"=pseudomonadata_no_class_annotated_counts)
head(major_taxa_counts_tab)
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4 InCGS1-5
## Actinomycetota 4147 4207 5093 3993 4713 4861
## Verrucomicrobiota 287 509 290 419 462 278
## Acidobacteriota 1809 1230 1230 2189 1577 1759
## Bacteroidota 46 48 198 88 83 17
## Chloroflexota 1864 1029 2261 1492 1944 883
## RCP2-54 222 116 91 267 148 222
## InCGS1-6 InCGS1-7 InCGS1-8 InCGS1-9 InCGS2-1 InCGS2-10
## Actinomycetota 3778 4663 2750 4102 4039 4497
## Verrucomicrobiota 490 427 567 422 218 277
## Acidobacteriota 1982 1451 2733 2090 1758 1565
## Bacteroidota 40 166 29 55 39 95
## Chloroflexota 1263 1833 1831 1530 1424 1145
## RCP2-54 270 177 430 267 232 164
## InCGS2-2 InCGS2-3 InCGS2-4 InCGS2-5 InCGS2-6 InCGS2-7
## Actinomycetota 3658 4881 4607 3780 3356 4733
## Verrucomicrobiota 558 385 783 298 361 161
## Acidobacteriota 2004 1493 1504 1586 1933 926
## Bacteroidota 20 144 41 65 74 270
## Chloroflexota 1426 1916 2121 1801 1143 1101
## RCP2-54 289 95 171 121 300 96
## InCGS2-8 InCGS2-9 OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3
## Actinomycetota 3676 4870 4094 5117 4006 6147
## Verrucomicrobiota 601 117 365 276 456 96
## Acidobacteriota 1759 1310 1881 970 1870 776
## Bacteroidota 79 81 119 361 71 325
## Chloroflexota 1922 691 944 1245 945 1235
## RCP2-54 216 108 224 59 202 23
## OuCGS1-4 OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9
## Actinomycetota 4366 4259 4040 3470 3098 5469
## Verrucomicrobiota 348 328 403 175 396 268
## Acidobacteriota 1720 1427 1616 867 2058 727
## Bacteroidota 111 301 114 1659 72 402
## Chloroflexota 1083 1127 1145 898 1154 1073
## RCP2-54 152 140 194 72 273 82
## OuCGS2-1 OuCGS2-10 OuCGS2-2 OuCGS2-3 OuCGS2-4 OuCGS2-5
## Actinomycetota 4872 4090 4497 3941 4101 4201
## Verrucomicrobiota 262 371 160 696 458 301
## Acidobacteriota 1441 1854 1367 1901 1606 1960
## Bacteroidota 89 49 137 31 35 19
## Chloroflexota 1401 1297 1266 535 1201 703
## RCP2-54 109 248 96 374 161 227
## OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## Actinomycetota 3482 2713 5793 3971
## Verrucomicrobiota 623 292 244 279
## Acidobacteriota 1978 2644 1517 1488
## Bacteroidota 15 15 44 55
## Chloroflexota 1716 591 454 1058
## RCP2-54 314 373 106 151
write.csv(major_taxa_counts_tab, "majortaxa_count_rar_5.12.2025.csv")
#Check that all sequences are accounted for
identical(colSums(major_taxa_counts_tab), rowSums(asv_counts))
## [1] TRUE
[1] TRUE
#Convert totals to relative abundance
major_taxa_proportions_tab <- apply(major_taxa_counts_tab, 2, function(x) x/sum(x)*100)
#colSums(major_taxa_proportions_tab)
write.csv(major_taxa_proportions_tab, "majortaxa_relative abundance_rar_5.12.2025col.csv")
### Taxa Barplot
phylatrim <- as.data.frame(major_taxa_proportions_tab[c("Actinomycetota","Verrucomicrobiota","Acidobacteriota","Alphaproteobacteria","Gammaproteobacteria","Bacteroidota","Chloroflexota", "Bacillota","Myxococcota", "Bdellovibrionota", "Gemmatimonadota"),])
#Create an "other" category for the phyla not retained
filtered_proportions <- colSums(major_taxa_proportions_tab) -
colSums(phylatrim)
phylatrim <- rbind(phylatrim, "Other"=filtered_proportions)
phylatrim
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4
## Actinomycetota 31.8289969 32.2895080 39.0897229 30.6470182 36.1731522
## Verrucomicrobiota 2.2027784 3.9066697 2.2258040 3.2159030 3.5459360
## Acidobacteriota 13.8844117 9.4404789 9.4404789 16.8009824 12.1037685
## Alphaproteobacteria 26.8861770 35.6589147 21.2372400 24.0693837 22.4959705
## Gammaproteobacteria 2.9702970 4.8583928 2.2641799 3.7685164 2.5021107
## Bacteroidota 0.3530586 0.3684089 1.5196869 0.6754164 0.6370404
## Chloroflexota 14.3065469 7.8977665 17.3535958 11.4513777 14.9205618
## Bacillota 1.8573950 1.8650702 2.5097859 2.1413769 2.5558370
## Myxococcota 1.7499424 0.7137923 2.0569499 1.7729680 1.2971065
## Bdellovibrionota 0.3991097 0.3070074 0.1995548 0.4835367 0.1688541
## Gemmatimonadota 0.3377082 0.1151278 0.2609563 0.4605112 0.4298104
## Other 3.2235782 2.5788625 1.8420447 4.5130094 3.1698519
## InCGS1-5 InCGS1-6 InCGS1-7 InCGS1-8 InCGS1-9
## Actinomycetota 37.3090797 28.9968532 35.7893929 21.1067618 31.4836135
## Verrucomicrobiota 2.1337017 3.7608412 3.2773045 4.3518305 3.2389285
## Acidobacteriota 13.5006524 15.2122189 11.1366951 20.9762837 16.0411390
## Alphaproteobacteria 29.1350065 29.8334485 25.6351217 28.0911812 24.1077596
## Gammaproteobacteria 4.6895387 5.1039988 3.0623993 3.0393737 5.0272469
## Bacteroidota 0.1304782 0.3070074 1.2740809 0.2225804 0.4221352
## Chloroflexota 6.7771893 9.6937601 14.0686162 14.0532658 11.7430348
## Bacillota 1.1896538 0.8365953 0.9670735 0.5372630 2.3409318
## Myxococcota 1.0745261 1.1436027 1.2510553 0.8749712 0.9440479
## Bdellovibrionota 0.1458285 0.6063397 0.3377082 0.1765293 0.2456060
## Gemmatimonadota 0.2072300 0.4912119 0.3530586 0.4605112 0.5449382
## Other 3.7071149 4.0141223 2.8474941 6.1094482 3.8606186
## InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4
## Actinomycetota 31.0000768 34.5153120 28.0758308 37.4625835 35.3595825
## Verrucomicrobiota 1.6731906 2.1260266 4.2827539 2.9549467 6.0096707
## Acidobacteriota 13.4929772 12.0116663 15.3810730 11.4590529 11.5434799
## Alphaproteobacteria 31.5066390 29.0275539 30.3783867 20.6539259 21.5442474
## Gammaproteobacteria 3.3233556 3.3770819 3.3080052 2.5942129 1.8420447
## Bacteroidota 0.2993323 0.7291427 0.1535037 1.1052268 0.3146826
## Chloroflexota 10.9294650 8.7880881 10.9448154 14.7056566 16.2790698
## Bacillota 1.7883184 3.8222427 0.7982194 4.7893161 2.1720777
## Myxococcota 1.3892087 1.4352598 1.7038913 1.6655154 0.9286975
## Bdellovibrionota 0.3607337 0.5526134 0.3760841 0.3453834 0.3146826
## Gemmatimonadota 0.2225804 0.2149052 0.1842045 0.1765293 0.2839819
## Other 4.0141223 3.4001075 4.4132320 2.0876506 3.4077826
## InCGS2-5 InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9
## Actinomycetota 29.0122035 25.75792463 36.32665592 28.2139842 37.3781564
## Verrucomicrobiota 2.2872055 2.77074219 1.23570497 4.6127869 0.8979968
## Acidobacteriota 12.1728452 14.83613478 7.10722235 13.5006524 10.0544938
## Alphaproteobacteria 28.2370097 31.93644946 34.49996162 28.9201013 36.9560212
## Gammaproteobacteria 3.1391511 3.91434492 5.48008289 2.6939903 3.4077826
## Bacteroidota 0.4988871 0.56796377 2.07230025 0.6063397 0.6216901
## Chloroflexota 13.8230102 8.77273774 8.45037992 14.7517077 5.3035536
## Bacillota 4.9735206 5.07329803 2.17975286 0.8596208 0.6754164
## Myxococcota 1.3278072 1.33548239 0.70611712 1.3354824 1.6271395
## Bdellovibrionota 0.5219127 0.68309156 0.13815335 0.7291427 0.3684089
## Gemmatimonadota 0.2379308 0.06907668 0.06907668 0.1458285 0.2225804
## Other 3.7685164 4.28275386 1.73459206 3.6303630 2.4867603
## OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3 OuCGS1-4
## Actinomycetota 31.4222120 39.2739274 30.7467956 47.1793691 33.50986261
## Verrucomicrobiota 2.8014429 2.1183514 3.4998849 0.7368179 2.67096477
## Acidobacteriota 14.4370251 7.4449305 14.3525981 5.9559444 13.20132013
## Alphaproteobacteria 31.1228797 27.7611482 31.4145368 26.1800599 30.30163481
## Gammaproteobacteria 4.0371479 4.3748561 3.4001075 1.4813109 3.48453450
## Bacteroidota 0.9133471 2.7707422 0.5449382 2.4944355 0.85194566
## Chloroflexota 7.2453757 9.5556067 7.2530509 9.4788549 8.31222657
## Bacillota 1.9878732 1.8266943 3.3156804 2.1413769 1.79599355
## Myxococcota 1.8497199 1.7576176 1.4429350 1.7652928 2.02624914
## Bdellovibrionota 0.4144601 0.5219127 0.6907668 0.2916571 0.66774119
## Gemmatimonadota 0.2763067 0.3377082 0.1842045 0.4144601 0.09977742
## Other 3.4922097 2.2565047 3.1545015 1.8804206 3.07774964
## OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9
## Actinomycetota 32.6886177 31.0077519 26.6328958 23.7777266 41.9755929
## Verrucomicrobiota 2.5174610 3.0931000 1.3431576 3.0393737 2.0569499
## Acidobacteriota 10.9524906 12.4031008 6.6543864 15.7955330 5.5798603
## Alphaproteobacteria 30.2555837 33.5482385 26.8708266 33.7861693 24.2996393
## Gammaproteobacteria 3.1852022 2.7553918 5.2421521 4.4823087 3.4691841
## Bacteroidota 2.3102310 0.8749712 12.7331338 0.5526134 3.0854248
## Chloroflexota 8.6499348 8.7880881 6.8923171 8.8571648 8.2354747
## Bacillota 3.4615089 1.6194643 9.0490444 1.4889861 2.0339243
## Myxococcota 1.9725228 1.7269169 2.2488295 1.8190191 2.7784174
## Bdellovibrionota 0.3377082 0.5679638 0.3991097 0.6447156 0.4605112
## Gemmatimonadota 0.4067849 0.2763067 0.2149052 0.4758615 1.5196869
## Other 3.2619541 3.3387060 1.7192417 5.2805281 4.5053343
## OuCGS2-1 OuCGS2-10 OuCGS2-2 OuCGS2-3 OuCGS2-4
## Actinomycetota 37.39350679 31.3915112 34.5153120 30.24790851 31.4759383
## Verrucomicrobiota 2.01089876 2.8474941 1.2280298 5.34192954 3.5152352
## Acidobacteriota 11.05994320 14.2297951 10.4919794 14.59052882 12.3263489
## Alphaproteobacteria 29.29618543 32.6425666 32.4813877 34.99884872 34.1008519
## Gammaproteobacteria 3.89899455 2.6325888 3.7224653 2.93959629 2.9779722
## Bacteroidota 0.68309156 0.3760841 1.0515005 0.23793077 0.2686315
## Chloroflexota 10.75293576 9.9547164 9.7167856 4.10622458 9.2178985
## Bacillota 0.33770819 0.5142375 0.4298104 0.36073375 1.0668509
## Myxococcota 0.97474864 1.0207998 2.4867603 1.91112134 1.3124568
## Bdellovibrionota 0.73681787 0.4067849 0.6447156 0.29933226 0.3837593
## Gemmatimonadota 0.09210223 0.1918797 0.2225804 0.09210223 0.1535037
## Other 2.76306700 3.7915419 3.0086730 4.87374319 3.2005526
## OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## Actinomycetota 32.24345690 26.7249981 20.8227800 44.46235321 30.47816410
## Verrucomicrobiota 2.31023102 4.7816410 2.2411543 1.87274541 2.14137693
## Acidobacteriota 15.04336480 15.1815182 20.2931921 11.64325735 11.42067695
## Alphaproteobacteria 36.20385294 29.0505795 39.0052959 30.28628444 38.23010208
## Gammaproteobacteria 4.11389976 3.1007752 5.1577251 3.26962929 3.82224269
## Bacteroidota 0.14582854 0.1151278 0.1151278 0.33770819 0.42213524
## Chloroflexota 5.39565584 13.1706194 4.5360350 3.48453450 8.12034692
## Bacillota 0.36073375 1.2280298 1.1973290 0.54493821 0.92869752
## Myxococcota 1.09755162 1.5964387 1.5120117 0.76751861 1.20500422
## Bdellovibrionota 0.16117891 0.3991097 0.2916571 0.62169008 0.54493821
## Gemmatimonadota 0.08442705 0.1995548 0.2225804 0.09210223 0.09210223
## Other 2.83981887 4.4516080 4.6051117 2.61723847 2.59421291
#Add taxa names as a column
phylatrim$Major_Taxa <- row.names(phylatrim)
#transform into long format
phylalong <- gather(phylatrim, Sample, Proportion, -Major_Taxa)
phylalong$Sample <- gsub("X","",phylalong$Sample)
phylamet<-data.frame("Sample"=row.names(samdftest),
"treatment"=samdftest$treatment,
stringsAsFactors=T)
#merge metadata with major taxa data
phylalong <- merge(phylalong, phylamet)
#Summarize by depth and hydration
phyla_summary <-
phylalong %>% # the names of the new data frame and the data frame to be summarised
group_by(treatment, Major_Taxa) %>% # the grouping variable
summarise(mean_prop = mean(Proportion)) # calculates the mean of each group
## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.
pal12 <- rev(cols25(n=12)) #set n= to the number of taxa you are plotting
#make stacked bar plot
phylum_bar <- ggplot(phyla_summary, aes(x = treatment, y = mean_prop, fill = Major_Taxa))+
geom_bar(stat = "identity", col=I("black")) +
scale_fill_manual(values=pal12)+
guides(fill=guide_legend(ncol=1))+
#facet_wrap(~treatment, labeller = labeller(Substrate=sublab), nrow=1, scales="free_x") +
labs(x='Treatment',y="Relative abundance (%)",
fill="Major taxa")+
#scale_facet_discrete(labels=c('B'='Epilithic biofilms', 'L'='Leaf litter', 'S'='Benthic sediment', 'W'='Surface water'))+
theme(axis.text.x = element_text(angle=60, hjust=1), legend.position = "right")+
theme(text=element_text(size=18), #change font size of all text
axis.text=element_text(size=16), #change font size of axis text
axis.title=element_text(size=16), #change font size of axis titles
plot.title=element_text(size=16), #change font size of plot title
legend.text=element_text(size=16), #change font size of legend text
strip.text.x = element_text(size = 16),
legend.title=element_text(size=16)) #change font size of legend title
phylum_bar
#by family
#make taxonomy object by genus
family_counts_tab <- otu_table(tax_glom(pseqtest, taxrank="Family"), taxa_are_rows = FALSE)
family_counts_tab <- t(family_counts_tab)
#make vector of genus names to set as row names
family_tax_vec <- as.vector(tax_table(tax_glom(pseqtest, taxrank="Family"))[,5])
rownames(family_counts_tab) <- as.vector(family_tax_vec)
#determine the number of unclassified seqs at the family level
unclassified_family_counts <- colSums(t(asv_counts)) - colSums(family_counts_tab)
unclassified_family_counts
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4 InCGS1-5 InCGS1-6 InCGS1-7
## 3468 1998 3693 4072 3454 3102 3440 3540
## InCGS1-8 InCGS1-9 InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4 InCGS2-5
## 4913 3958 3047 2707 3550 2533 3299 3052
## InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9 OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3
## 2725 1865 3404 2058 3411 3452 3023 2555
## OuCGS1-4 OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9 OuCGS2-1 OuCGS2-10
## 2998 3092 2913 2206 3521 4575 2340 3276
## OuCGS2-2 OuCGS2-3 OuCGS2-4 OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## 2637 3250 2585 2709 3584 3723 2238 2132
#Add a row of "unclassified" to the family count table
family_and_unidentified_counts_tab <- rbind(family_counts_tab,
"Unclassified_family"=unclassified_family_counts)
write.csv(family_and_unidentified_counts_tab, "family raw abundance_rar.csv")
#Check that all seqs are accounted for.
identical(colSums(family_and_unidentified_counts_tab), rowSums(asv_counts))
## [1] TRUE
#Convert totals to relative abundance
family_proportions_tab <- apply(family_and_unidentified_counts_tab, 2, function(x) x/sum(x)*100)
write.csv(family_proportions_tab, "family relative abundance_rar.csv")
#Merge metadata, phylum, and family abundances
#phy_and_fam <- merge(t(major_taxa_counts_tab), t(family_and_unidentified_counts_tab),
# by="row.names", all=TRUE)
phy_and_fam <- merge(t(major_taxa_proportions_tab), t(family_proportions_tab),
by="row.names", all=TRUE)
phyfam <- phy_and_fam[,-1]
rownames(phyfam) <- phy_and_fam[,1]
tax_merge <- merge(samdftest, phyfam,
by="row.names", all=TRUE)
taxmerge <- tax_merge[,-1]
rownames(taxmerge) <- tax_merge[,1]
write.csv(tax_merge, "cogon_5.12.2025_taxmerge.csv")
#summary stats about major taxa
colnames(taxmerge)
## [1] "dada2_input" "filtered"
## [3] "dada_f" "dada_r"
## [5] "merged" "nonchim"
## [7] "final_perc_reads_retained" "treatment"
## [9] "Actinomycetota" "Verrucomicrobiota"
## [11] "Acidobacteriota" "Bacteroidota"
## [13] "Chloroflexota" "RCP2-54"
## [15] "Bacillota" "Candidatus Eremiobacterota"
## [17] "Myxococcota" "Bdellovibrionota"
## [19] "Gemmatimonadota" "Chlamydiota"
## [21] "Thermodesulfobacteriota" "NB1-j"
## [23] "Methylomirabilota" "Cyanobacteriota"
## [25] "Planctomycetota" "Dependentiae"
## [27] "Elusimicrobiota" "MBNT15"
## [29] "Entotheonellaeota" "FCPU426"
## [31] "GAL15" "Candidatus Kryptonia"
## [33] "Latescibacterota" "Nitrospirota"
## [35] "Patescibacteria" "Thermoproteota"
## [37] "Fibrobacterota" "Candidatus Kapabacteria"
## [39] "Deinococcota" "Unclassified_Phylum"
## [41] "Alphaproteobacteria" "Gammaproteobacteria"
## [43] "Caulobacteraceae" "Mycobacteriaceae"
## [45] "Solirubrobacteraceae" "Incertae Sedis"
## [47] "Xanthobacteraceae" "Xiphinematobacteraceae"
## [49] "Unknown Family" "Acetobacteraceae"
## [51] "Chitinophagaceae" "Ktedonobacteraceae"
## [53] "Acidothermaceae" "Beijerinckiaceae"
## [55] "Koribacteraceae" "JG30-KF-AS9"
## [57] "Solibacteraceae" "Gaiellaceae"
## [59] "Acidimicrobiaceae" "Acidobacteriaceae (Subgroup 1)"
## [61] "Bacillaceae" "Burkholderiaceae"
## [63] "Sphingomonadaceae" "Bryobacteraceae"
## [65] "Alicyclobacillaceae" "A21b"
## [67] "Sporolactobacillaceae" "Geodermatophilaceae"
## [69] "Anaeromyxobacteraceae" "Polyangiaceae"
## [71] "Elsteraceae" "Streptomycetaceae"
## [73] "Dongiaceae" "Sutterellaceae"
## [75] "Labraceae" "Actinospicaceae"
## [77] "Micromonosporaceae" "Rhodomicrobiaceae"
## [79] "Frankiaceae" "Reyranellaceae"
## [81] "Paenibacillaceae" "Inquilinaceae"
## [83] "Hyphomicrobiaceae" "Planococcaceae"
## [85] "Micropepsaceae" "Pseudonocardiaceae"
## [87] "Chthoniobacteraceae" "Magnetospirillaceae"
## [89] "Gemmatimonadaceae" "Oligoflexaceae"
## [91] "Rhizobiaceae" "Parachlamydiaceae"
## [93] "67-14" "cvE6"
## [95] "Thermomonosporaceae" "Methyloligellaceae"
## [97] "Salisediminibacteriaceae" "Haliangiaceae"
## [99] "Amb-16S-1323" "Nocardioidaceae"
## [101] "WX65" "Myxococcaceae"
## [103] "Nocardiaceae" "Obscuribacteraceae"
## [105] "Thermoanaerobaculaceae" "WD2101 soil group"
## [107] "URHD0088" "Cytophagaceae"
## [109] "Pedosphaeraceae" "Microbacteriaceae"
## [111] "Lysobacteraceae" "KF-JG30-B3"
## [113] "Catenulisporaceae" "Comamonadaceae"
## [115] "Kineosporiaceae" "Nitrosomonadaceae"
## [117] "Vermiphilaceae" "Rhodanobacteraceae"
## [119] "Babeliaceae" "Hyphomicrobiales Incertae Sedis"
## [121] "Sphingobacteriaceae" "A0839"
## [123] "Pyrinomonadaceae" "Oxalobacteraceae"
## [125] "Rickettsiellaceae" "Incertae Sedis.1"
## [127] "Geminicoccaceae" "Enterobacteriaceae"
## [129] "Micrococcaceae" "Flavobacteriaceae"
## [131] "Streptosporangiaceae" "Nevskiaceae"
## [133] "JG30-KF-CM45" "Hymenobacteraceae"
## [135] "Thermoactinomycetaceae" "Bacteriovoracaceae"
## [137] "Pseudobdellovibrionaceae" "Vampirovibrionaceae"
## [139] "Vicinamibacteraceae" "Weeksellaceae"
## [141] "Intrasporangiaceae" "Simkaniaceae"
## [143] "Phaselicystaceae" "Caldilineaceae"
## [145] "Legionellaceae" "Blastocatellaceae"
## [147] "Pseudomonadaceae" "Methylacidiphilaceae"
## [149] "Entotheonellaceae" "Peptostreptococcaceae"
## [151] "Defluviicoccaceae" "BSV26"
## [153] "Chlamydiaceae" "Erwiniaceae"
## [155] "Opitutaceae" "Nostocaceae"
## [157] "Anaerolineaceae" "Sandaracinaceae"
## [159] "Isosphaeraceae" "Propionibacteriaceae"
## [161] "Sporichthyaceae" "Nitrospiraceae"
## [163] "Symbiobacteraceae" "Coleofasciculaceae"
## [165] "Ilumatobacteraceae" "27F-1492R"
## [167] "Spirosomataceae" "Iamiaceae"
## [169] "env.OPS 17" "BIrii41"
## [171] "Rhodospirillaceae" "Mycoplasmataceae"
## [173] "Parvibaculaceae" "Coxiellaceae"
## [175] "Fibrobacteraceae" "Nakamurellaceae"
## [177] "Vulgatibacteraceae" "Mitochondria"
## [179] "Sporomusaceae" "TRA3-20"
## [181] "Clostridiaceae" "Erysipelotrichaceae"
## [183] "AKYH767" "Promicromonosporaceae"
## [185] "Lactobacillaceae" "Azospirillaceae"
## [187] "Microscillaceae" "Woeseiaceae"
## [189] "Acidiferrobacteraceae" "Criblamydiaceae"
## [191] "37-13" "Longimicrobiaceae"
## [193] "Deinococcaceae" "Steroidobacteraceae"
## [195] "FFCH9454" "Devosiaceae"
## [197] "Terrimicrobiaceae" "Nitrosotaleaceae"
## [199] "Leptolyngbyaceae" "Oscillospiraceae"
## [201] "Unclassified_family"
# Read data
asv <- read.table("16S_dada2_ASVs_04.23.2025.count_table", sep = "\t", row.names = 1, header = TRUE)
meta1 <- read.csv("fin_read_track_16S_04162025.csv")
meta2 <- read.csv("16s cogon grass C N P.csv")
colnames(meta1)
## [1] "X" "dada2_input"
## [3] "filtered" "dada_f"
## [5] "dada_r" "merged"
## [7] "nonchim" "final_perc_reads_retained"
## [9] "treatment"
colnames(meta2)
## [1] "X" "Total.N.." "Total.C.." "Total.P.."
# Merge metadata
meta <- merge(meta1[, c("X", "treatment")], meta2, by = "X")
# Bray-Curtis dissimilarity
bray <- vegdist(asv, method = "bray")
table(meta$treatment)
##
## InCG OuCG
## 20 20
# Run PERMANOVA
adonis2(bray ~ treatment + `Total.N..` + `Total.C..` + `Total.P..`, data = meta, permutations = 999)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray ~ treatment + Total.N.. + Total.C.. + Total.P.., data = meta, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Model 4 1.3298 0.15411 1.5941 0.006 **
## Residual 35 7.2994 0.84589
## Total 39 8.6293 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(bray ~ treatment + Total.N.. + Total.C.. + Total.P.., data = meta, by = "margin")
## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray ~ treatment + Total.N.. + Total.C.. + Total.P.., data = meta, by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## treatment 1 0.4186 0.04850 2.0069 0.010 **
## Total.N.. 1 0.1406 0.01630 0.6743 0.879
## Total.C.. 1 0.4467 0.05176 2.1417 0.006 **
## Total.P.. 1 0.4399 0.05097 2.1091 0.009 **
## Residual 35 7.2994 0.84589
## Total 39 8.6293 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nmds <- metaMDS(asv, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1020735
## Run 1 stress 0.1020737
## ... Procrustes: rmse 0.0001165666 max resid 0.0005141103
## ... Similar to previous best
## Run 2 stress 0.1022275
## ... Procrustes: rmse 0.0050929 max resid 0.02568961
## Run 3 stress 0.1020737
## ... Procrustes: rmse 0.0005960203 max resid 0.002681769
## ... Similar to previous best
## Run 4 stress 0.1020743
## ... Procrustes: rmse 0.0003347005 max resid 0.001504045
## ... Similar to previous best
## Run 5 stress 0.1020734
## ... New best solution
## ... Procrustes: rmse 0.0004088529 max resid 0.00183878
## ... Similar to previous best
## Run 6 stress 0.102074
## ... Procrustes: rmse 0.0002894765 max resid 0.001299915
## ... Similar to previous best
## Run 7 stress 0.1022246
## ... Procrustes: rmse 0.004681454 max resid 0.02326084
## Run 8 stress 0.1020737
## ... Procrustes: rmse 0.0001800918 max resid 0.0008104657
## ... Similar to previous best
## Run 9 stress 0.102074
## ... Procrustes: rmse 0.0006067891 max resid 0.002730031
## ... Similar to previous best
## Run 10 stress 0.1020737
## ... Procrustes: rmse 0.0001519553 max resid 0.0006838149
## ... Similar to previous best
## Run 11 stress 0.1020736
## ... Procrustes: rmse 9.099708e-05 max resid 0.0004094525
## ... Similar to previous best
## Run 12 stress 0.1020734
## ... New best solution
## ... Procrustes: rmse 2.242329e-05 max resid 0.0001005555
## ... Similar to previous best
## Run 13 stress 0.1020739
## ... Procrustes: rmse 0.0005930036 max resid 0.002669431
## ... Similar to previous best
## Run 14 stress 0.1022269
## ... Procrustes: rmse 0.005047094 max resid 0.0254849
## Run 15 stress 0.1020737
## ... Procrustes: rmse 0.0005055131 max resid 0.002275027
## ... Similar to previous best
## Run 16 stress 0.1020737
## ... Procrustes: rmse 0.0001705018 max resid 0.0007648836
## ... Similar to previous best
## Run 17 stress 0.1020738
## ... Procrustes: rmse 0.0002463124 max resid 0.00110858
## ... Similar to previous best
## Run 18 stress 0.1020737
## ... Procrustes: rmse 0.0001720728 max resid 0.0007744372
## ... Similar to previous best
## Run 19 stress 0.1020738
## ... Procrustes: rmse 0.0002637384 max resid 0.001187033
## ... Similar to previous best
## Run 20 stress 0.1020737
## ... Procrustes: rmse 0.0005107592 max resid 0.002299596
## ... Similar to previous best
## *** Best solution repeated 8 times
plot(nmds, type = "n")
points(nmds, display = "sites", col = as.factor(meta$treatment), pch = 19)
legend("topright", legend = levels(as.factor(meta$treatment)), col = 1:2, pch = 19)
alpha <- estimate_richness(pseqtest, split=TRUE, measures=NULL)
rownames(meta)<- meta[,1]
meta<- meta[,-1]
row.names(alpha)<- row.names(meta)
alpha <- merge(meta, alpha, by= 'row.names', all=TRUE)
row.names(alpha)<- alpha[,1]
alpha<- alpha[,-1]
alpha.full<- lm(Observed ~ Total.P..* treatment + Total.P.. * treatment , data = alpha)
alpha.full.an<- anova(alpha.full)
alpha.full.an
## Analysis of Variance Table
##
## Response: Observed
## Df Sum Sq Mean Sq F value Pr(>F)
## Total.P.. 1 67416 67416 0.9742 0.3302
## treatment 1 25697 25697 0.3713 0.5461
## Total.P..:treatment 1 207367 207367 2.9967 0.0920 .
## Residuals 36 2491175 69199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p.alpha<- ggplot(alpha, aes(Total.P.., Observed, color = treatment))+
geom_point()+
geom_line()+
stat_smooth(method = lm)+
ggtitle("16S ASV Richness Over Time")
p.alpha
## `geom_smooth()` using formula = 'y ~ x'
treatanova<- aov(Total.P.. ~treatment, data = alpha)
summary(treatanova)
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 2.500e-08 2.500e-08 0.032 0.86
## Residuals 38 2.995e-05 7.882e-07
treatanova<- aov(Total.C.. ~treatment, data = alpha)
summary(treatanova)
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 1.756 1.7556 4.068 0.0508 .
## Residuals 38 16.401 0.4316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
treatanova<- aov(Total.N.. ~treatment, data = alpha)
summary(treatanova)
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.002723 0.0027225 3.489 0.0695 .
## Residuals 38 0.029655 0.0007804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###family Barplot
familytrim <- as.data.frame(family_proportions_tab[c("Solirubrobacteraceae","Xanthobacteraceae","Mycobacteriaceae","Acetobacteraceae","Solibacteraceae","Ktedonobacteraceae","Xiphinematobacteraceae", "Beijerinckiaceae","Bryobacteraceae", "Burkholderiaceae"),])
#Create an "other" category for the family not retained
filtered_familyproportions <- colSums(family_proportions_tab) -
colSums(familytrim)
familytrim <- rbind(familytrim, "Other"=filtered_familyproportions)
familytrim
## InCGS1-1 InCGS1-10 InCGS1-2 InCGS1-3 InCGS1-4
## Solirubrobacteraceae 15.849259 7.099547 20.4773966 12.3493745 20.9225574
## Xanthobacteraceae 12.050042 14.030240 9.1871978 10.3922020 10.9908665
## Mycobacteriaceae 4.543710 10.207998 2.3332566 3.8759690 2.5404866
## Acetobacteraceae 3.691765 7.636810 3.9987720 2.9012204 4.0217975
## Solibacteraceae 2.962622 1.435260 1.7729680 1.8343695 1.8113439
## Ktedonobacteraceae 2.770742 3.000998 2.9779722 1.8036687 1.9725228
## Xiphinematobacteraceae 1.665515 3.492210 1.4659605 2.0415995 2.6095633
## Beijerinckiaceae 1.335482 1.980198 0.8749712 1.0284749 1.1973290
## Bryobacteraceae 1.066851 1.181979 0.8903216 0.9210223 1.1743035
## Burkholderiaceae 0.660066 3.131476 0.4374856 1.2280298 0.4988871
## Other 53.403945 46.803285 55.5836979 61.6240694 52.2603423
## InCGS1-5 InCGS1-6 InCGS1-7 InCGS1-8 InCGS1-9
## Solirubrobacteraceae 22.0968608 14.291197 17.9829611 9.1795226 13.032466
## Xanthobacteraceae 14.5291273 15.020339 10.5764065 15.0510400 12.034692
## Mycobacteriaceae 4.9044439 5.472408 5.1577251 3.6917645 6.915343
## Acetobacteraceae 3.7301405 2.195103 4.0831990 1.2126794 1.742267
## Solibacteraceae 1.3968839 2.371633 1.9034462 2.8167933 2.072300
## Ktedonobacteraceae 1.5196869 2.402333 2.2334792 2.8014429 1.803669
## Xiphinematobacteraceae 1.8036687 3.269629 2.4100084 3.1545015 2.479085
## Beijerinckiaceae 1.5427124 1.243380 1.1743035 1.6962161 1.089876
## Bryobacteraceae 0.8903216 1.166628 1.2894313 1.1282524 1.151278
## Burkholderiaceae 1.4429350 1.297106 0.7291427 0.3530586 2.003224
## Other 46.1432190 51.270243 52.4598972 58.9147287 55.675800
## InCGS2-1 InCGS2-10 InCGS2-2 InCGS2-3 InCGS2-4
## Solirubrobacteraceae 12.4645023 17.5685010 11.3208995 19.3721698 23.3325658
## Xanthobacteraceae 18.2669430 16.0948653 15.9643871 11.1980966 10.9141147
## Mycobacteriaceae 7.1379231 3.9066697 6.0326963 1.5273620 4.1522757
## Acetobacteraceae 2.2181288 3.9834216 2.3409318 1.9955484 3.3233556
## Solibacteraceae 2.4330340 1.4122342 1.5350372 2.6172385 1.0207998
## Ktedonobacteraceae 4.2290276 3.8606186 3.8836442 6.8923171 5.6642874
## Xiphinematobacteraceae 1.3661831 1.5196869 3.8836442 2.7477166 5.3649551
## Beijerinckiaceae 0.9593983 0.9824238 1.9187965 0.8826464 1.2740809
## Bryobacteraceae 0.9286975 1.7652928 1.1743035 1.1589531 1.7115665
## Burkholderiaceae 0.6370404 0.5526134 0.3914345 0.5372630 0.3607337
## Other 49.3591220 48.3536726 51.5542252 51.0706885 42.8812649
## InCGS2-5 InCGS2-6 InCGS2-7 InCGS2-8 InCGS2-9
## Solirubrobacteraceae 13.063167 12.065393 18.3974211 13.0017653 15.818559
## Xanthobacteraceae 13.316448 17.783406 10.8834139 13.9304628 15.335022
## Mycobacteriaceae 4.781641 4.152276 4.9888710 5.7794152 7.176299
## Acetobacteraceae 5.249827 2.908896 7.3681787 2.8781948 9.348377
## Solibacteraceae 1.742267 2.632589 1.1512779 1.5120117 1.473636
## Ktedonobacteraceae 4.766291 4.344155 3.4077826 6.7388134 1.757618
## Xiphinematobacteraceae 1.903446 2.256505 1.1973290 4.1906516 0.782869
## Beijerinckiaceae 1.734592 1.872745 5.3803055 1.6348146 2.394658
## Bryobacteraceae 1.918797 1.404559 0.8058945 1.4122342 1.389209
## Burkholderiaceae 1.059176 1.097552 2.0646251 0.4374856 1.128252
## Other 50.464349 49.481925 44.3549006 48.4841507 43.395502
## OuCGS1-1 OuCGS1-10 OuCGS1-2 OuCGS1-3 OuCGS1-4
## Solirubrobacteraceae 10.8603884 8.5424822 12.0884181 24.0310078 10.906439
## Xanthobacteraceae 16.4402487 11.1980966 15.1584926 6.4625067 13.454601
## Mycobacteriaceae 4.6051117 4.4669583 6.1324737 4.6741883 6.278302
## Acetobacteraceae 1.9111213 3.6150127 3.3924323 11.7046588 3.791542
## Solibacteraceae 2.1567273 1.0822012 3.2005526 1.0668509 2.532811
## Ktedonobacteraceae 1.2357050 1.6117891 2.4483844 1.8804206 1.558063
## Xiphinematobacteraceae 1.7652928 1.2126794 2.9242459 0.5986645 2.103001
## Beijerinckiaceae 1.5273620 1.2203546 1.9571725 3.0393737 2.033924
## Bryobacteraceae 0.8135697 0.6293653 0.9593983 0.9286975 1.189654
## Burkholderiaceae 1.3738583 1.8497199 0.6754164 0.3914345 1.166628
## Other 57.3106148 64.5713409 51.0630133 45.2221966 54.985033
## OuCGS1-5 OuCGS1-6 OuCGS1-7 OuCGS1-8 OuCGS1-9
## Solirubrobacteraceae 10.5457057 9.056720 5.6949881 8.0052191 3.1084504
## Xanthobacteraceae 14.2144447 14.091642 7.7519380 16.1792923 8.7420370
## Mycobacteriaceae 3.9450457 7.421905 3.6764142 5.2191266 1.8650702
## Acetobacteraceae 3.6840893 4.489984 2.6249137 2.3793077 1.1896538
## Solibacteraceae 1.7422672 2.647939 0.8212449 2.5481618 0.5602886
## Ktedonobacteraceae 1.9341469 2.486760 1.1205772 2.0723003 0.3300330
## Xiphinematobacteraceae 1.5734132 2.509786 1.0131246 2.0492747 0.9824238
## Beijerinckiaceae 2.0415995 2.041600 5.7870903 1.8880958 1.1743035
## Bryobacteraceae 0.8903216 1.212679 0.5142375 0.9824238 0.2763067
## Burkholderiaceae 0.9900990 1.028475 2.2104536 1.3815335 0.8365953
## Other 58.4388671 53.012511 68.7850180 57.2952644 80.9348377
## OuCGS2-1 OuCGS2-10 OuCGS2-2 OuCGS2-3 OuCGS2-4
## Solirubrobacteraceae 14.736357 13.784634 16.6321283 12.5259038 11.2594980
## Xanthobacteraceae 11.696984 11.835137 11.0522680 17.4994244 14.2144447
## Mycobacteriaceae 9.125796 6.347379 4.0524983 9.6323586 10.5457057
## Acetobacteraceae 7.114898 6.339704 8.7343618 2.8398189 5.4110062
## Solibacteraceae 1.135928 1.565738 1.5657380 1.4429350 0.9517231
## Ktedonobacteraceae 7.421905 4.881418 2.4176836 1.9341469 5.0042214
## Xiphinematobacteraceae 1.673191 2.233479 0.8212449 4.8583928 3.3924323
## Beijerinckiaceae 2.195103 2.325581 3.2542789 2.2027784 3.3847571
## Bryobacteraceae 2.087651 1.366183 1.4045591 1.6348146 1.4045591
## Burkholderiaceae 1.235705 1.166628 1.2894313 0.8212449 0.9440479
## Other 41.576483 48.154118 48.7758078 44.6081817 43.4876046
## OuCGS2-5 OuCGS2-6 OuCGS2-7 OuCGS2-8 OuCGS2-9
## Solirubrobacteraceae 18.8963082 11.1059943 7.5447080 22.7722772 9.8319134
## Xanthobacteraceae 17.4687236 13.5697291 22.7415765 13.4546013 15.7955330
## Mycobacteriaceae 6.9920946 6.7464886 5.5107836 9.6323586 9.0336941
## Acetobacteraceae 5.5568348 2.5481618 2.1644025 6.9383683 8.4734055
## Solibacteraceae 1.5350372 1.4506102 2.0262491 0.9977742 1.5350372
## Ktedonobacteraceae 2.1183514 6.2322511 1.8880958 1.1205772 4.3364802
## Xiphinematobacteraceae 1.9187965 4.6281372 1.5657380 1.5887635 1.9725228
## Beijerinckiaceae 2.2104536 2.0492747 1.6194643 2.4560596 1.8190191
## Bryobacteraceae 0.8212449 1.0438253 1.3278072 1.6578402 0.8826464
## Burkholderiaceae 1.1819787 0.8903216 0.7751938 0.6140149 1.5657380
## Other 41.3001765 49.7352061 52.8359813 38.7673651 44.7540103
#Add taxa names as a column
familytrim$Major_Family <- row.names(familytrim)
#transform into long format
familylong <- gather(familytrim, Sample, Proportion, -Major_Family)
familylong$Sample <- gsub("X","",familylong$Sample)
familymet<-data.frame("Sample"=row.names(samdftest),
"treatment"=samdftest$treatment,
stringsAsFactors=T)
#merge metadata with major family data
familylong <- merge(familylong, familymet)
#Summarize by depth and hydration
family_summary <-
familylong %>% # the names of the new data frame and the data frame to be summarised
group_by(treatment, Major_Family) %>% # the grouping variable
summarise(mean_prop = mean(Proportion)) # calculates the mean of each group
## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.
pal10 <- rev(cols25(n=10)) #set n= to the number of taxa you are plotting
#make stacked bar plot
family_bar <- ggplot(family_summary, aes(x = treatment, y = mean_prop, fill = Major_Family))+
geom_bar(stat = "identity", col=I("black")) +
scale_fill_manual(values=pal12)+
guides(fill=guide_legend(ncol=1))+
#facet_wrap(~treatment, labeller = labeller(Substrate=sublab), nrow=1, scales="free_x") +
labs(x='Treatment',y="Relative abundance (%)",
fill="Major Family")+
#scale_facet_discrete(labels=c('B'='Epilithic biofilms', 'L'='Leaf litter', 'S'='Benthic sediment', 'W'='Surface water'))+
theme(axis.text.x = element_text(angle=60, hjust=1), legend.position = "right")+
theme(text=element_text(size=18), #change font size of all text
axis.text=element_text(size=16), #change font size of axis text
axis.title=element_text(size=16), #change font size of axis titles
plot.title=element_text(size=16), #change font size of plot title
legend.text=element_text(size=16), #change font size of legend text
strip.text.x = element_text(size = 16),
legend.title=element_text(size=16)) #change font size of legend title
family_bar
