Gene Marker Analysis Pipeline
Marjan Ghotbi
2025-04-23
1 Introduction
High-throughput Microvolume Extraction (HiMEx) methods were tested using: - Physical extraction through thermal shock (freeze-thaw, FT), - Combination of physical and chemical extraction (FT + proteinase K, FTP), - Further enhanced disruption with surfactant (FTP + IGEPAL, FTPIG).
2 prokaryotes amplicon analysis
Performs: Metadata cleanup and matching Rarefaction Beta diversity (ordination, betadisper) Differential abundance via DESeq2 Alpha diversity (Shannon, modeling with lme4)
3 Load packages
# Load libraries
# install.packages("devtools")
devtools::install_github("mghotbi/DspikeIn", build_vignettes = TRUE, dependencies = TRUE)## Skipping install of 'DspikeIn' from a github remote, the SHA1 (7fb23ccc) has not changed since last install.
## Use `force = TRUE` to force installation## starting httpd help server ...## done##
## Thank you for using DspikeIn! 🎉
## For support, 📧 contact: Mitra Ghotbi (mitra.ghotbi@gmail.com)
## GitHub Repository: https://github.com/mghotbi/DspikeIn
## 🐛 Found a bug or have a suggestion? Open an issue:
## https://github.com/mghotbi/DspikeIn/issues## Loading required package: permute## Loading required package: lattice## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: generics##
## Attaching package: 'generics'## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, union##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
## unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:phyloseq':
##
## distance## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## The following object is masked from 'package:phyloseq':
##
## sampleNames## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:S4Vectors':
##
## expand##
## Attaching package: 'lme4'## The following object is masked from 'package:generics':
##
## refit##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:SummarizedExperiment':
##
## coverage## The following object is masked from 'package:GenomicRanges':
##
## coverage## The following objects are masked from 'package:IRanges':
##
## coverage, transform## The following object is masked from 'package:S4Vectors':
##
## transform## The following object is masked from 'package:vegan':
##
## diversity## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:base':
##
## transform##
## Attaching package: 'dplyr'## The following object is masked from 'package:Biobase':
##
## combine## The following object is masked from 'package:matrixStats':
##
## count## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union## The following object is masked from 'package:GenomeInfoDb':
##
## intersect## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, setequal, union## The following object is masked from 'package:generics':
##
## explain## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.4 ✔ tidyr 1.3.1
## ✔ readr 2.1.5## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ microbiome::alpha() masks ggplot2::alpha()
## ✖ dplyr::collapse() masks IRanges::collapse()
## ✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count() masks matrixStats::count()
## ✖ dplyr::desc() masks IRanges::desc()
## ✖ tidyr::expand() masks Matrix::expand(), S4Vectors::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks S4Vectors::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ BiocGenerics::Position() masks ggplot2::Position(), base::Position()
## ✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename() masks S4Vectors::rename()
## ✖ lubridate::second() masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ dplyr::slice() masks IRanges::slice()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Registering fonts with R4 Load + prep data
# ───────────────────────────────────────────────────────────────
# Step 1: Load phyloseq objs
# ───────────────────────────────────────────────────────────────
library(here)## here() starts at /Users/mitraghotbi/Documents/MethodComparison5 QC & Contaminant Detection
6 Plot: Library Size per Sample (Colored by Method)
step 1: decontamination
# Custom color palette for Methods
method_palette <- c(
"conventional" = "#3D52A4",
"FTPIG" = "#A6D38D",
"FTP" = "#52B848",
"FT" = "#8CCCF0",
"Spike" = "#5C9DD5",
"blank" = "steelblue2",
"neg" = "red"
)
ggplot(df, aes(x = Index, y = LibrarySize, color = Method)) +
geom_point(size = 4) +
scale_color_manual(values = method_palette) +
theme_minimal() +
labs(
title = "Library Size per Sample",
x = "Sample Index",
y = "Library Size",
color = "Method"
)
7 Identify Blanks and Check Their Library Size
# Subset to blanks
blank <- subset_samples(physeq_16S, Sample_blank == "Blank")
# Recalculate and visualize
df2 <- microbiome::meta(blank)
df2$LibrarySize <- sample_sums(blank)
df2$SampleID <- rownames(df2)
df2 <- df2[order(df2$LibrarySize), ]
df2$Index <- seq_len(nrow(df2))
# Plot library sizes for blanks
ggplot(df2, aes(x = Index, y = LibrarySize, color = Sample_blank)) +
geom_point(size = 3) +
geom_text(aes(label = SampleID), hjust = 1.1, vjust = -0.5, size = 3) +
theme_minimal() +
labs(
title = "Library Size by Sample (Blanks)",
x = "Index",
y = "Library Size"
)
8 Identify Contaminants with decontam (Prevalence Method)
# Add logical vector for blanks to sample_data
sample_data(physeq_16S)$is.neg <- sample_data(physeq_16S)$Sample_blank == "Blank"
# Run decontam with default threshold
contamdf.prev <- isContaminant(physeq_16S, method = "prevalence", neg = "is.neg")## Warning in .is_contaminant(seqtab, conc = conc, neg = neg, method = method, :
## Removed 8 samples with zero total counts (or frequency).
## Warning in .is_contaminant(seqtab, conc = conc, neg = neg, method = method, :
## Removed 8 samples with zero total counts (or frequency).##
## FALSE TRUE
## 4619 35# Stricter threshold
contamdf.prev05 <- isContaminant(physeq_16S, method = "prevalence", neg = "is.neg", threshold = 0.05)## Warning in .is_contaminant(seqtab, conc = conc, neg = neg, method = method, :
## Removed 8 samples with zero total counts (or frequency).
## Warning in .is_contaminant(seqtab, conc = conc, neg = neg, method = method, :
## Removed 8 samples with zero total counts (or frequency).##
## FALSE TRUE
## 4625 29# Extract contaminant ASV IDs
contaminants <- rownames(contamdf.prev05[contamdf.prev05$contaminant == TRUE, ])
contaminants_prevalence <- contamdf.prev05[contamdf.prev05$contaminant == TRUE, ]
# Save list of contaminants
write.csv(contaminants_prevalence, "contaminants_prevalence_0.05_rarefaction4000.csv")9 Curated ASVs to Remove (Manual + Decontam + spiked bactertia removal)
10 spiked bacteria codes:
“1ae1bd31e0995f907062cc851d803bd7”,“39e27cbe03e2b71dbc8a5a1c8b8e7171”,“ac3a3e27aa0a24e16938ecd9de1c1060”,“1b468de6094a372d7fbc7d58b16f81c6”
ASVs_to_remove <- c(
# decontam + manual additions
"ac164649e33237f799c1ffee4f344d5d",
"92fb114641e27921122496428e2ef0dd",
"3a8a7409bdc33da14a01ec3e3313b9a7",
"3fcfe36ad281b523c2a5b9793512f5ad",
"863728e1cec6befd5ba02d15baef4c36",
"7c081f19c943c8d11819cdce94984719",
"1c16af05dfb4ec260717d24d0b9b8274",
"9076f0dc1a57c84a55fd88e445511903",
"1af47a3186d12f284691296661bc7310",
"f2d099ef5556f02539ded15181a3d994",
"fd55e4491e6b3cf2259827d9f78367d6",
"01cc2516697cf1b59d3bdcdc881c233e",
"637ddfd0682928e2e0c4e0001883425e",
"c728db0a39ffcd40e58dd84831880fde",
"e1aee885aa820fc1bbc8eea6c27cdc3d",
"ac70ecd361a0c67c771380f7aac46cc7",
"79f37fee0660e917bd1debe546718bad",
"4db5ac9fbbb89c5fc979475084c8c596",
"645dd8b575c7e7c0952933bbe90b4bb9",
"dd4ac7657b2665699cffd5a9b939a260",
"ac9b52e609fb05e7ede04ae1e6a6ff2f",
"e157cf6c27f73286aa20b0107e08e7d1",
"1ae1bd31e0995f907062cc851d803bd7",
"39e27cbe03e2b71dbc8a5a1c8b8e7171",
"ac3a3e27aa0a24e16938ecd9de1c1060",
"1b468de6094a372d7fbc7d58b16f81c6"
)
physeq_16S_clean <- prune_taxa(!taxa_names(physeq_16S) %in% ASVs_to_remove, physeq_16S)
physeq_16S_clean## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4628 taxa and 131 samples ]
## sample_data() Sample Data: [ 131 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 4628 taxa by 7 taxonomic ranks ]11 Step 2: Filter non-informative samples
# Start with samples that have at least one read
physeq2 <- prune_samples(sample_sums(physeq_16S_clean) > 0, physeq_16S_clean)
# Remove blanks
physeq3 <- subset_samples(physeq2, Type != "Blank")
# Remove synthetic spike-ins
physeq4 <-subset_samples(physeq3, Type != "Spike")
physeq4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4628 taxa and 84 samples ]
## sample_data() Sample Data: [ 84 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 4628 taxa by 7 taxonomic ranks ]## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4628 taxa and 84 samples ]
## sample_data() Sample Data: [ 84 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 4628 taxa by 7 taxonomic ranks ]## Sample_code
## GPT101a.0_2.500mL.04Oct22.spiked GPT101a-0_2-500mL-04Oct22-spiked
## GPT101b.0_2.500mL.04Oct22.spiked GPT101b-0_2-500mL-04Oct22-spiked
## GPT102a.0_2.500mL.07Oct22.spiked GPT102a-0_2-500mL-07Oct22-spiked
## GPT102b.0_2.500mL.07Oct22.spiked GPT102b-0_2-500mL-07Oct22-spiked
## GPT114a.0_2.500mL.15Nov22.spiked GPT114a-0_2-500mL-15Nov22-spiked
## GPT114c.0_2.500mL.15Nov22.spiked GPT114c-0_2-500mL-15Nov22-spiked
## GPT83a.0_2.500mL.05Aug22.spiked GPT83a-0_2-500mL-05Aug22-spiked
## GPT85b.0_2.500mL.12Aug22.spiked GPT85b-0_2-500mL-12Aug22-spiked
## GPT92a.0_2.500mL.02Sept22.spiked GPT92a-0_2-500mL-02Sept22-spiked
## GPT92c.0_2.500mL.02Sept22.spiked GPT92c-0_2-500mL-02Sept22-spiked
## Microvol.GPT100.30Sept22.MG.3 Microvol-GPT100-30Sept22-MG-3
## Microvol.GPT83.05Aug22.MG Microvol-GPT83-05Aug22-MG
## Microvol.GPT85.12Aug22.MG Microvol-GPT85-12Aug22-MG
## GPT101.02.500mL.4Oct22.spiked.ED GPT101-02-500mL-4Oct22-spiked-ED
## GPT102.02.500mL.7Oct22.spiked.ED GPT102-02-500mL-7Oct22-spiked-ED
## GPT83.02.500mL.05Aug22.spiked.ED GPT83-02-500mL-05Aug22-spiked-ED
## GPT85.02.500mL.12Aug22.spiked.ED GPT85-02-500mL-12Aug22-spiked-ED
## GPT92.02.500mL.02Sept22.spiked.ED GPT92-02-500mL-02Sept22-spiked-ED
## GPT.85.12.Aug.22 GPT-85-12-Aug-22
## mv.gpt83.1 gpt83-FT
## mv.gpt83.2 gpt83-FT
## mv.gpt83.3 gpt83-FTR
## mv.gpt83.4 gpt83-FTR
## mv.gpt83.5 gpt83-FTPR
## mv.gpt83.6 gpt83-FTPR
## mv.gpt85.1 gpt85-FT
## mv.gpt85.2 gpt85-FT
## mv.gpt85.3 gpt85-FTR
## mv.gpt85.4 gpt85-FTR
## mv.gpt85.5 gpt85-FTPR
## mv.gpt85.6 gpt85-FTPR
## gpt.101.13 gpt-101-1000-FTP
## gpt.101.14 gpt-101-400-FTP
## gpt.101.15 gpt-101-200-FTP
## gpt.101.16 gpt-101-100-FTP
## gpt.101.21 gpt-101-1000-FTPIG20
## gpt.101.22 gpt-101-400-FTPIG20
## gpt.101.23 gpt-101-200-FTPIG20
## gpt.101.24 gpt-101-100-FTPIG20
## gpt.101.5 gpt-101-1000-FT
## gpt.101.6 gpt-101-400-FT
## gpt.101.7 gpt-101-200-FT
## gpt.101.8 gpt-101-100-FT
## gpt.102.1 gpt-102-1000-FT
## gpt.102.11 gpt-102-200-FTP
## gpt.102.12 gpt-102-100-FTP
## gpt.102.17 gpt-102-1000-FTPIG20
## gpt.102.19 gpt-102-200-FTPIG20
## gpt.102.2 gpt-102-400-FT
## gpt.102.20 gpt-102-100-FTPIG20
## gpt.102.3 gpt-102-200-FT
## gpt.102.4 gpt-102-100-FT
## GPT100.0.2.1000mL.30Sept22.spiked.7 GPT100-0.2-50mL-30Sept22-spiked-7
## GPT100.0.2.50mL.30Sept22.spiked.1 GPT100-0.2-50mL-30Sept22-spiked-1
## GPT100.0.2.50mL.30Sept22.spiked.2 GPT100-0.2-50mL-30Sept22-spiked-2
## gpt100R.129 gpt100R-FTP
## gpt100R.130 gpt100R-FTP
## gpt100R.131 gpt100R-FTP
## gpt100R.132 gpt100R-FTP
## gpt100R.133 gpt100R-FTPIG
## gpt100R.134 gpt100R-FTPIG
## gpt100R.135 gpt100R-FTPIG
## gpt100R.136 gpt100R-FTPIG-sp104
## GPT101.0.2.500mL.04Oct22.14 GPT101-0.2-500mL-04Oct22-14
## GPT102.0.2.500mL.07Oct22.15 GPT102-0.2-500mL-07Oct22-15
## GPT114.0.2.3.500mL.15Nov22.spiked.56 GPT114-0.2-3-500mL-15Nov22-spiked-56
## GPT114.0.2.500mL.15Nov22.spiked.55 GPT114-0.2-500mL-15Nov22-spiked-55
## gpt114.137 gpt114-FTP
## gpt114.138 gpt114-FTP
## gpt114.139 gpt114-FTP
## gpt114.140 gpt114-FTP
## gpt114.141 gpt114-FTPIG
## gpt114.142 gpt114-FTPIG
## gpt114.143 gpt114-FTPIG
## gpt114.144 gpt114-FTPIG-sp104
## GPT83.0.2.500mL.5Aug22.9 GPT83-0.2-500mL-5Aug22-9
## gpt95R.121 gpt95R-FTP
## gpt95R.122 gpt95R-FTP
## gpt95R.123 gpt95R-FTP
## gpt95R.124 gpt95R-FTP
## gpt95R.125 gpt95R-FTPIG
## gpt95R.126 gpt95R-FTPIG
## gpt95R.127 gpt95R-FTPIG
## gpt95R.128 gpt95R-FTPIG-sp104
## Code Method
## GPT101a.0_2.500mL.04Oct22.spiked 101_con_1 conventional
## GPT101b.0_2.500mL.04Oct22.spiked 101_con_2 conventional
## GPT102a.0_2.500mL.07Oct22.spiked 102_con_1 conventional
## GPT102b.0_2.500mL.07Oct22.spiked 102_con_2 conventional
## GPT114a.0_2.500mL.15Nov22.spiked 114_con_1 conventional
## GPT114c.0_2.500mL.15Nov22.spiked 114_con_2 conventional
## GPT83a.0_2.500mL.05Aug22.spiked 83_Con_1 conventional
## GPT85b.0_2.500mL.12Aug22.spiked 85_con_1 conventional
## GPT92a.0_2.500mL.02Sept22.spiked 92_con_1 conventional
## GPT92c.0_2.500mL.02Sept22.spiked 92_con_2 conventional
## Microvol.GPT100.30Sept22.MG.3 100_con_3 conventional
## Microvol.GPT83.05Aug22.MG 83_Con_2 conventional
## Microvol.GPT85.12Aug22.MG 85_con_2 conventional
## GPT101.02.500mL.4Oct22.spiked.ED 101_con_3 conventional
## GPT102.02.500mL.7Oct22.spiked.ED 102_con_3 conventional
## GPT83.02.500mL.05Aug22.spiked.ED 83_Con_3 conventional
## GPT85.02.500mL.12Aug22.spiked.ED 85_Con_3 conventional
## GPT92.02.500mL.02Sept22.spiked.ED 92_Con_3 conventional
## GPT.85.12.Aug.22 85_Con_4 conventional
## mv.gpt83.1 83_FT_1 FT
## mv.gpt83.2 83_FT_2 FT
## mv.gpt83.3 83_FT_3 FT
## mv.gpt83.4 83_FTP_1 FTP
## mv.gpt83.5 83_FTP_2 FTP
## mv.gpt83.6 83_FTP_3 FTP
## mv.gpt85.1 85_FT_1 FT
## mv.gpt85.2 85_FT_2 FT
## mv.gpt85.3 85_FT_3 FT
## mv.gpt85.4 85_FTP_1 FTP
## mv.gpt85.5 85_FTP_2 FTP
## mv.gpt85.6 85_FTP_3 FTP
## gpt.101.13 101_FTP_1 FTP
## gpt.101.14 101_FTP_2 FTP
## gpt.101.15 101_FTP_3 FTP
## gpt.101.16 101_FTP_4 FTP
## gpt.101.21 101_FTPIG_1 FTPIG
## gpt.101.22 101_FTPIG_2 FTPIG
## gpt.101.23 101_FTPIG_3 FTPIG
## gpt.101.24 101_FTPIG_4 FTPIG
## gpt.101.5 101_FT_1 FT
## gpt.101.6 101_FT_2 FT
## gpt.101.7 101_FT_3 FT
## gpt.101.8 101_FT_4 FT
## gpt.102.1 102_FT_4 FT
## gpt.102.11 102_FTP_2 FTP
## gpt.102.12 102_FTP_3 FTP
## gpt.102.17 102_FTPIG_1 FTPIG
## gpt.102.19 102_FTPIG_3 FTPIG
## gpt.102.2 102_FT_1 FT
## gpt.102.20 102_FTPIG_4 FTPIG
## gpt.102.3 102_FT_2 FT
## gpt.102.4 102_FT_3 FT
## GPT100.0.2.1000mL.30Sept22.spiked.7 100_con_1 conventional
## GPT100.0.2.50mL.30Sept22.spiked.1 100_con_2 conventional
## GPT100.0.2.50mL.30Sept22.spiked.2 100_con_4 conventional
## gpt100R.129 100_FTP_1 FTP
## gpt100R.130 100_FTP_2 FTP
## gpt100R.131 100_FTP_3 FTP
## gpt100R.132 100_FTP_4 FTP
## gpt100R.133 100_FTPIG_1 FTPIG
## gpt100R.134 100_FTPIG_2 FTPIG
## gpt100R.135 100_FTPIG_3 FTPIG
## gpt100R.136 100_FTPIG_4_sp FTPIG
## GPT101.0.2.500mL.04Oct22.14 101_con_4 conventional
## GPT102.0.2.500mL.07Oct22.15 102_con_4 conventional
## GPT114.0.2.3.500mL.15Nov22.spiked.56 114_con_3 conventional
## GPT114.0.2.500mL.15Nov22.spiked.55 114_con_4 conventional
## gpt114.137 114_FTP_1 FTP
## gpt114.138 114_FTP_2 FTP
## gpt114.139 114_FTP_3 FTP
## gpt114.140 114_FTP_4 FTP
## gpt114.141 114_FTPIG_1 FTPIG
## gpt114.142 114_FTPIG_2 FTPIG
## gpt114.143 114_FTPIG_3 FTPIG
## gpt114.144 114_FTPIG_4_sp FTPIG
## GPT83.0.2.500mL.5Aug22.9 83_Con_4 conventional
## gpt95R.121 92_FTP_1 FTP
## gpt95R.122 92_FTP_2 FTP
## gpt95R.123 92_FTP_3 FTP
## gpt95R.124 92_FTP_4 FTP
## gpt95R.125 92_FTPIG_1 FTPIG
## gpt95R.126 92_FTPIG_2 FTPIG
## gpt95R.127 92_FTPIG_3 FTPIG
## gpt95R.128 92_FTPIG_4 FTPIG
## Extracted_volume Size Sample_blank
## GPT101a.0_2.500mL.04Oct22.spiked 50000 bulk Sample
## GPT101b.0_2.500mL.04Oct22.spiked 50000 bulk Sample
## GPT102a.0_2.500mL.07Oct22.spiked 50000 bulk Sample
## GPT102b.0_2.500mL.07Oct22.spiked 50000 bulk Sample
## GPT114a.0_2.500mL.15Nov22.spiked 50000 bulk Sample
## GPT114c.0_2.500mL.15Nov22.spiked 50000 bulk Sample
## GPT83a.0_2.500mL.05Aug22.spiked 50000 bulk Sample
## GPT85b.0_2.500mL.12Aug22.spiked 50000 bulk Sample
## GPT92a.0_2.500mL.02Sept22.spiked 50000 bulk Sample
## GPT92c.0_2.500mL.02Sept22.spiked 50000 bulk Sample
## Microvol.GPT100.30Sept22.MG.3 50000 bulk Sample
## Microvol.GPT83.05Aug22.MG 50000 bulk Sample
## Microvol.GPT85.12Aug22.MG 50000 bulk Sample
## GPT101.02.500mL.4Oct22.spiked.ED 50000 bulk Sample
## GPT102.02.500mL.7Oct22.spiked.ED 50000 bulk Sample
## GPT83.02.500mL.05Aug22.spiked.ED 50000 bulk Sample
## GPT85.02.500mL.12Aug22.spiked.ED 50000 bulk Sample
## GPT92.02.500mL.02Sept22.spiked.ED 50000 bulk Sample
## GPT.85.12.Aug.22 50000 bulk Sample
## mv.gpt83.1 1000 mv Sample
## mv.gpt83.2 1000 mv Sample
## mv.gpt83.3 1000 mv Sample
## mv.gpt83.4 1000 mv Sample
## mv.gpt83.5 1000 mv Sample
## mv.gpt83.6 1000 mv Sample
## mv.gpt85.1 1000 mv Sample
## mv.gpt85.2 1000 mv Sample
## mv.gpt85.3 1000 mv Sample
## mv.gpt85.4 1000 mv Sample
## mv.gpt85.5 1000 mv Sample
## mv.gpt85.6 1000 mv Sample
## gpt.101.13 1000 mv Sample
## gpt.101.14 400 mv Sample
## gpt.101.15 200 mv Sample
## gpt.101.16 100 mv Sample
## gpt.101.21 1000 mv Sample
## gpt.101.22 400 mv Sample
## gpt.101.23 200 mv Sample
## gpt.101.24 100 mv Sample
## gpt.101.5 1000 mv Sample
## gpt.101.6 400 mv Sample
## gpt.101.7 200 mv Sample
## gpt.101.8 100 mv Sample
## gpt.102.1 1000 mv Sample
## gpt.102.11 200 mv Sample
## gpt.102.12 100 mv Sample
## gpt.102.17 1000 mv Sample
## gpt.102.19 200 mv Sample
## gpt.102.2 400 mv Sample
## gpt.102.20 100 mv Sample
## gpt.102.3 200 mv Sample
## gpt.102.4 100 mv Sample
## GPT100.0.2.1000mL.30Sept22.spiked.7 50000 bulk Sample
## GPT100.0.2.50mL.30Sept22.spiked.1 50000 bulk Sample
## GPT100.0.2.50mL.30Sept22.spiked.2 50000 bulk Sample
## gpt100R.129 100 mv Sample
## gpt100R.130 100 mv Sample
## gpt100R.131 100 mv Sample
## gpt100R.132 100 mv Sample
## gpt100R.133 100 mv Sample
## gpt100R.134 100 mv Sample
## gpt100R.135 100 mv Sample
## gpt100R.136 100 mv Sample
## GPT101.0.2.500mL.04Oct22.14 50000 bulk Sample
## GPT102.0.2.500mL.07Oct22.15 50000 bulk Sample
## GPT114.0.2.3.500mL.15Nov22.spiked.56 50000 bulk Sample
## GPT114.0.2.500mL.15Nov22.spiked.55 50000 bulk Sample
## gpt114.137 100 mv Sample
## gpt114.138 100 mv Sample
## gpt114.139 100 mv Sample
## gpt114.140 100 mv Sample
## gpt114.141 100 mv Sample
## gpt114.142 100 mv Sample
## gpt114.143 100 mv Sample
## gpt114.144 100 mv Sample
## GPT83.0.2.500mL.5Aug22.9 50000 bulk Sample
## gpt95R.121 100 mv Sample
## gpt95R.122 100 mv Sample
## gpt95R.123 100 mv Sample
## gpt95R.124 100 mv Sample
## gpt95R.125 100 mv Sample
## gpt95R.126 100 mv Sample
## gpt95R.127 100 mv Sample
## gpt95R.128 100 mv Sample
## Purification Type KFTno Month Temp
## GPT101a.0_2.500mL.04Oct22.spiked conventional KFT 101 October 15
## GPT101b.0_2.500mL.04Oct22.spiked conventional KFT 101 October 15
## GPT102a.0_2.500mL.07Oct22.spiked conventional KFT 102 October 15
## GPT102b.0_2.500mL.07Oct22.spiked conventional KFT 102 October 15
## GPT114a.0_2.500mL.15Nov22.spiked conventional KFT 114 Novemver 12
## GPT114c.0_2.500mL.15Nov22.spiked conventional KFT 114 Novemver 12
## GPT83a.0_2.500mL.05Aug22.spiked conventional KFT 83 August 19
## GPT85b.0_2.500mL.12Aug22.spiked conventional KFT 85 August 21
## GPT92a.0_2.500mL.02Sept22.spiked conventional KFT 92 September 20
## GPT92c.0_2.500mL.02Sept22.spiked conventional KFT 92 September 20
## Microvol.GPT100.30Sept22.MG.3 conventional KFT 100 September 16
## Microvol.GPT83.05Aug22.MG conventional KFT 83 August 19
## Microvol.GPT85.12Aug22.MG conventional KFT 85 August 19
## GPT101.02.500mL.4Oct22.spiked.ED conventional KFT 101 October 15
## GPT102.02.500mL.7Oct22.spiked.ED conventional KFT 102 October 15
## GPT83.02.500mL.05Aug22.spiked.ED conventional KFT 83 August 19
## GPT85.02.500mL.12Aug22.spiked.ED conventional KFT 85 August 21
## GPT92.02.500mL.02Sept22.spiked.ED conventional KFT 92 September 20
## GPT.85.12.Aug.22 conventional KFT 85 August 21
## mv.gpt83.1 Direct KFT 83 August 19
## mv.gpt83.2 Direct KFT 83 August 19
## mv.gpt83.3 Direct KFT 83 August 19
## mv.gpt83.4 Direct KFT 83 August 19
## mv.gpt83.5 Direct KFT 83 August 19
## mv.gpt83.6 Direct KFT 83 August 19
## mv.gpt85.1 Direct KFT 85 August 21
## mv.gpt85.2 Direct KFT 85 August 21
## mv.gpt85.3 Direct KFT 85 August 21
## mv.gpt85.4 Direct KFT 85 August 21
## mv.gpt85.5 Direct KFT 85 August 21
## mv.gpt85.6 Direct KFT 85 August 21
## gpt.101.13 Direct KFT 101 October 15
## gpt.101.14 Direct KFT 101 October 15
## gpt.101.15 Direct KFT 101 October 15
## gpt.101.16 Direct KFT 101 October 15
## gpt.101.21 Direct KFT 101 October 15
## gpt.101.22 Direct KFT 101 October 15
## gpt.101.23 Direct KFT 101 October 15
## gpt.101.24 Direct KFT 101 October 15
## gpt.101.5 Direct KFT 101 October 15
## gpt.101.6 Direct KFT 101 October 15
## gpt.101.7 Direct KFT 101 October 15
## gpt.101.8 Direct KFT 101 October 15
## gpt.102.1 Direct KFT 102 October 15
## gpt.102.11 Direct KFT 102 October 15
## gpt.102.12 Direct KFT 102 October 15
## gpt.102.17 Direct KFT 102 October 15
## gpt.102.19 Direct KFT 102 October 15
## gpt.102.2 Direct KFT 102 October 15
## gpt.102.20 Direct KFT 102 October 15
## gpt.102.3 Direct KFT 102 October 15
## gpt.102.4 Direct KFT 102 October 15
## GPT100.0.2.1000mL.30Sept22.spiked.7 conventional KFT 100 September 16
## GPT100.0.2.50mL.30Sept22.spiked.1 conventional KFT 100 September 16
## GPT100.0.2.50mL.30Sept22.spiked.2 conventional KFT 100 September 16
## gpt100R.129 Direct KFT 100 September 16
## gpt100R.130 Direct KFT 100 September 16
## gpt100R.131 Direct KFT 100 September 16
## gpt100R.132 Direct KFT 100 September 16
## gpt100R.133 Direct KFT 100 September 16
## gpt100R.134 Direct KFT 100 September 16
## gpt100R.135 Direct KFT 100 September 16
## gpt100R.136 Direct KFT 100 September 16
## GPT101.0.2.500mL.04Oct22.14 conventional KFT 101 October 15
## GPT102.0.2.500mL.07Oct22.15 conventional KFT 102 October 15
## GPT114.0.2.3.500mL.15Nov22.spiked.56 conventional KFT 114 Novemver 12
## GPT114.0.2.500mL.15Nov22.spiked.55 conventional KFT 114 Novemver 12
## gpt114.137 Direct KFT 114 Novemver 12
## gpt114.138 Direct KFT 114 Novemver 12
## gpt114.139 Direct KFT 114 Novemver 12
## gpt114.140 Direct KFT 114 Novemver 12
## gpt114.141 Direct KFT 114 Novemver 12
## gpt114.142 Direct KFT 114 Novemver 12
## gpt114.143 Direct KFT 114 Novemver 12
## gpt114.144 Direct KFT 114 Novemver 12
## GPT83.0.2.500mL.5Aug22.9 conventional KFT 83 August 19
## gpt95R.121 Direct KFT 92 September 18
## gpt95R.122 Direct KFT 92 September 18
## gpt95R.123 Direct KFT 92 September 18
## gpt95R.124 Direct KFT 92 September 18
## gpt95R.125 Direct KFT 92 September 18
## gpt95R.126 Direct KFT 92 September 18
## gpt95R.127 Direct KFT 92 September 18
## gpt95R.128 Direct KFT 92 September 18
## Salinity is.neg
## GPT101a.0_2.500mL.04Oct22.spiked 15 FALSE
## GPT101b.0_2.500mL.04Oct22.spiked 15 FALSE
## GPT102a.0_2.500mL.07Oct22.spiked 9 FALSE
## GPT102b.0_2.500mL.07Oct22.spiked 9 FALSE
## GPT114a.0_2.500mL.15Nov22.spiked 21 FALSE
## GPT114c.0_2.500mL.15Nov22.spiked 21 FALSE
## GPT83a.0_2.500mL.05Aug22.spiked 16 FALSE
## GPT85b.0_2.500mL.12Aug22.spiked 16 FALSE
## GPT92a.0_2.500mL.02Sept22.spiked 15 FALSE
## GPT92c.0_2.500mL.02Sept22.spiked 15 FALSE
## Microvol.GPT100.30Sept22.MG.3 14 FALSE
## Microvol.GPT83.05Aug22.MG 16 FALSE
## Microvol.GPT85.12Aug22.MG 15 FALSE
## GPT101.02.500mL.4Oct22.spiked.ED 15 FALSE
## GPT102.02.500mL.7Oct22.spiked.ED 9 FALSE
## GPT83.02.500mL.05Aug22.spiked.ED 16 FALSE
## GPT85.02.500mL.12Aug22.spiked.ED 15 FALSE
## GPT92.02.500mL.02Sept22.spiked.ED 15 FALSE
## GPT.85.12.Aug.22 15 FALSE
## mv.gpt83.1 16 FALSE
## mv.gpt83.2 16 FALSE
## mv.gpt83.3 16 FALSE
## mv.gpt83.4 16 FALSE
## mv.gpt83.5 16 FALSE
## mv.gpt83.6 16 FALSE
## mv.gpt85.1 15 FALSE
## mv.gpt85.2 15 FALSE
## mv.gpt85.3 15 FALSE
## mv.gpt85.4 15 FALSE
## mv.gpt85.5 15 FALSE
## mv.gpt85.6 15 FALSE
## gpt.101.13 15 FALSE
## gpt.101.14 15 FALSE
## gpt.101.15 15 FALSE
## gpt.101.16 15 FALSE
## gpt.101.21 15 FALSE
## gpt.101.22 15 FALSE
## gpt.101.23 15 FALSE
## gpt.101.24 15 FALSE
## gpt.101.5 15 FALSE
## gpt.101.6 15 FALSE
## gpt.101.7 15 FALSE
## gpt.101.8 15 FALSE
## gpt.102.1 9 FALSE
## gpt.102.11 9 FALSE
## gpt.102.12 9 FALSE
## gpt.102.17 9 FALSE
## gpt.102.19 9 FALSE
## gpt.102.2 9 FALSE
## gpt.102.20 9 FALSE
## gpt.102.3 9 FALSE
## gpt.102.4 9 FALSE
## GPT100.0.2.1000mL.30Sept22.spiked.7 14 FALSE
## GPT100.0.2.50mL.30Sept22.spiked.1 14 FALSE
## GPT100.0.2.50mL.30Sept22.spiked.2 14 FALSE
## gpt100R.129 14 FALSE
## gpt100R.130 14 FALSE
## gpt100R.131 14 FALSE
## gpt100R.132 14 FALSE
## gpt100R.133 14 FALSE
## gpt100R.134 14 FALSE
## gpt100R.135 14 FALSE
## gpt100R.136 14 FALSE
## GPT101.0.2.500mL.04Oct22.14 15 FALSE
## GPT102.0.2.500mL.07Oct22.15 9 FALSE
## GPT114.0.2.3.500mL.15Nov22.spiked.56 21 FALSE
## GPT114.0.2.500mL.15Nov22.spiked.55 21 FALSE
## gpt114.137 21 FALSE
## gpt114.138 21 FALSE
## gpt114.139 21 FALSE
## gpt114.140 21 FALSE
## gpt114.141 21 FALSE
## gpt114.142 21 FALSE
## gpt114.143 21 FALSE
## gpt114.144 21 FALSE
## GPT83.0.2.500mL.5Aug22.9 16 FALSE
## gpt95R.121 15 FALSE
## gpt95R.122 15 FALSE
## gpt95R.123 15 FALSE
## gpt95R.124 15 FALSE
## gpt95R.125 15 FALSE
## gpt95R.126 15 FALSE
## gpt95R.127 15 FALSE
## gpt95R.128 15 FALSE12 compute the Spearman correlation
13 normalization
rle Tightly clustered points near y = 0: That shows consistent scaling across samples. Symmetry around the zero line: No bias toward high or low values in specific samples.
DESeq DESeq2 computes per-sample scaling factors that adjust for sequencing depth and compositional biases. Put samples on the same “scale” without rarefying, retaining more power. The distribution of values across samples is reasonably uniform. No major skews or samples standing out with inflated variance.
clr Values near 0 mean the taxon abundance is close to the geometric mean of that sample. Most of the values cluster close to 0: expected in CLR — the transformation centers data per sample. Some deep negative outliers -> taxa nearly absent in those samples (low relative abundance).
library(DspikeIn)
# Choose grouping variable — if none, can leave it NULL or use 'Size'
group_var <- "Size"
# CLR normalization (best for compositional comparison)
result_clr <- normalization_set(physeq4, method = "clr", groups = group_var)## Removing 214 features with zero counts across all samples.# DESeq normalization
result_deseq <- normalization_set(physeq4, method = "DESeq", groups = group_var)## Removing 214 features with zero counts across all samples.## converting counts to integer mode## using supplied model matrix## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 323 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testing## Removing 214 features with zero counts across all samples.# Visualize for sanity check
library(ggplot2)
boxplot(otu_table(result_clr$dat.normed), main = "CLR")
# Visualizing library size differences.
library(dplyr)
library(tidyr)
library(ggplot2)
library(microbiome)
# Step 1: Use DESeq2-normalized phyloseq object
physeq7 <- result_clr$dat.normed
# Step 2: Extract metadata + add normalized library sizes
meta_raw <- microbiome::meta(physeq7)
meta_raw$LibrarySize <- sample_sums(physeq7)
# Step 3: Create paired dataset (mean per mv/bulk per KFTno)
paired_df <- meta_raw %>%
filter(Size %in% c("mv", "bulk")) %>%
select(KFTno, Size, LibrarySize) %>%
group_by(KFTno, Size) %>%
summarise(LibrarySize = mean(LibrarySize), .groups = "drop") %>%
pivot_wider(names_from = Size, values_from = LibrarySize) %>%
filter(!is.na(mv), !is.na(bulk))
# Step 4: Correlation test
cor_result <- cor.test(paired_df$mv, paired_df$bulk, method = "spearman")
print(cor_result)##
## Spearman's rank correlation rho
##
## data: paired_df$mv and paired_df$bulk
## S = 10, p-value = 0.03413
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8214286# Step 5: Visualization of MV vs Bulk read counts
paired_df_long <- paired_df %>%
pivot_longer(cols = c("mv", "bulk"), names_to = "Extraction", values_to = "LibrarySize")
ggplot(paired_df_long, aes(x = as.factor(KFTno), y = LibrarySize, fill = Extraction)) +
geom_bar(stat = "identity", position = "dodge") +
theme_minimal() +
labs(title = "Normalized Library Size Comparison (MV vs Bulk)",
x = "KFT Number", y = "Normalized Read Count", fill = "Extraction") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +scale_fill_manual(values = DspikeIn::color_palette$mix_MG)
14 Rarefaction
Step 3: Rarefy dataset for fair comparisons visually assess sample richness and sequencing depth sufficiency
otu.rare <- otu_table(physeq4)
otu.rare <- as.data.frame(t(otu.rare))
sample_names <- rownames(otu.rare)
library(vegan)
# Fix: Set cex to a small positive number (like 0.5), or omit it
vegan::rarecurve(otu.rare, step = 100, cex = 0.5, label = TRUE)
15 Choosing rarefaction depth depending on sequencing depth distribution
16 ordination
# Create ordination object if not already done
ord <- ordinate(physeq_rare_4000, method = "PCoA", distance = "bray")
# Ensure consistent factor levels for plotting
physeq_rare_4000@sam_data$KFTno <- factor(
physeq_rare_4000@sam_data$KFTno,
levels = c(83, 85, 92, 100, 101, 102, 114)
)
# Plot
pcoa_plot <- plot_ordination(physeq_rare_4000, ord, type = "samples", color = "KFTno", shape = "Method") +
scale_colour_manual(values = c(
"83" = "#CC79A7",
"85" = "pink",
"92" = "blue",
"100" = "red",
"101" = "#56B4E9",
"102" = "#009E73",
"114" = "#0072B2"
)) +
geom_point(size = 5, alpha = 0.7) +
theme_minimal() +
labs(
title = "PCoA of Bray-Curtis Distances",
color = "KFT Number",
shape = "Method"
) +
theme(
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)
)
pcoa_plot
17 whether DNA input volume impacts microbial community structure??
# Convert variables to factors for plotting
physeq_rare_4000@sam_data$Extracted_volume <- factor(
physeq_rare_4000@sam_data$Extracted_volume,
levels = c(100, 200, 400, 1000, 50000)
)
physeq_rare_4000@sam_data$KFTno <- factor(
physeq_rare_4000@sam_data$KFTno,
levels = c(83, 85, 92, 100, 101, 102, 114)
)
# Run PCoA ordination
ord <- ordinate(physeq_rare_4000, method = "PCoA", distance = "bray")
# Create the plot
pcoa_plot <- plot_ordination(physeq_rare_4000, ord, type = "samples",
color = "KFTno", shape = "Method") +
geom_point(aes(size = Extracted_volume), alpha = 0.7) +
# Method shapes
scale_shape_manual(values = c(
"conventional" = 19, # Circle
"FT" = 15, # Square
"FTP" = 17, # Triangle
"FTPIG" = 8 # Star
)) +
# KFT colors
scale_color_manual(values = c(
"83" = "#CC33CC",
"85" = "#FFA07A",
"92" = "#4682B4",
"100" = "red",
"101" = "#56B4E9",
"102" = "#009E73",
"114" = "#0072B2"
)) +
# Size mapped to Extracted Volume
scale_size_manual(
name = "Extracted Volume (µL)",
values = c("100" = 3, "200" = 4, "400" = 5, "1000" = 6, "50000" = 7),
breaks = c("100", "200", "400", "1000", "50000"),
labels = c("100 µL", "200 µL", "400 µL", "1000 µL", "50,000 µL")
) +
# Confidence ellipses per group
stat_ellipse(aes(group = KFTno, color = KFTno),
type = "t", linetype = "solid", size = 0.8, level = 0.95, alpha = 0.3) +
# Theme and labels
labs(
title = "PCoA of Bray-Curtis Distances",
color = "KFT Sample Number",
shape = "Extraction Method",
size = "Extracted Volume"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "right",
panel.grid = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12)
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
# barplot
library(DspikeIn)
library(ggplot2)
# Ensure sample codes are readable
physeq_rare_4000@sam_data$KFTno <- factor(
physeq_rare_4000@sam_data$KFTno,
levels = c(83, 85, 92, 100, 101, 102, 114)
)
bp_ab <- taxa_barplot(
physeq_rare_4000,
target_glom = "Order",
treatment_variable = "Code",
abundance_type = "relative",
x_angle = 0,
fill_variable = "Order",
facet_variable = "KFTno",
top_n_taxa = 30,
palette = color_palette$mix_MG
)
# Customize theme and improve layout
bp_ab$barplot <- bp_ab$barplot +
facet_grid(~KFTno, scales = "free_x") +
theme_minimal(base_size = 14) +
theme(
legend.position = "right",
legend.box = "vertical",
legend.spacing.x = unit(0.3, "cm"),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 10),
axis.text.y = element_text(size = 12),
strip.text.x = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12)
) +
labs(
title = "Relative Abundance of Top 30 Orders per Sample (Faceted by KFTno)",
y = "Relative Abundance (%)",
x = "Sample Code",
fill = "Order"
)
# Print the adjusted plot
print(bp_ab$barplot)
18 function for betadispers, Tukey HSD and Adonis
run_beta_div_tests <- function(physeq_obj, group_var = "Method", strata_var = "KFTno") {
library(vegan)
# # First, install devtools if you don't have it
# if (!requireNamespace("devtools", quietly = TRUE)) {
# install.packages("devtools")
# }
#
# # Then install pairwiseAdonis correctly
# devtools::install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(pairwiseAdonis)
library(microbiome)
meta_df <- microbiome::meta(physeq_obj)
dist_bray <- phyloseq::distance(physeq_obj, method = "bray", weighted = TRUE)
message("Running betadisper for dispersion homogeneity")
betadispmod <- betadisper(dist_bray, meta_df[[group_var]])
print(anova(betadispmod))
message("Running permutest for dispersion differences")
perm_result <- permutest(betadispmod, pairwise = TRUE, permutations = 9999, p.adjust = "fdr")
print(perm_result)
if (length(unique(meta_df[[group_var]])) > 2) {
message("Running Tukey HSD post-hoc on dispersion")
hsd <- TukeyHSD(betadispmod)
hsd_df <- as.data.frame(hsd$group)
hsd_df$Significance <- cut(hsd_df$`p adj`,
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("***", "**", "*", "ns"))
print(hsd_df)
}
message("Running PERMANOVA (adonis2)")
formula <- as.formula(paste("dist_bray ~", group_var))
adonis_res <- adonis2(formula, data = meta_df, strata = meta_df[[strata_var]], permutations = 9999)
print(adonis_res)
message("Running pairwise PERMANOVA with Holm correction")
pairwise_res <- pairwise.adonis2(formula,
data = meta_df,
strata = meta_df[[strata_var]],
p.adjust = "holm")
print(pairwise_res)
}19 subset
# Set rarefied object
physeq_rare <- physeq_rare_4000
# Subset 1: Early September to Mid-November
subset_2_30Sep_15Nov <- subset_samples(physeq_rare, KFTno %in% c(92, 100, 114))
# Subset 2: Mid-August
subset_5_12Aug <- subset_samples(physeq_rare, KFTno %in% c(83, 85))
# Subset 3: Early October
subset_4_7Oct <- subset_samples(physeq_rare, KFTno %in% c(101, 102))## Loading required package: cluster## Running betadisper for dispersion homogeneity## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.038767 0.0193836 7.3854 0.009257 **
## Residuals 11 0.028871 0.0026246
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running permutest for dispersion differences##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.038767 0.0193836 7.3854 9999 0.01 **
## Residuals 11 0.028871 0.0026246
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## conventional FT FTP
## conventional 0.1022000 0.0048
## FT 0.1024159 0.3372
## FTP 0.0018204 0.3592578## Running Tukey HSD post-hoc on dispersion## diff lwr upr p adj Significance
## FT-conventional -0.08100842 -0.1788487 0.01683189 0.108580473 ns
## FTP-conventional -0.11632047 -0.2001059 -0.03253504 0.008280828 **
## FTP-FT -0.03531205 -0.1363611 0.06573698 0.625304688 ns## Running PERMANOVA (adonis2)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = formula, data = meta_df, permutations = 9999, strata = meta_df[[strata_var]])
## Df SumOfSqs R2 F Pr(>F)
## Model 2 0.27925 0.33024 2.7119 0.0045 **
## Residual 11 0.56636 0.66976
## Total 13 0.84561 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running pairwise PERMANOVA with Holm correction## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## $parent_call
## [1] "dist_bray ~ Method , strata = 1 , permutations 999"
##
## $conventional_vs_FT
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.10425 0.17997 1.5363 1
## Residual 7 0.47502 0.82003
## Total 8 0.57927 1.00000
##
## $conventional_vs_FTP
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.20630 0.30274 3.9077 1
## Residual 9 0.47513 0.69726
## Total 10 0.68143 1.00000
##
## $FT_vs_FTP
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.087752 0.32462 2.8839 1
## Residual 6 0.182572 0.67538
## Total 7 0.270324 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"# Run on September–November subset by Volume
run_beta_div_tests(subset_2_30Sep_15Nov, group_var = "Extracted_volume")## Running betadisper for dispersion homogeneity## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.029158 0.0291578 5.0564 0.03155 *
## Residuals 32 0.184528 0.0057665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running permutest for dispersion differences##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.029158 0.0291578 5.0564 9999 0.031 *
## Residuals 32 0.184528 0.0057665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## 100 50000
## 100 0.0319
## 50000 0.031551## Running PERMANOVA (adonis2)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = formula, data = meta_df, permutations = 9999, strata = meta_df[[strata_var]])
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.5041 0.08761 3.0729 1e-04 ***
## Residual 32 5.2493 0.91239
## Total 33 5.7534 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running pairwise PERMANOVA with Holm correction## $parent_call
## [1] "dist_bray ~ Extracted_volume , strata = 3 , permutations 999"
##
## $`50000_vs_100`
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.5041 0.08761 3.0729 1
## Residual 32 5.2493 0.91239
## Total 33 5.7534 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"## Running betadisper for dispersion homogeneity## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.020563 0.0068545 3.1972 0.04559 *
## Residuals 20 0.042878 0.0021439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running permutest for dispersion differences##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.020563 0.0068545 3.1972 9999 0.0401 *
## Residuals 20 0.042878 0.0021439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## conventional FT FTP FTPIG
## conventional 0.063000 0.219600 0.0134
## FT 0.070695 0.865400 0.2133
## FTP 0.215256 0.854162 0.2469
## FTPIG 0.018077 0.210638 0.239523## Running Tukey HSD post-hoc on dispersion## diff lwr upr p adj Significance
## FT-conventional 0.036298547 -0.037583355 0.11018045 0.52851132 ns
## FTP-conventional 0.032056266 -0.041825637 0.10593817 0.62519658 ns
## FTPIG-conventional 0.077348220 0.007357608 0.14733883 0.02702579 *
## FTP-FT -0.004242282 -0.086206893 0.07772233 0.99887787 ns
## FTPIG-FT 0.041049673 -0.037425468 0.11952481 0.47645712 ns
## FTPIG-FTP 0.045291955 -0.033183187 0.12376710 0.39287912 ns## Running PERMANOVA (adonis2)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = formula, data = meta_df, permutations = 9999, strata = meta_df[[strata_var]])
## Df SumOfSqs R2 F Pr(>F)
## Model 3 0.71508 0.43502 5.1332 1e-04 ***
## Residual 20 0.92870 0.56498
## Total 23 1.64378 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Running pairwise PERMANOVA with Holm correction## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## $parent_call
## [1] "dist_bray ~ Method , strata = 1 , permutations 999"
##
## $conventional_vs_FTP
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.36298 0.48079 10.186 1
## Residual 11 0.39198 0.51921
## Total 12 0.75496 1.00000
##
## $conventional_vs_FTPIG
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.36136 0.39449 7.818 1
## Residual 12 0.55466 0.60551
## Total 13 0.91601 1.00000
##
## $conventional_vs_FT
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.34701 0.46826 9.6869 1
## Residual 11 0.39405 0.53174
## Total 12 0.74105 1.00000
##
## $FTP_vs_FTPIG
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.06910 0.11445 1.1631 1
## Residual 9 0.53466 0.88555
## Total 10 0.60375 1.00000
##
## $FTP_vs_FT
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.07509 0.16718 1.6059 1
## Residual 8 0.37405 0.83282
## Total 9 0.44913 1.00000
##
## $FTPIG_vs_FT
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.12383 0.18747 2.0765 1
## Residual 9 0.53672 0.81253
## Total 10 0.66056 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"20 Beta Dispersion with Tukey HSD seprately & Visualization
# ──────────────────────────────────────────────────────────────
# Beta Dispersion Test (Bray-Curtis)
# ──────────────────────────────────────────────────────────────
bray_dist <- phyloseq::distance(physeq_rare_4000, method = "bray")
metadata <- sample_data(physeq_rare_4000)
# Beta dispersion model
betadisp_mod <- betadisper(bray_dist, group = metadata$Method)
# ANOVA and Tukey HSD
anova(betadisp_mod)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.08195 0.0273178 2.8684 0.04279 *
## Residuals 68 0.64761 0.0095237
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1HSD <- TukeyHSD(betadisp_mod)
HSD_results <- as.data.frame(HSD$group)
HSD_results$Significance <- cut(HSD_results$`p adj`,
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("***", "**", "*", "ns"))
print(HSD_results)## diff lwr upr p adj Significance
## FT-conventional -0.11257913 -0.21750819 -0.007650068 0.03074004 *
## FTP-conventional -0.04935408 -0.12521767 0.026509507 0.32477826 ns
## FTPIG-conventional -0.03497003 -0.11511092 0.045170868 0.66055352 ns
## FTP-FT 0.06322505 -0.04288969 0.169339791 0.40287352 ns
## FTPIG-FT 0.07760910 -0.03160453 0.186822735 0.24988719 ns
## FTPIG-FTP 0.01438405 -0.06730311 0.096071217 0.96667674 ns
21 Visualization of Dispersion
# ──────────────────────────────────────────────
# Visualization of Dispersion
# ──────────────────────────────────────────────
# Extract scores and distances
betscores <- as.data.frame(scores(betadisp_mod, display = "sites"))
betscores$Distance <- betadisp_mod$distances
betscores$Method <- physeq_rare_4000@sam_data$Method
betscores$KFTno <- physeq_rare_4000@sam_data$KFTno
betscores$Extracted_volume <- as.factor(physeq_rare_4000@sam_data$Extracted_volume)
betscores$Method <- factor(betscores$Method, levels = c("conventional", "FT", "FTP", "FTPIG"))
# Plotdispersplot <- ggplot(betscores, aes(x = Method, y = Distance)) +
geom_boxplot(outlier.shape = NA, fill = "white", color = "black") +
geom_jitter(aes(shape = Method, color = KFTno), size = 4) + # fixed
scale_shape_manual(values = c(
"conventional" = 19,
"FT" = 15,
"FTP" = 17,
"FTPIG" = 9
)) +
scale_colour_manual(values = c(
"100" = "red",
"101" = "#56B4E9",
"102" = "#009E73",
"11" = "#F0E442",
"114" = "#0072B2",
"31" = "#D55E00",
"83" = "#CC33CC",
"84" = "#CC79A7",
"85" = "#FFA07A",
"86" = "#FDE725",
"90" = "orange",
"92" = "blue",
"95" = "darkgreen"
)) +
labs(shape = "Extraction Method", color = "Sample (KFTno)") +
theme_bw() +
ylab("Distance to Centroid") +
ggtitle("Beta Dispersion") +
theme(
text = element_text(family = "sans", size = 12), # fixed font
strip.text.x = element_text(face = "bold"),
strip.text.y = element_text(face = "bold"),
strip.background = element_rect(colour = "black", fill = "gray98")
)
dispersplot
22 PERMANOVA
library(phyloseq)
library(vegan)
# Install devtools if not already installed
if (!requireNamespace("devtools", quietly = TRUE)) {
install.packages("devtools")
}
# Load devtools
library(devtools)## Loading required package: usethis##
## Attaching package: 'devtools'## The following object is masked from 'package:emmeans':
##
## test## The following object is masked from 'package:permute':
##
## check# Install the pairwiseAdonis package
install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")## Skipping install of 'pairwiseAdonis' from a github remote, the SHA1 (cb190f76) has not changed since last install.
## Use `force = TRUE` to force installationlibrary(pairwiseAdonis)
library(cluster)
# Extract metadata and distance
dist_bray <- phyloseq::distance(subset_5_12Aug, method = "bray")
metadata <- microbiome::meta(subset_5_12Aug)
# PERMANOVA for Method
adonis_method <- adonis2(dist_bray ~ Method, data = metadata, strata = metadata$KFTno, permutations = 9999)
print(adonis_method)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = dist_bray ~ Method, data = metadata, permutations = 9999, strata = metadata$KFTno)
## Df SumOfSqs R2 F Pr(>F)
## Model 2 0.27925 0.33024 2.7119 0.0046 **
## Residual 11 0.56636 0.66976
## Total 13 0.84561 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Pairwise PERMANOVA for Method
pairwise_method <- pairwise.adonis2(
dist_bray ~ Method,
data = metadata,
strata = metadata$KFTno,
p.adjust = "holm"
)## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## $parent_call
## [1] "dist_bray ~ Method , strata = 1 , permutations 999"
##
## $conventional_vs_FT
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.10425 0.17997 1.5363 1
## Residual 7 0.47502 0.82003
## Total 8 0.57927 1.00000
##
## $conventional_vs_FTP
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.20630 0.30274 3.9077 1
## Residual 9 0.47513 0.69726
## Total 10 0.68143 1.00000
##
## $FT_vs_FTP
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.087752 0.32462 2.8839 1
## Residual 6 0.182572 0.67538
## Total 7 0.270324 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"# PERMANOVA for Volume
adonis_volume <- adonis2(dist_bray ~ Extracted_volume, data = metadata, strata = metadata$KFTno, permutations = 9999)
print(adonis_volume)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = dist_bray ~ Extracted_volume, data = metadata, permutations = 9999, strata = metadata$KFTno)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.19150 0.22646 3.5132 0.0049 **
## Residual 12 0.65411 0.77354
## Total 13 0.84561 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Pairwise PERMANOVA for Volume
pairwise_volume <- pairwise.adonis2(
dist_bray ~ Extracted_volume,
data = metadata,
strata = metadata$KFTno,
p.adjust = "holm"
)## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.## $parent_call
## [1] "dist_bray ~ Extracted_volume , strata = 1 , permutations 999"
##
## $`50000_vs_1000`
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.19150 0.22646 3.5132 1
## Residual 12 0.65411 0.77354
## Total 13 0.84561 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"23 DESeq2 differential abundance
# DESeq2 Analysis Pipeline
# ───────────────────────────────────────────────
# Step 1: Subset your phyloseq object to compare only FTPIG vs conventional
# ───────────────────────────────────────────────
Subset_Dseq <- subset_samples(physeq_rare_4000, Method %in% c("conventional", "FTPIG"))
Subset_Dseq <- prune_samples(sample_sums(Subset_Dseq) > 0, Subset_Dseq) # remove empty
# Optional sanity check
table(sample_data(Subset_Dseq)$Method)##
## conventional FTPIG
## 24 18# ───────────────────────────────────────────────
# Step 2: Convert to DESeq2 object and run DESeq
# ───────────────────────────────────────────────
dds <- phyloseq_to_deseq2(Subset_Dseq, ~ Method)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 80 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testing# ───────────────────────────────────────────────
# Step 3: Variance Stabilizing Transformation (VST)
# ───────────────────────────────────────────────
vsd <- varianceStabilizingTransformation(dds, blind = TRUE)
# Create new phyloseq object from VST data
otu_vst <- otu_table(assay(vsd), taxa_are_rows = TRUE)
ps_vst <- phyloseq(otu_vst, tax_table(Subset_Dseq), sample_data(Subset_Dseq))
# ───────────────────────────────────────────────
# Step 4: Run differential abundance with visualization via DspikeIn
# ───────────────────────────────────────────────
results_DESeq2 <- perform_and_visualize_DA(
obj = ps_vst,
method = "DESeq2",
group_var = "Method",
contrast = c("FTPIG", "conventional"),
output_csv_path = "results_DESeq2_all.csv",
target_glom = "Genus",
significance_level = 0.05
)## converting counts to integer mode## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testing## baseMean logFC lfcSE stat pvalue padj FDR
## 1 1.000000 4.172172e-17 0.4658283 8.956460e-17 1.0000000 1 1
## 2 2.166667 6.454029e-02 0.3214107 2.008032e-01 0.8408525 1 1
## 3 1.000000 4.172172e-17 0.4658283 8.956460e-17 1.0000000 1 1
## 4 1.023810 -5.889367e-02 0.4621417 -1.274364e-01 0.8985950 1 1
## 5 2.023810 -2.974729e-02 0.3324153 -8.948833e-02 0.9286938 1 1
## 6 3.976190 -2.017783e-02 0.2386861 -8.453712e-02 0.9326294 1 1
## Significance group OTU Kingdom
## 1 Not Significant conventional 0022f817a6bf1db59594dc65d5944954 Bacteria
## 2 Not Significant conventional 019c7df12e27631b56f1488d691c9cc6 Bacteria
## 3 Not Significant conventional 021f91de49c1c31e0ef11831bde954a0 Bacteria
## 4 Not Significant FTPIG 031c33b1499e82bceed61a3870fc5412 Bacteria
## 5 Not Significant FTPIG 03f9b4536cfa51fc6adddc65b22bbfc8 Bacteria
## 6 Not Significant FTPIG 047c3e377c58d16788eb0ca27fcc9d8b Bacteria
## Phylum Class Order
## 1 Proteobacteria Gammaproteobacteria Burkholderiales_597441
## 2 Planctomycetota Planctomycetia Pirellulales
## 3 Firmicutes_A Clostridia_258483 Acetivibrionales
## 4 Myxococcota_A_473307 UBA4151 UBA4151
## 5 Bacteroidota Bacteroidia Chitinophagales
## 6 Omnitrophota Koll11 UBA10015
## Family Genus Species
## 1 Gallionellaceae Sideroxydans <NA>
## 2 Pirellulaceae Rhodopirellula <NA>
## 3 DSM-8532 Thermoclostridium_A <NA>
## 4 UBA4151 CACGMS01 <NA>
## 5 Saprospiraceae Portibacter <NA>
## 6 Kpj58rc 2-02-FULL-52-10 <NA>
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1 taxa and 42 samples ]
## sample_data() Sample Data: [ 42 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 1 taxa by 7 taxonomic ranks ]# Save filtered results
write.csv(results_DESeq2$results, "results_DESeq2_significant.csv")
# ───────────────────────────────────────────────
# inspect taxonomy of significant taxa
# ───────────────────────────────────────────────
tax_data <- as.data.frame(tax_table(results_DESeq2$obj_significant))
head(tax_data)## Kingdom Phylum Class
## c497da3b39f30aceede6bec3b03cd100 Bacteria Proteobacteria Alphaproteobacteria
## Order Family
## c497da3b39f30aceede6bec3b03cd100 Pelagibacterales Pelagibacteraceae
## Genus Species
## c497da3b39f30aceede6bec3b03cd100 Pelagibacter_A_533952 <NA>24 Alpha diversity (Shannon)
# ───────────────────────────────────────────────
# Step 1: Estimate Shannon Diversity from rarefied phyloseq object
# ───────────────────────────────────────────────
diversity <- estimate_richness(physeq_rare_4000, measures = "Shannon")
diversity_data <- cbind(diversity, microbiome::meta(physeq_rare_4000))
# ───────────────────────────────────────────────
# Step 2: Test normality of Shannon diversity
# ───────────────────────────────────────────────
shapiro_result <- shapiro.test(diversity_data$Shannon)
print(shapiro_result)##
## Shapiro-Wilk normality test
##
## data: diversity_data$Shannon
## W = 0.98372, p-value = 0.4784# scale for statistical comparison or plot normalization
diversity_data$Shannon_scaled <- scale(diversity_data$Shannon)
# Build the plot
bp_shannon <- ggplot(diversity_data, aes(x = KFTno, y = Shannon, fill = as.factor(KFTno))) +
geom_boxplot(outlier.shape = NA, alpha = 0.7, position = position_dodge(width = 0.8)) +
geom_jitter(
aes(
color = as.factor(KFTno),
shape = Method,
fill = as.factor(KFTno),
size = Extracted_volume
),
stroke = 1.5,
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.8)
) +
facet_grid(~KFTno, scales = "free_x") +
theme_minimal() +
theme(
legend.position = "right",
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 14),
axis.text.y = element_text(size = 14)
) +
labs(
x = "KFT Sample Number",
y = "Shannon Diversity Index",
title = "Shannon Diversity Boxplots with Extracted Volume Mapping"
) +
# Fill colors for KFTno
scale_fill_manual(values = c(
"100" = "red",
"101" = "#56B4E9",
"102" = "#009E73",
"11" = "#F0E442",
"114" = "#0072B2",
"31" = "#D55E00",
"83" = "#CC33CC",
"84" = "#CC79A7",
"85" = "#FFA07A",
"86" = "#FDE725",
"90" = "orange",
"92" = "blue",
"95" = "darkgreen"
)) +
# Outline color (black for all)
scale_colour_manual(values = rep("black", length(unique(diversity_data$KFTno)))) +
# Shapes for Method
scale_shape_manual(values = c(
"conventional" = 21,
"FT" = 22,
"FTP" = 24,
"FTPIG" = 23,
"FTPIG-X10" = 25,
"FTPIG-ch-X10" = 21
)) +
# Size scale for Extracted Volume
scale_size_manual(
name = "Extracted Volume (µL)",
values = c("100" = 3, "200" = 5, "400" = 7, "1000" = 9, "50000" = 12),
breaks = c("100", "200", "400", "1000", "50000"),
labels = c("100 µL", "200 µL", "400 µL", "1000 µL", "50000 µL")
)
# Print the plot
print(bp_shannon)
25 LMM models
Fixed effects: KFTno, Method Random effect: (1 | Month:Temp) (nested environmental covariates)
# Perform pairwise comparisons for Method
# Load necessary libraries
library(lme4)
library(lmerTest)
library(emmeans)
# Ensure categorical variables are factors
diversity_data$KFTno <- as.factor(diversity_data$KFTno)
diversity_data$Method <- as.factor(diversity_data$Method)
diversity_data$Month <- as.factor(diversity_data$Month)
diversity_data$Temp <- as.factor(diversity_data$Temp)
diversity_data$Purification <- as.factor(diversity_data$Purification)
diversity_data$Extracted_volume <- as.factor(diversity_data$Extracted_volume)
# ================================
# Model definitions
# ================================
# Null model: Random effect of KFTno only
lmm4 <- lmer(Shannon ~ 1 + (1 | KFTno), data = diversity_data)
# Model with Method, KFTno, and Month:Temp interaction as random effect
lmm6 <- lmer(
Shannon ~ KFTno + Method + (1 | Month:Temp),
data = diversity_data
)## boundary (singular) fit: see help('isSingular')# Model with extra random effect for purification
lmm5 <- lmer(
Shannon ~ KFTno + Method + (1 | Month:Temp) + (1 | Purification),
data = diversity_data
)## boundary (singular) fit: see help('isSingular')## Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
## eigenvalue close to zero: 1.4e-09# Model with Extracted Volume + Month:Temp
lmm3 <- lmer(
Shannon ~ KFTno + Method + (1 | Extracted_volume) + (1 | Month:Temp),
data = diversity_data
)## boundary (singular) fit: see help('isSingular')# ================================
# Model comparison
# ================================
AIC(lmm3, lmm4, lmm5, lmm6)## df AIC
## lmm3 13 -11.35319
## lmm4 3 -15.12762
## lmm5 13 -10.14720
## lmm6 12 -12.14720## df BIC
## lmm3 13 18.243466
## lmm4 3 -8.297622
## lmm5 13 19.449456
## lmm6 12 15.172789# ================================
# Summary & ANOVA of best model (e.g., lmm6)
# ================================
summary(lmm6)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ KFTno + Method + (1 | Month:Temp)
## Data: diversity_data
##
## REML criterion at convergence: -36.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8325 -0.3925 0.1126 0.6218 1.9190
##
## Random effects:
## Groups Name Variance Std.Dev.
## Month:Temp (Intercept) 0.0000 0.0000
## Residual 0.0226 0.1503
## Number of obs: 72, groups: Month:Temp, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.195862 0.064810 62.000000 64.741 < 2e-16 ***
## KFTno85 -0.002113 0.081429 62.000000 -0.026 0.97939
## KFTno92 0.108286 0.078521 62.000000 1.379 0.17283
## KFTno100 0.238932 0.078521 62.000000 3.043 0.00343 **
## KFTno101 0.173290 0.075071 62.000000 2.308 0.02433 *
## KFTno102 0.216348 0.077244 62.000000 2.801 0.00679 **
## KFTno114 0.484704 0.077057 62.000000 6.290 3.58e-08 ***
## MethodFT 0.265639 0.063704 62.000000 4.170 9.64e-05 ***
## MethodFTP 0.247171 0.044768 62.000000 5.521 7.05e-07 ***
## MethodFTPIG 0.133806 0.048660 62.000000 2.750 0.00780 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) KFTn85 KFTn92 KFT100 KFT101 KFT102 KFT114 MthdFT MthFTP
## KFTno85 -0.692
## KFTno92 -0.734 0.586
## KFTno100 -0.734 0.586 0.667
## KFTno101 -0.757 0.622 0.660 0.660
## KFTno102 -0.749 0.601 0.647 0.647 0.672
## KFTno114 -0.756 0.596 0.674 0.674 0.670 0.657
## MethodFT -0.235 -0.072 0.072 0.072 -0.074 -0.030 0.083
## MethodFTP -0.281 -0.043 -0.078 -0.078 -0.023 0.006 -0.053 0.308
## MethodFTPIG -0.147 -0.027 -0.198 -0.198 -0.136 -0.135 -0.175 0.259 0.452
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## KFTno 1.55944 0.25991 6 62 11.498 1.318e-08 ***
## Method 0.84173 0.28058 3 62 12.412 1.856e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ================================
# Pairwise comparisons
# ================================
# 1. Method overall
pairwise_method <- emmeans(lmm6, pairwise ~ Method, adjust = "holm")
summary(pairwise_method)## $emmeans
## Method emmean SE df lower.CL upper.CL
## conventional 4.37 0.0334 0.74 3.36 5.38
## FT 4.64 0.0563 8.15 4.51 4.77
## FTP 4.62 0.0377 0.84 3.86 5.38
## FTPIG 4.50 0.0446 1.49 4.23 4.78
##
## Results are averaged over the levels of: KFTno
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## conventional - FT -0.2656 0.0662 62.0 -4.011 0.0008
## conventional - FTP -0.2472 0.0551 34.9 -4.483 0.0005
## conventional - FTPIG -0.1338 0.0606 32.7 -2.209 0.1030
## FT - FTP 0.0185 0.0671 61.9 0.275 0.7841
## FT - FTPIG 0.1318 0.0717 62.0 1.838 0.1417
## FTP - FTPIG 0.1134 0.0492 61.2 2.306 0.0980
##
## Results are averaged over the levels of: KFTno
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 6 tests# 2. Method within each KFTno (interaction)
pairwise_KFT_method <- emmeans(lmm6, pairwise ~ KFTno:Method, adjust = "holm")
summary(pairwise_KFT_method)## $emmeans
## KFTno Method emmean SE df lower.CL upper.CL
## 83 conventional 4.20 0.0663 0.46 -1.24e+01 2.08e+01
## 85 conventional 4.19 0.0616 0.32 -1.98e+02 2.07e+02
## 92 conventional 4.30 0.0597 0.81 2.90e+00 5.70e+00
## 100 conventional 4.43 0.0592 0.21 -2.00e+04 2.00e+04
## 101 conventional 4.37 0.0536 0.16 -2.11e+06 2.11e+06
## 102 conventional 4.41 0.0550 0.18 -1.52e+05 1.52e+05
## 114 conventional 4.68 0.0557 0.17 -5.02e+05 5.02e+05
## 83 FT 4.46 0.0796 1.07 3.60e+00 5.30e+00
## 85 FT 4.46 0.0700 0.64 9.00e-01 8.10e+00
## 92 FT 4.57 0.0875 2.61 4.30e+00 4.90e+00
## 100 FT 4.70 0.0760 0.87 3.30e+00 6.10e+00
## 101 FT 4.63 0.0623 0.40 -3.43e+01 4.35e+01
## 102 FT 4.68 0.0670 0.54 -2.80e+00 1.22e+01
## 114 FT 4.95 0.0745 0.81 3.30e+00 6.60e+00
## 83 FTP 4.44 0.0701 0.54 -3.00e+00 1.19e+01
## 85 FTP 4.44 0.0619 0.34 -1.43e+02 1.52e+02
## 92 FTP 4.55 0.0780 0.49 -9.70e+00 1.88e+01
## 100 FTP 4.68 0.0527 0.20 -5.26e+04 5.27e+04
## 101 FTP 4.62 0.0533 0.19 -8.56e+04 8.57e+04
## 102 FTP 4.66 0.0571 0.25 -2.49e+03 2.50e+03
## 114 FTP 4.93 0.0518 0.18 -2.30e+05 2.30e+05
## 83 FTPIG 4.33 0.0784 0.82 2.60e+00 6.00e+00
## 85 FTPIG 4.33 0.0716 0.59 -1.10e+00 9.80e+00
## 92 FTPIG 4.44 0.0816 0.50 -9.50e+00 1.83e+01
## 100 FTPIG 4.57 0.0543 0.21 -2.13e+04 2.13e+04
## 101 FTPIG 4.50 0.0564 0.22 -1.02e+04 1.02e+04
## 102 FTPIG 4.55 0.0582 0.25 -2.68e+03 2.69e+03
## 114 FTPIG 4.81 0.0538 0.19 -6.42e+04 6.42e+04
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## KFTno83 conventional - KFTno85 conventional 0.00211 0.0815 0.29 0.026
## KFTno83 conventional - KFTno92 conventional -0.10829 0.0879 0.42 -1.231
## KFTno83 conventional - KFTno100 conventional -0.23893 0.0793 0.25 -3.013
## KFTno83 conventional - KFTno101 conventional -0.17329 0.0753 0.21 -2.301
## KFTno83 conventional - KFTno102 conventional -0.21635 0.0774 0.24 -2.794
## KFTno83 conventional - KFTno114 conventional -0.48470 0.0776 0.23 -6.248
## KFTno83 conventional - KFTno83 FT -0.26564 0.0662 61.95 -4.011
## KFTno83 conventional - KFTno85 FT -0.26353 0.1010 0.65 -2.613
## KFTno83 conventional - KFTno92 FT -0.37392 0.1190 1.20 -3.131
## KFTno83 conventional - KFTno100 FT -0.50457 0.1050 0.80 -4.813
## KFTno83 conventional - KFTno101 FT -0.43893 0.0955 0.54 -4.595
## KFTno83 conventional - KFTno102 FT -0.48199 0.0995 0.63 -4.845
## KFTno83 conventional - KFTno114 FT -0.75034 0.1040 0.78 -7.188
## KFTno83 conventional - KFTno83 FTP -0.24717 0.0551 34.92 -4.483
## KFTno83 conventional - KFTno85 FTP -0.24506 0.0958 0.45 -2.557
## KFTno83 conventional - KFTno92 FTP -0.35546 0.1130 0.53 -3.145
## KFTno83 conventional - KFTno100 FTP -0.48610 0.0898 0.38 -5.412
## KFTno83 conventional - KFTno101 FTP -0.42046 0.0903 0.36 -4.655
## KFTno83 conventional - KFTno102 FTP -0.46352 0.0935 0.42 -4.958
## KFTno83 conventional - KFTno114 FTP -0.73188 0.0901 0.37 -8.125
## KFTno83 conventional - KFTno83 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno83 conventional - KFTno85 FTPIG -0.13169 0.0994 0.49 -1.324
## KFTno83 conventional - KFTno92 FTPIG -0.24209 0.1130 0.44 -2.144
## KFTno83 conventional - KFTno100 FTPIG -0.37274 0.0874 0.32 -4.263
## KFTno83 conventional - KFTno101 FTPIG -0.30710 0.0889 0.32 -3.454
## KFTno83 conventional - KFTno102 FTPIG -0.35015 0.0909 0.34 -3.851
## KFTno83 conventional - KFTno114 FTPIG -0.61851 0.0879 0.32 -7.036
## KFTno85 conventional - KFTno92 conventional -0.11040 0.0843 0.30 -1.310
## KFTno85 conventional - KFTno100 conventional -0.24104 0.0733 0.19 -3.290
## KFTno85 conventional - KFTno101 conventional -0.17540 0.0684 0.14 -2.565
## KFTno85 conventional - KFTno102 conventional -0.21846 0.0710 0.17 -3.076
## KFTno85 conventional - KFTno114 conventional -0.48682 0.0716 0.17 -6.800
## KFTno85 conventional - KFTno83 FT -0.26775 0.1090 0.86 -2.457
## KFTno85 conventional - KFTno85 FT -0.26564 0.0662 61.95 -4.011
## KFTno85 conventional - KFTno92 FT -0.37604 0.1200 1.15 -3.124
## KFTno85 conventional - KFTno100 FT -0.50668 0.1040 0.78 -4.849
## KFTno85 conventional - KFTno101 FT -0.44104 0.0948 0.51 -4.654
## KFTno85 conventional - KFTno102 FT -0.48410 0.0990 0.60 -4.891
## KFTno85 conventional - KFTno114 FT -0.75246 0.1040 0.76 -7.223
## KFTno85 conventional - KFTno83 FTP -0.24928 0.1010 0.51 -2.471
## KFTno85 conventional - KFTno85 FTP -0.24717 0.0551 34.92 -4.483
## KFTno85 conventional - KFTno92 FTP -0.35757 0.1120 0.45 -3.181
## KFTno85 conventional - KFTno100 FTP -0.48822 0.0874 0.33 -5.584
## KFTno85 conventional - KFTno101 FTP -0.42257 0.0875 0.30 -4.829
## KFTno85 conventional - KFTno102 FTP -0.46563 0.0910 0.35 -5.115
## KFTno85 conventional - KFTno114 FTP -0.73399 0.0878 0.32 -8.356
## KFTno85 conventional - KFTno83 FTPIG -0.13592 0.1040 0.53 -1.312
## KFTno85 conventional - KFTno85 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno85 conventional - KFTno92 FTPIG -0.24420 0.1120 0.35 -2.181
## KFTno85 conventional - KFTno100 FTPIG -0.37485 0.0845 0.26 -4.435
## KFTno85 conventional - KFTno101 FTPIG -0.30921 0.0856 0.25 -3.612
## KFTno85 conventional - KFTno102 FTPIG -0.35227 0.0880 0.28 -4.004
## KFTno85 conventional - KFTno114 FTPIG -0.62062 0.0852 0.26 -7.287
## KFTno92 conventional - KFTno100 conventional -0.13065 0.0817 0.17 -1.599
## KFTno92 conventional - KFTno101 conventional -0.06500 0.0782 0.17 -0.831
## KFTno92 conventional - KFTno102 conventional -0.10806 0.0793 0.19 -1.363
## KFTno92 conventional - KFTno114 conventional -0.37642 0.0794 0.16 -4.739
## KFTno92 conventional - KFTno83 FT -0.15735 0.0999 1.05 -1.576
## KFTno92 conventional - KFTno85 FT -0.15524 0.0921 0.72 -1.685
## KFTno92 conventional - KFTno92 FT -0.26564 0.0662 61.95 -4.011
## KFTno92 conventional - KFTno100 FT -0.39629 0.0961 0.75 -4.125
## KFTno92 conventional - KFTno101 FT -0.33064 0.0861 0.50 -3.839
## KFTno92 conventional - KFTno102 FT -0.37370 0.0897 0.61 -4.167
## KFTno92 conventional - KFTno114 FT -0.64206 0.0951 0.74 -6.754
## KFTno92 conventional - KFTno83 FTP -0.13889 0.0937 0.90 -1.483
## KFTno92 conventional - KFTno85 FTP -0.13677 0.0874 0.67 -1.565
## KFTno92 conventional - KFTno92 FTP -0.24717 0.0551 34.92 -4.483
## KFTno92 conventional - KFTno100 FTP -0.37782 0.0804 0.44 -4.702
## KFTno92 conventional - KFTno101 FTP -0.31218 0.0812 0.48 -3.843
## KFTno92 conventional - KFTno102 FTP -0.35523 0.0839 0.56 -4.236
## KFTno92 conventional - KFTno114 FTP -0.62359 0.0800 0.44 -7.795
## KFTno92 conventional - KFTno83 FTPIG -0.02552 0.1000 1.20 -0.254
## KFTno92 conventional - KFTno85 FTPIG -0.02341 0.0949 0.95 -0.247
## KFTno92 conventional - KFTno92 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno92 conventional - KFTno100 FTPIG -0.26445 0.0818 0.49 -3.233
## KFTno92 conventional - KFTno101 FTPIG -0.19881 0.0837 0.56 -2.376
## KFTno92 conventional - KFTno102 FTPIG -0.24187 0.0850 0.60 -2.847
## KFTno92 conventional - KFTno114 FTPIG -0.51022 0.0817 0.50 -6.246
## KFTno100 conventional - KFTno101 conventional 0.06564 0.0636 0.11 1.032
## KFTno100 conventional - KFTno102 conventional 0.02258 0.0657 0.12 0.344
## KFTno100 conventional - KFTno114 conventional -0.24577 0.0628 0.10 -3.911
## KFTno100 conventional - KFTno83 FT -0.02671 0.1020 0.58 -0.262
## KFTno100 conventional - KFTno85 FT -0.02459 0.0927 0.40 -0.265
## KFTno100 conventional - KFTno92 FT -0.13499 0.1140 0.62 -1.189
## KFTno100 conventional - KFTno100 FT -0.26564 0.0662 61.95 -4.011
## KFTno100 conventional - KFTno101 FT -0.20000 0.0847 0.29 -2.361
## KFTno100 conventional - KFTno102 FT -0.24305 0.0889 0.35 -2.735
## KFTno100 conventional - KFTno114 FT -0.51141 0.0923 0.43 -5.544
## KFTno100 conventional - KFTno83 FTP -0.00824 0.1030 0.47 -0.080
## KFTno100 conventional - KFTno85 FTP -0.00613 0.0958 0.35 -0.064
## KFTno100 conventional - KFTno92 FTP -0.11652 0.1140 0.32 -1.023
## KFTno100 conventional - KFTno100 FTP -0.24717 0.0551 34.92 -4.483
## KFTno100 conventional - KFTno101 FTP -0.18153 0.0882 0.26 -2.057
## KFTno100 conventional - KFTno102 FTP -0.22459 0.0912 0.30 -2.463
## KFTno100 conventional - KFTno114 FTP -0.49294 0.0855 0.24 -5.768
## KFTno100 conventional - KFTno83 FTPIG 0.10513 0.1110 0.62 0.949
## KFTno100 conventional - KFTno85 FTPIG 0.10724 0.1050 0.50 1.026
## KFTno100 conventional - KFTno92 FTPIG -0.00316 0.1180 0.34 -0.027
## KFTno100 conventional - KFTno100 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno100 conventional - KFTno101 FTPIG -0.06816 0.0927 0.30 -0.736
## KFTno100 conventional - KFTno102 FTPIG -0.11122 0.0943 0.32 -1.179
## KFTno100 conventional - KFTno114 FTPIG -0.37958 0.0893 0.27 -4.251
## KFTno101 conventional - KFTno102 conventional -0.04306 0.0617 61.01 -0.698
## KFTno101 conventional - KFTno114 conventional -0.31141 0.0619 0.10 -5.029
## KFTno101 conventional - KFTno83 FT -0.09235 0.1050 0.71 -0.881
## KFTno101 conventional - KFTno85 FT -0.09024 0.0956 0.49 -0.944
## KFTno101 conventional - KFTno92 FT -0.20063 0.1170 0.87 -1.722
## KFTno101 conventional - KFTno100 FT -0.33128 0.0984 0.60 -3.368
## KFTno101 conventional - KFTno101 FT -0.26564 0.0662 61.95 -4.011
## KFTno101 conventional - KFTno102 FT -0.30870 0.0930 61.89 -3.320
## KFTno101 conventional - KFTno114 FT -0.57705 0.0982 0.59 -5.876
## KFTno101 conventional - KFTno83 FTP -0.07388 0.0963 0.38 -0.767
## KFTno101 conventional - KFTno85 FTP -0.07177 0.0882 0.27 -0.814
## KFTno101 conventional - KFTno92 FTP -0.18217 0.1080 0.28 -1.683
## KFTno101 conventional - KFTno100 FTP -0.31281 0.0799 0.20 -3.917
## KFTno101 conventional - KFTno101 FTP -0.24717 0.0551 34.92 -4.483
## KFTno101 conventional - KFTno102 FTP -0.29023 0.0843 58.93 -3.442
## KFTno101 conventional - KFTno114 FTP -0.55859 0.0805 0.20 -6.938
## KFTno101 conventional - KFTno83 FTPIG 0.03948 0.1040 0.50 0.380
## KFTno101 conventional - KFTno85 FTPIG 0.04160 0.0968 0.38 0.430
## KFTno101 conventional - KFTno92 FTPIG -0.06880 0.1120 0.29 -0.614
## KFTno101 conventional - KFTno100 FTPIG -0.19945 0.0827 0.22 -2.412
## KFTno101 conventional - KFTno101 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno101 conventional - KFTno102 FTPIG -0.17686 0.0867 56.65 -2.039
## KFTno101 conventional - KFTno114 FTPIG -0.44522 0.0836 0.22 -5.328
## KFTno102 conventional - KFTno114 conventional -0.26836 0.0640 0.11 -4.194
## KFTno102 conventional - KFTno83 FT -0.04929 0.1040 0.70 -0.473
## KFTno102 conventional - KFTno85 FT -0.04718 0.0952 0.49 -0.496
## KFTno102 conventional - KFTno92 FT -0.15758 0.1150 0.85 -1.366
## KFTno102 conventional - KFTno100 FT -0.28822 0.0975 0.58 -2.958
## KFTno102 conventional - KFTno101 FT -0.22258 0.0880 61.84 -2.529
## KFTno102 conventional - KFTno102 FT -0.26564 0.0662 61.95 -4.011
## KFTno102 conventional - KFTno114 FT -0.53400 0.0972 0.57 -5.492
## KFTno102 conventional - KFTno83 FTP -0.03082 0.0966 0.40 -0.319
## KFTno102 conventional - KFTno85 FTP -0.02871 0.0888 0.28 -0.323
## KFTno102 conventional - KFTno92 FTP -0.13911 0.1080 0.29 -1.290
## KFTno102 conventional - KFTno100 FTP -0.26976 0.0799 0.21 -3.375
## KFTno102 conventional - KFTno101 FTP -0.20411 0.0812 59.35 -2.514
## KFTno102 conventional - KFTno102 FTP -0.24717 0.0551 34.92 -4.483
## KFTno102 conventional - KFTno114 FTP -0.51553 0.0805 0.21 -6.403
## KFTno102 conventional - KFTno83 FTPIG 0.08254 0.1050 0.55 0.785
## KFTno102 conventional - KFTno85 FTPIG 0.08465 0.0984 0.42 0.860
## KFTno102 conventional - KFTno92 FTPIG -0.02574 0.1130 0.32 -0.228
## KFTno102 conventional - KFTno100 FTPIG -0.15639 0.0840 0.24 -1.861
## KFTno102 conventional - KFTno101 FTPIG -0.09075 0.0862 58.05 -1.052
## KFTno102 conventional - KFTno102 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno102 conventional - KFTno114 FTPIG -0.40216 0.0848 0.24 -4.740
## KFTno114 conventional - KFTno83 FT 0.21907 0.0996 0.55 2.200
## KFTno114 conventional - KFTno85 FT 0.22118 0.0904 0.37 2.448
## KFTno114 conventional - KFTno92 FT 0.11078 0.1110 0.59 0.997
## KFTno114 conventional - KFTno100 FT -0.01987 0.0903 0.41 -0.220
## KFTno114 conventional - KFTno101 FT 0.04577 0.0824 0.26 0.555
## KFTno114 conventional - KFTno102 FT 0.00272 0.0866 0.32 0.031
## KFTno114 conventional - KFTno114 FT -0.26564 0.0662 61.95 -4.011
## KFTno114 conventional - KFTno83 FTP 0.23753 0.1000 0.43 2.375
## KFTno114 conventional - KFTno85 FTP 0.23965 0.0928 0.32 2.582
## KFTno114 conventional - KFTno92 FTP 0.12925 0.1110 0.29 1.166
## KFTno114 conventional - KFTno100 FTP -0.00140 0.0817 0.21 -0.017
## KFTno114 conventional - KFTno101 FTP 0.06424 0.0852 0.23 0.754
## KFTno114 conventional - KFTno102 FTP 0.02119 0.0882 0.27 0.240
## KFTno114 conventional - KFTno114 FTP -0.24717 0.0551 34.92 -4.483
## KFTno114 conventional - KFTno83 FTPIG 0.35090 0.1080 0.57 3.251
## KFTno114 conventional - KFTno85 FTPIG 0.35301 0.1020 0.45 3.472
## KFTno114 conventional - KFTno92 FTPIG 0.24261 0.1150 0.31 2.105
## KFTno114 conventional - KFTno100 FTPIG 0.11197 0.0852 0.24 1.313
## KFTno114 conventional - KFTno101 FTPIG 0.17761 0.0896 0.27 1.983
## KFTno114 conventional - KFTno102 FTPIG 0.13455 0.0913 0.28 1.474
## KFTno114 conventional - KFTno114 FTPIG -0.13381 0.0606 32.72 -2.209
## KFTno83 FT - KFTno85 FT 0.00211 0.0815 0.29 0.026
## KFTno83 FT - KFTno92 FT -0.10829 0.0879 0.42 -1.231
## KFTno83 FT - KFTno100 FT -0.23893 0.0793 0.25 -3.013
## KFTno83 FT - KFTno101 FT -0.17329 0.0753 0.21 -2.301
## KFTno83 FT - KFTno102 FT -0.21635 0.0774 0.24 -2.794
## KFTno83 FT - KFTno114 FT -0.48470 0.0776 0.23 -6.248
## KFTno83 FT - KFTno83 FTP 0.01847 0.0671 61.91 0.275
## KFTno83 FT - KFTno85 FTP 0.02058 0.1070 0.86 0.192
## KFTno83 FT - KFTno92 FTP -0.08982 0.1100 0.87 -0.817
## KFTno83 FT - KFTno100 FTP -0.22046 0.0960 0.57 -2.296
## KFTno83 FT - KFTno101 FTP -0.15482 0.1030 0.73 -1.507
## KFTno83 FT - KFTno102 FTP -0.19788 0.1030 0.75 -1.914
## KFTno83 FT - KFTno114 FTP -0.46624 0.0953 0.55 -4.890
## KFTno83 FT - KFTno83 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno83 FT - KFTno85 FTPIG 0.13395 0.1110 0.95 1.212
## KFTno83 FT - KFTno92 FTPIG 0.02355 0.1100 0.77 0.214
## KFTno83 FT - KFTno100 FTPIG -0.10710 0.0938 0.51 -1.141
## KFTno83 FT - KFTno101 FTPIG -0.04146 0.1020 0.69 -0.408
## KFTno83 FT - KFTno102 FTPIG -0.08451 0.1010 0.67 -0.836
## KFTno83 FT - KFTno114 FTPIG -0.35287 0.0933 0.50 -3.781
## KFTno85 FT - KFTno92 FT -0.11040 0.0843 0.30 -1.310
## KFTno85 FT - KFTno100 FT -0.24104 0.0733 0.19 -3.290
## KFTno85 FT - KFTno101 FT -0.17540 0.0684 0.14 -2.565
## KFTno85 FT - KFTno102 FT -0.21846 0.0710 0.17 -3.076
## KFTno85 FT - KFTno114 FT -0.48682 0.0716 0.17 -6.800
## KFTno85 FT - KFTno83 FTP 0.01636 0.1040 0.73 0.157
## KFTno85 FT - KFTno85 FTP 0.01847 0.0671 61.91 0.275
## KFTno85 FT - KFTno92 FTP -0.09193 0.1050 0.63 -0.872
## KFTno85 FT - KFTno100 FTP -0.22258 0.0891 0.42 -2.497
## KFTno85 FT - KFTno101 FTP -0.15693 0.0960 0.55 -1.635
## KFTno85 FT - KFTno102 FTP -0.19999 0.0969 0.57 -2.064
## KFTno85 FT - KFTno114 FTP -0.46835 0.0885 0.40 -5.290
## KFTno85 FT - KFTno83 FTPIG 0.12972 0.1070 0.78 1.217
## KFTno85 FT - KFTno85 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno85 FT - KFTno92 FTPIG 0.02144 0.1050 0.53 0.204
## KFTno85 FT - KFTno100 FTPIG -0.10921 0.0863 0.36 -1.265
## KFTno85 FT - KFTno101 FTPIG -0.04357 0.0943 0.50 -0.462
## KFTno85 FT - KFTno102 FTPIG -0.08663 0.0941 0.49 -0.921
## KFTno85 FT - KFTno114 FTPIG -0.35498 0.0859 0.35 -4.131
## KFTno92 FT - KFTno100 FT -0.13065 0.0817 0.17 -1.599
## KFTno92 FT - KFTno101 FT -0.06500 0.0782 0.17 -0.831
## KFTno92 FT - KFTno102 FT -0.10806 0.0793 0.19 -1.363
## KFTno92 FT - KFTno114 FT -0.37642 0.0794 0.16 -4.739
## KFTno92 FT - KFTno83 FTP 0.12675 0.1110 1.58 1.139
## KFTno92 FT - KFTno85 FTP 0.12887 0.1100 1.45 1.171
## KFTno92 FT - KFTno92 FTP 0.01847 0.0671 61.91 0.275
## KFTno92 FT - KFTno100 FTP -0.11218 0.0988 0.78 -1.136
## KFTno92 FT - KFTno101 FTP -0.04654 0.1060 1.16 -0.441
## KFTno92 FT - KFTno102 FTP -0.08959 0.1060 1.17 -0.849
## KFTno92 FT - KFTno114 FTP -0.35795 0.0976 0.77 -3.668
## KFTno92 FT - KFTno83 FTPIG 0.24012 0.1170 2.02 2.053
## KFTno92 FT - KFTno85 FTPIG 0.24223 0.1160 1.91 2.086
## KFTno92 FT - KFTno92 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno92 FT - KFTno100 FTPIG 0.00119 0.1000 0.90 0.012
## KFTno92 FT - KFTno101 FTPIG 0.06683 0.1080 1.33 0.622
## KFTno92 FT - KFTno102 FTPIG 0.02377 0.1060 1.29 0.223
## KFTno92 FT - KFTno114 FTPIG -0.24459 0.0990 0.89 -2.470
## KFTno100 FT - KFTno101 FT 0.06564 0.0636 0.11 1.032
## KFTno100 FT - KFTno102 FT 0.02258 0.0657 0.12 0.344
## KFTno100 FT - KFTno114 FT -0.24577 0.0628 0.10 -3.911
## KFTno100 FT - KFTno83 FTP 0.25740 0.1110 0.90 2.314
## KFTno100 FT - KFTno85 FTP 0.25951 0.1090 0.84 2.389
## KFTno100 FT - KFTno92 FTP 0.14911 0.1120 0.68 1.328
## KFTno100 FT - KFTno100 FTP 0.01847 0.0671 61.91 0.275
## KFTno100 FT - KFTno101 FTP 0.08411 0.1020 0.68 0.821
## KFTno100 FT - KFTno102 FTP 0.04105 0.1030 0.68 0.399
## KFTno100 FT - KFTno114 FTP -0.22730 0.0927 0.46 -2.452
## KFTno100 FT - KFTno83 FTPIG 0.37077 0.1190 1.13 3.127
## KFTno100 FT - KFTno85 FTPIG 0.37288 0.1160 1.08 3.201
## KFTno100 FT - KFTno92 FTPIG 0.26248 0.1170 0.75 2.247
## KFTno100 FT - KFTno100 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno100 FT - KFTno101 FTPIG 0.19748 0.1060 0.76 1.857
## KFTno100 FT - KFTno102 FTPIG 0.15442 0.1060 0.73 1.461
## KFTno100 FT - KFTno114 FTPIG -0.11394 0.0963 0.52 -1.183
## KFTno101 FT - KFTno102 FT -0.04306 0.0617 61.01 -0.698
## KFTno101 FT - KFTno114 FT -0.31141 0.0619 0.10 -5.029
## KFTno101 FT - KFTno83 FTP 0.19176 0.0990 0.58 1.937
## KFTno101 FT - KFTno85 FTP 0.19387 0.0957 0.52 2.026
## KFTno101 FT - KFTno92 FTP 0.08347 0.1000 0.42 0.831
## KFTno101 FT - KFTno100 FTP -0.04717 0.0812 0.28 -0.581
## KFTno101 FT - KFTno101 FTP 0.01847 0.0671 61.91 0.275
## KFTno101 FT - KFTno102 FTP -0.02459 0.0901 61.68 -0.273
## KFTno101 FT - KFTno114 FTP -0.29295 0.0808 0.27 -3.627
## KFTno101 FT - KFTno83 FTPIG 0.30512 0.1060 0.75 2.867
## KFTno101 FT - KFTno85 FTPIG 0.30724 0.1040 0.70 2.962
## KFTno101 FT - KFTno92 FTPIG 0.19684 0.1050 0.46 1.880
## KFTno101 FT - KFTno100 FTPIG 0.06619 0.0840 0.31 0.788
## KFTno101 FT - KFTno101 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno101 FT - KFTno102 FTPIG 0.08878 0.0925 61.89 0.960
## KFTno101 FT - KFTno114 FTPIG -0.17958 0.0839 0.31 -2.142
## KFTno102 FT - KFTno114 FT -0.26836 0.0640 0.11 -4.194
## KFTno102 FT - KFTno83 FTP 0.23482 0.1020 0.65 2.313
## KFTno102 FT - KFTno85 FTP 0.23693 0.0986 0.59 2.403
## KFTno102 FT - KFTno92 FTP 0.12653 0.1020 0.48 1.237
## KFTno102 FT - KFTno100 FTP -0.00412 0.0840 0.32 -0.049
## KFTno102 FT - KFTno101 FTP 0.06153 0.0922 61.53 0.667
## KFTno102 FT - KFTno102 FTP 0.01847 0.0671 61.91 0.275
## KFTno102 FT - KFTno114 FTP -0.24989 0.0835 0.31 -2.992
## KFTno102 FT - KFTno83 FTPIG 0.34818 0.1100 0.87 3.172
## KFTno102 FT - KFTno85 FTPIG 0.35029 0.1070 0.81 3.262
## KFTno102 FT - KFTno92 FTPIG 0.23990 0.1070 0.56 2.233
## KFTno102 FT - KFTno100 FTPIG 0.10925 0.0880 0.38 1.242
## KFTno102 FT - KFTno101 FTPIG 0.17489 0.0967 61.74 1.808
## KFTno102 FT - KFTno102 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno102 FT - KFTno114 FTPIG -0.13652 0.0877 0.37 -1.556
## KFTno114 FT - KFTno83 FTP 0.50317 0.1090 0.85 4.601
## KFTno114 FT - KFTno85 FTP 0.50528 0.1070 0.80 4.728
## KFTno114 FT - KFTno92 FTP 0.39489 0.1100 0.65 3.589
## KFTno114 FT - KFTno100 FTP 0.26424 0.0912 0.44 2.897
## KFTno114 FT - KFTno101 FTP 0.32988 0.1010 0.64 3.273
## KFTno114 FT - KFTno102 FTP 0.28682 0.1010 0.65 2.836
## KFTno114 FT - KFTno114 FTP 0.01847 0.0671 61.91 0.275
## KFTno114 FT - KFTno83 FTPIG 0.61654 0.1170 1.08 5.283
## KFTno114 FT - KFTno85 FTPIG 0.61865 0.1150 1.03 5.394
## KFTno114 FT - KFTno92 FTPIG 0.50825 0.1140 0.71 4.440
## KFTno114 FT - KFTno100 FTPIG 0.37761 0.0944 0.50 3.999
## KFTno114 FT - KFTno101 FTPIG 0.44325 0.1050 0.72 4.240
## KFTno114 FT - KFTno102 FTPIG 0.40019 0.1040 0.70 3.854
## KFTno114 FT - KFTno114 FTPIG 0.13183 0.0717 62.00 1.838
## KFTno83 FTP - KFTno85 FTP 0.00211 0.0815 0.29 0.026
## KFTno83 FTP - KFTno92 FTP -0.10829 0.0879 0.42 -1.231
## KFTno83 FTP - KFTno100 FTP -0.23893 0.0793 0.25 -3.013
## KFTno83 FTP - KFTno101 FTP -0.17329 0.0753 0.21 -2.301
## KFTno83 FTP - KFTno102 FTP -0.21635 0.0774 0.24 -2.794
## KFTno83 FTP - KFTno114 FTP -0.48470 0.0776 0.23 -6.248
## KFTno83 FTP - KFTno83 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno83 FTP - KFTno85 FTPIG 0.11548 0.0956 0.56 1.209
## KFTno83 FTP - KFTno92 FTPIG 0.00508 0.0974 0.62 0.052
## KFTno83 FTP - KFTno100 FTPIG -0.12557 0.0875 0.38 -1.436
## KFTno83 FTP - KFTno101 FTPIG -0.05992 0.0849 0.35 -0.706
## KFTno83 FTP - KFTno102 FTPIG -0.10298 0.0855 0.36 -1.205
## KFTno83 FTP - KFTno114 FTPIG -0.37134 0.0861 0.36 -4.311
## KFTno85 FTP - KFTno92 FTP -0.11040 0.0843 0.30 -1.310
## KFTno85 FTP - KFTno100 FTP -0.24104 0.0733 0.19 -3.290
## KFTno85 FTP - KFTno101 FTP -0.17540 0.0684 0.14 -2.565
## KFTno85 FTP - KFTno102 FTP -0.21846 0.0710 0.17 -3.076
## KFTno85 FTP - KFTno114 FTP -0.48682 0.0716 0.17 -6.800
## KFTno85 FTP - KFTno83 FTPIG 0.11125 0.0948 0.53 1.174
## KFTno85 FTP - KFTno85 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno85 FTP - KFTno92 FTPIG 0.00297 0.0937 0.47 0.032
## KFTno85 FTP - KFTno100 FTPIG -0.12768 0.0816 0.29 -1.565
## KFTno85 FTP - KFTno101 FTPIG -0.06204 0.0783 0.25 -0.792
## KFTno85 FTP - KFTno102 FTPIG -0.10510 0.0793 0.26 -1.326
## KFTno85 FTP - KFTno114 FTPIG -0.37345 0.0803 0.28 -4.650
## KFTno92 FTP - KFTno100 FTP -0.13065 0.0817 0.17 -1.599
## KFTno92 FTP - KFTno101 FTP -0.06500 0.0782 0.17 -0.831
## KFTno92 FTP - KFTno102 FTP -0.10806 0.0793 0.19 -1.363
## KFTno92 FTP - KFTno114 FTP -0.37642 0.0794 0.16 -4.739
## KFTno92 FTP - KFTno83 FTPIG 0.22165 0.1040 1.01 2.132
## KFTno92 FTP - KFTno85 FTPIG 0.22376 0.1010 0.84 2.210
## KFTno92 FTP - KFTno92 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno92 FTP - KFTno100 FTPIG -0.01728 0.0933 0.45 -0.185
## KFTno92 FTP - KFTno101 FTPIG 0.04836 0.0911 0.46 0.531
## KFTno92 FTP - KFTno102 FTPIG 0.00530 0.0909 0.47 0.058
## KFTno92 FTP - KFTno114 FTPIG -0.26305 0.0915 0.43 -2.876
## KFTno100 FTP - KFTno101 FTP 0.06564 0.0636 0.11 1.032
## KFTno100 FTP - KFTno102 FTP 0.02258 0.0657 0.12 0.344
## KFTno100 FTP - KFTno114 FTP -0.24577 0.0628 0.10 -3.911
## KFTno100 FTP - KFTno83 FTPIG 0.35230 0.0988 0.60 3.566
## KFTno100 FTP - KFTno85 FTPIG 0.35441 0.0944 0.51 3.754
## KFTno100 FTP - KFTno92 FTPIG 0.24401 0.0974 0.43 2.505
## KFTno100 FTP - KFTno100 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno100 FTP - KFTno101 FTPIG 0.17901 0.0814 0.28 2.199
## KFTno100 FTP - KFTno102 FTPIG 0.13595 0.0817 0.29 1.663
## KFTno100 FTP - KFTno114 FTPIG -0.13241 0.0800 0.27 -1.655
## KFTno101 FTP - KFTno102 FTP -0.04306 0.0617 61.01 -0.698
## KFTno101 FTP - KFTno114 FTP -0.31141 0.0619 0.10 -5.029
## KFTno101 FTP - KFTno83 FTPIG 0.28666 0.0947 0.52 3.026
## KFTno101 FTP - KFTno85 FTPIG 0.28877 0.0897 0.43 3.218
## KFTno101 FTP - KFTno92 FTPIG 0.17837 0.0936 0.41 1.906
## KFTno101 FTP - KFTno100 FTPIG 0.04772 0.0793 0.26 0.602
## KFTno101 FTP - KFTno101 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno101 FTP - KFTno102 FTPIG 0.07031 0.0775 61.12 0.907
## KFTno101 FTP - KFTno114 FTPIG -0.19805 0.0782 0.25 -2.532
## KFTno102 FTP - KFTno114 FTP -0.26836 0.0640 0.11 -4.194
## KFTno102 FTP - KFTno83 FTPIG 0.32971 0.0975 0.59 3.381
## KFTno102 FTP - KFTno85 FTPIG 0.33183 0.0929 0.49 3.570
## KFTno102 FTP - KFTno92 FTPIG 0.22143 0.0956 0.48 2.315
## KFTno102 FTP - KFTno100 FTPIG 0.09078 0.0823 0.31 1.103
## KFTno102 FTP - KFTno101 FTPIG 0.15642 0.0803 61.04 1.949
## KFTno102 FTP - KFTno102 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno102 FTP - KFTno114 FTPIG -0.15499 0.0812 0.29 -1.909
## KFTno114 FTP - KFTno83 FTPIG 0.59807 0.0972 0.57 6.151
## KFTno114 FTP - KFTno85 FTPIG 0.60018 0.0929 0.48 6.460
## KFTno114 FTP - KFTno92 FTPIG 0.48978 0.0953 0.41 5.139
## KFTno114 FTP - KFTno100 FTPIG 0.35914 0.0795 0.27 4.515
## KFTno114 FTP - KFTno101 FTPIG 0.42478 0.0799 0.27 5.317
## KFTno114 FTP - KFTno102 FTPIG 0.38172 0.0802 0.27 4.762
## KFTno114 FTP - KFTno114 FTPIG 0.11337 0.0492 61.18 2.306
## KFTno83 FTPIG - KFTno85 FTPIG 0.00211 0.0815 0.29 0.026
## KFTno83 FTPIG - KFTno92 FTPIG -0.10829 0.0879 0.42 -1.231
## KFTno83 FTPIG - KFTno100 FTPIG -0.23893 0.0793 0.25 -3.013
## KFTno83 FTPIG - KFTno101 FTPIG -0.17329 0.0753 0.21 -2.301
## KFTno83 FTPIG - KFTno102 FTPIG -0.21635 0.0774 0.24 -2.794
## KFTno83 FTPIG - KFTno114 FTPIG -0.48470 0.0776 0.23 -6.248
## KFTno85 FTPIG - KFTno92 FTPIG -0.11040 0.0843 0.30 -1.310
## KFTno85 FTPIG - KFTno100 FTPIG -0.24104 0.0733 0.19 -3.290
## KFTno85 FTPIG - KFTno101 FTPIG -0.17540 0.0684 0.14 -2.565
## KFTno85 FTPIG - KFTno102 FTPIG -0.21846 0.0710 0.17 -3.076
## KFTno85 FTPIG - KFTno114 FTPIG -0.48682 0.0716 0.17 -6.800
## KFTno92 FTPIG - KFTno100 FTPIG -0.13065 0.0817 0.17 -1.599
## KFTno92 FTPIG - KFTno101 FTPIG -0.06500 0.0782 0.17 -0.831
## KFTno92 FTPIG - KFTno102 FTPIG -0.10806 0.0793 0.19 -1.363
## KFTno92 FTPIG - KFTno114 FTPIG -0.37642 0.0794 0.16 -4.739
## KFTno100 FTPIG - KFTno101 FTPIG 0.06564 0.0636 0.11 1.032
## KFTno100 FTPIG - KFTno102 FTPIG 0.02258 0.0657 0.12 0.344
## KFTno100 FTPIG - KFTno114 FTPIG -0.24577 0.0628 0.10 -3.911
## KFTno101 FTPIG - KFTno102 FTPIG -0.04306 0.0617 61.01 -0.698
## KFTno101 FTPIG - KFTno114 FTPIG -0.31141 0.0619 0.10 -5.029
## KFTno102 FTPIG - KFTno114 FTPIG -0.26836 0.0640 0.11 -4.194
## p.value
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 0.5490
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 0.3888
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0612
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.0287
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: holm method for 378 tests
# Shannon diversity across methods and timepoints with both boxplots and
violin plots
library(ggplot2)
ggplot(diversity_data, aes(x = Method, y = Shannon)) +
geom_boxplot(outlier.shape = NA, position = position_dodge(0.8), alpha = 0.15) +
# color = KFTno
geom_jitter(aes(color = as.factor(KFTno)), width = 0.2, size = 3, alpha = 0.7) +
# Mean points
stat_summary(fun = mean, geom = "point", shape = 21, size = 3.5, color = "black", fill = "white") +
stat_summary(
fun.data = function(x) {
mean_val <- mean(x)
se_val <- sd(x) / sqrt(length(x))
data.frame(y = mean_val, ymin = mean_val - se_val, ymax = mean_val + se_val)
},
geom = "errorbar", width = 0.25, color = "black"
) +
# Color palette for KFTno
scale_colour_manual(values = c(
"100" = "red",
"101" = "#56B4E9",
"102" = "#009E73",
"11" = "#F0E442",
"114" = "#0072B2",
"31" = "#D55E00",
"83" = "#CC33CC",
"84" = "#CC79A7",
"85" = "#FFA07A",
"86" = "#FDE725",
"90" = "orange",
"92" = "blue",
"95" = "darkgreen"
)) +
labs(
title = "Shannon Diversity Index by Method and KFT Sample",
x = "Extraction Method",
y = "Shannon Diversity Index",
color = "KFT Sample"
) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(size = 14, face = "bold"),
legend.title = element_text(size = 12)
)
26 session info
## R version 4.5.0 (2025-04-11)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] devtools_2.4.5 usethis_3.1.0
## [3] pairwiseAdonis_0.4.1 cluster_2.1.8.1
## [5] here_1.0.1 extrafont_0.19
## [7] lubridate_1.9.4 forcats_1.0.0
## [9] stringr_1.5.1 purrr_1.0.4
## [11] readr_2.1.5 tidyr_1.3.1
## [13] tibble_3.2.1 tidyverse_2.0.0
## [15] dplyr_1.1.4 microbiome_1.30.0
## [17] emmeans_1.11.0 lmerTest_3.1-3
## [19] lme4_1.1-37 Matrix_1.7-3
## [21] DESeq2_1.48.0 SummarizedExperiment_1.38.0
## [23] Biobase_2.68.0 MatrixGenerics_1.20.0
## [25] matrixStats_1.5.0 GenomicRanges_1.60.0
## [27] GenomeInfoDb_1.44.0 IRanges_2.42.0
## [29] S4Vectors_0.46.0 BiocGenerics_0.54.0
## [31] generics_0.1.3 vegan_2.6-10
## [33] lattice_0.22-7 permute_0.9-7
## [35] ggrepel_0.9.6 decontam_1.28.0
## [37] phyloseq_1.52.0 ggplot2_3.5.2
## [39] DspikeIn_0.99.0
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.6 httr_1.4.7
## [3] RColorBrewer_1.1-3 numDeriv_2016.8-1.1
## [5] profvis_0.4.0 tools_4.5.0
## [7] backports_1.5.0 R6_2.6.1
## [9] lazyeval_0.2.2 mgcv_1.9-3
## [11] rhdf5filters_1.20.0 urlchecker_1.0.1
## [13] withr_3.0.2 gridExtra_2.3
## [15] cli_3.6.5 textshaping_1.0.0
## [17] officer_0.6.8 labeling_0.4.3
## [19] sass_0.4.10 mvtnorm_1.3-3
## [21] randomForest_4.7-1.2 ggridges_0.5.6
## [23] askpass_1.2.1 systemfonts_1.2.2
## [25] yulab.utils_0.2.0 sessioninfo_1.2.3
## [27] limma_3.64.0 rstudioapi_0.17.1
## [29] gridGraphics_0.5-1 car_3.1-3
## [31] zip_2.3.2 biomformat_1.36.0
## [33] DECIPHER_3.4.0 abind_1.4-8
## [35] lifecycle_1.0.4 yaml_2.3.10
## [37] edgeR_4.6.1 carData_3.0-5
## [39] rhdf5_2.52.0 SparseArray_1.8.0
## [41] Rtsne_0.17 grid_4.5.0
## [43] promises_1.3.2 crayon_1.5.3
## [45] miniUI_0.1.2 pillar_1.10.2
## [47] knitr_1.50 boot_1.3-31
## [49] estimability_1.5.1 codetools_0.2-20
## [51] fastmatch_1.1-6 glue_1.8.0
## [53] ggfun_0.1.8 fontLiberation_0.1.0
## [55] data.table_1.17.0 remotes_2.5.0
## [57] vctrs_0.6.5 treeio_1.32.0
## [59] Rdpack_2.6.4 gtable_0.3.6
## [61] cachem_1.1.0 xfun_0.52
## [63] rbibutils_2.3 S4Arrays_1.8.0
## [65] mime_0.13 tidygraph_1.3.1
## [67] coda_0.19-4.1 reformulas_0.4.0
## [69] survival_3.8-3 SingleCellExperiment_1.30.0
## [71] iterators_1.0.14 statmod_1.5.0
## [73] ellipsis_0.3.2 nlme_3.1-168
## [75] pbkrtest_0.5.3 ggtree_3.16.0
## [77] fontquiver_0.2.1 rprojroot_2.0.4
## [79] bslib_0.9.0 colorspace_2.1-1
## [81] DBI_1.2.3 ade4_1.7-23
## [83] phangorn_2.12.1 tidyselect_1.2.1
## [85] extrafontdb_1.0 compiler_4.5.0
## [87] curl_6.2.2 flextable_0.9.7
## [89] xml2_1.3.8 fontBitstreamVera_0.1.1
## [91] DelayedArray_0.34.1 scales_1.4.0
## [93] quadprog_1.5-8 digest_0.6.37
## [95] minqa_1.2.8 rmarkdown_2.29
## [97] XVector_0.48.0 htmltools_0.5.8.1
## [99] pkgconfig_2.0.3 fastmap_1.2.0
## [101] rlang_1.1.6 htmlwidgets_1.6.4
## [103] UCSC.utils_1.4.0 shiny_1.10.0
## [105] farver_2.1.2 jquerylib_0.1.4
## [107] jsonlite_2.0.0 BiocParallel_1.42.0
## [109] magrittr_2.0.3 Formula_1.2-5
## [111] GenomeInfoDbData_1.2.14 ggplotify_0.1.2
## [113] patchwork_1.3.0 Rhdf5lib_1.30.0
## [115] Rcpp_1.0.14 ape_5.8-1
## [117] ggnewscale_0.5.1 viridis_0.6.5
## [119] gdtools_0.4.2 stringi_1.8.7
## [121] ggstar_1.0.4 ggalluvial_0.12.5
## [123] ggraph_2.2.1 MASS_7.3-65
## [125] plyr_1.8.9 pkgbuild_1.4.7
## [127] parallel_4.5.0 Biostrings_2.76.0
## [129] graphlayouts_1.2.2 splines_4.5.0
## [131] multtest_2.64.0 hms_1.1.3
## [133] msa_1.40.0 locfit_1.5-9.12
## [135] igraph_2.1.4 ggpubr_0.6.0
## [137] uuid_1.2-1 ggtreeExtra_1.18.0
## [139] ggsignif_0.6.4 reshape2_1.4.4
## [141] pkgload_1.4.0 evaluate_1.0.3
## [143] nloptr_2.2.1 tzdb_0.5.0
## [145] foreach_1.5.2 tweenr_2.0.3
## [147] httpuv_1.6.16 Rttf2pt1_1.3.12
## [149] openssl_2.3.2 polyclip_1.10-7
## [151] ggforce_0.4.2 broom_1.0.8
## [153] xtable_1.8-4 tidytree_0.4.6
## [155] rstatix_0.7.2 later_1.4.2
## [157] viridisLite_0.4.2 ragg_1.4.0
## [159] aplot_0.2.5 memoise_2.0.1
## [161] TreeSummarizedExperiment_2.16.0 timechange_0.3.0