Untitled
GrGladkov
2023-07-18
0.1 Preprocessing:
Basic dada2 pipeline for 16S:
truncLen=c(220, 160)
trimLeft=c(19,20)
maxEE=c(2,5)
truncQ=2
pool=pseudo
database - SILVA 138.1
Delete all bad taxa with parameters below:
non assigned at the phylum level +
chloroplasts/mitochondria +
low reads/richness(based on the plot results)
0.2 Import
Import libraries and the previously generated phyloseq object
library(readxl)
library(tidyverse)
library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(ANCOMBC)
library(heatmaply)
library(compositions)
library(WGCNA)
require(phyloseq)
ps.f2 <- readRDS("ps.f2")ps.liv <- prune_samples(sample_data(ps.f2)$Group %in% "Pimpirev", ps.f2)
ps.liv <- prune_taxa(taxa_sums(ps.liv) > 0, ps.liv)
ps.liv <- prune_samples(!sample_names(ps.liv) %in% paste0("Abakumov.Krio2.",seq(9, 12)), ps.liv)
ps.liv <- prune_taxa(taxa_sums(ps.liv) > 0, ps.liv)
ps.cr <- prune_samples(sample_data(ps.f2)$Cryoconite == "cryoconite", ps.f2)
ps.cr <- prune_taxa(taxa_sums(ps.cr) > 0, ps.cr)0.3 Functions
plot_rich_reads_samlenames_lm <- function(physeq, group = "Site", label = "Repeat"){
rish <- estimate_richness(physeq, measures = "Observed")
reads.sum <- as.data.frame(sample_sums(physeq))
reads.summary <- cbind(rish, reads.sum)
colnames(reads.summary) <- c("otus","reads")
reads.summary["Repeat"] <-unlist(purrr::map(stringr::str_split(rownames(physeq@sam_data), "\\.", 2), function(x) x[[2]]))
reads.summary["Site"] <- physeq@sam_data[[group]]
library(ggrepel)
require(ggforce)
p1 <- ggplot(data=reads.summary) +
geom_point(aes(y=otus, x=log2(reads), color=Site),size=3) +
geom_text_repel(aes(y=otus, x=log2(reads), label=paste0(Repeat))) +
theme_bw() +
geom_smooth(aes(y=otus, x=log2(reads), fill=Site, color=Site),method=lm, se=FALSE, ymin = 1) +
scale_x_continuous(sec.axis = sec_axis(sec.axis ~ 2**.))
return(p1)
}
beta_custom_norm_NMDS_elli_w <- function(ps, seed = 7888, normtype="vst", Color="What", Group="Repeat"){
require(phyloseq)
require(ggplot2)
require(ggpubr)
library(ggforce)
ps@otu_table[ps@otu_table < 0] <- 0
ordination.b <- ordinate(ps, "NMDS", "bray")
mds <- as.data.frame(ordination.b$points)
p <- plot_ordination(ps,
ordination.b,
type="sample",
color = Color,
title="NMDS - Bray-Curtis",
# title=NULL,
axes = c(1,2) ) +
theme_bw() +
theme(text = element_text(size = 10)) +
geom_point(size = 3) +
annotate("text",
x=min(mds$MDS1) + abs(min(mds$MDS1))/4,
y=max(mds$MDS2),
label=paste0("Stress -- ", round(ordination.b$stress, 3))) +
geom_mark_ellipse(aes_string(group = Group, label = Group),
label.fontsize = 10,
label.buffer = unit(2, "mm"),
label.minwidth = unit(5, "mm"),
con.cap = unit(0.1, "mm"),
con.colour='gray') +
theme(legend.position = "none") +
scale_colour_viridis_d(option = "magma",
aesthetics = "color",
begin = 0,
end = 0.8)
return(p)
}
plot_alpha_w_toc_mod <- function(ps, group, metric) {
require(phyloseq)
require(ggplot2)
ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
er <- estimate_richness(ps_a)
df_er <- cbind(ps_a@sam_data, er)
df_er <- df_er %>% select(c(group, metric))
stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
rstatix::tukey_hsd() %>%
filter(p.adj.signif != "ns")
y <- seq(max(er[[metric]]), length=length(stat.test$p.adj), by=max(er[[metric]]/20))
plot_richness(ps_a, x=group, measures=metric, color="Region") +
geom_boxplot() +
geom_point(size=1.2, alpha=0.3) +
ggpubr::stat_pvalue_manual(
stat.test,
label = "p.adj.signif",
y.position = y,
tip.length = 0.005,
bracket.nudge.y = 1.5,
vjust = 0.6,
size = 3) +
theme_light() +
scale_colour_viridis_d(option = "magma",
aesthetics = "color",
begin = 0,
end = 0.8) +
theme(axis.text.x = element_text(angle = 45, hjust=1),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
legend.position = "None") +
labs(y=paste(metric, "index"))
}
phyloseq_to_ampvis2 <- function(physeq) {
#check object for class
if(!any(class(physeq) %in% "phyloseq"))
stop("physeq object must be of class \"phyloseq\"", call. = FALSE)
#ampvis2 requires taxonomy and abundance table, phyloseq checks for the latter
if(is.null(physeq@tax_table))
stop("No taxonomy found in the phyloseq object and is required for ampvis2", call. = FALSE)
#OTUs must be in rows, not columns
if(phyloseq::taxa_are_rows(physeq))
abund <- as.data.frame(phyloseq::otu_table(physeq)@.Data)
else
abund <- as.data.frame(t(phyloseq::otu_table(physeq)@.Data))
#tax_table is assumed to have OTUs in rows too
tax <- phyloseq::tax_table(physeq)@.Data
#merge by rownames (OTUs)
otutable <- merge(
abund,
tax,
by = 0,
all.x = TRUE,
all.y = FALSE,
sort = FALSE
)
colnames(otutable)[1] <- "OTU"
#extract sample_data (metadata)
if(!is.null(physeq@sam_data)) {
metadata <- data.frame(
phyloseq::sample_data(physeq),
row.names = phyloseq::sample_names(physeq),
stringsAsFactors = FALSE,
check.names = FALSE
)
#check if any columns match exactly with rownames
#if none matched assume row names are sample identifiers
samplesCol <- unlist(lapply(metadata, function(x) {
identical(x, rownames(metadata))}))
if(any(samplesCol)) {
#error if a column matched and it's not the first
if(!samplesCol[[1]])
stop("Sample ID's must be in the first column in the sample metadata, please reorder", call. = FALSE)
} else {
#assume rownames are sample identifiers, merge at the end with name "SampleID"
if(any(colnames(metadata) %in% "SampleID"))
stop("A column in the sample metadata is already named \"SampleID\" but does not seem to contain sample ID's", call. = FALSE)
metadata$SampleID <- rownames(metadata)
#reorder columns so SampleID is the first
metadata <- metadata[, c(which(colnames(metadata) %in% "SampleID"), 1:(ncol(metadata)-1L)), drop = FALSE]
}
} else
metadata <- NULL
#extract phylogenetic tree, assumed to be of class "phylo"
if(!is.null(physeq@phy_tree)) {
tree <- phyloseq::phy_tree(physeq)
} else
tree <- NULL
#extract OTU DNA sequences, assumed to be of class "XStringSet"
if(!is.null(physeq@refseq)) {
#convert XStringSet to DNAbin using a temporary file (easiest)
fastaTempFile <- tempfile(pattern = "ampvis2_", fileext = ".fa")
Biostrings::writeXStringSet(physeq@refseq, filepath = fastaTempFile)
} else
fastaTempFile <- NULL
#load as normally with amp_load
ampvis2::amp_load(
otutable = otutable,
metadata = metadata,
tree = tree,
fasta = fastaTempFile
)
}
detachAllPackages <- function() {
basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")
package.list <- search()[ifelse(unlist(gregexpr("package:",search()))==1,TRUE,FALSE)]
package.list <- setdiff(package.list,basic.packages)
if (length(package.list)>0) for (package in package.list) detach(package, character.only=TRUE, force = TRUE)
}
get_list_from_ps <- function(physeq, amaz_fac) {
my_keys <- levels(sample_data(physeq)[[amaz_fac]])
mylist <- vector(mode="list", length=length(my_keys))
names(mylist) <- my_keys
for (i in levels(sample_data(physeq)[[amaz_fac]])) {
x <- prune_samples(sample_data(physeq)[[amaz_fac]] %in% i, physeq)
x <- prune_taxa(taxa_sums(x) > 0, x)
mylist[[i]] <- taxa_names(x)
}
return(mylist)
}
unregister <- function() {
env <- foreach:::.foreachGlobals
rm(list=ls(name=env), pos=env)
}
lsf.str()## beta_custom_norm_NMDS_elli_w : function (ps, seed = 7888, normtype = "vst", Color = "What", Group = "Repeat")
## detachAllPackages : function ()
## get_list_from_ps : function (physeq, amaz_fac)
## phyloseq_to_ampvis2 : function (physeq)
## plot_alpha_w_toc_mod : function (ps, group, metric)
## plot_rich_reads_samlenames_lm : function (physeq, group = "Site", label = "Repeat")
## unregister : function ()plot_rich_reads_samlenames_lm(ps.f2, group = "Group", label = "Group") 
0.4 Alpha
ps.f2.al <- ps.f2
ps.f2.al@sam_data <- ps.f2.al@sam_data %>%
data.frame() %>%
mutate(Group=forcats::fct_relevel(Group, c("Pimpirev","IGAN", "Mush" ,"Garabashi", "Shkelda","Mud", "Soil"))) %>%
sample_data()
p_a1 <- plot_alpha_w_toc_mod(ps.f2.al, group="Group", metric="Observed")
p_a2 <- plot_alpha_w_toc_mod(ps.f2.al, group="Group", metric="Shannon")
p_a3 <- plot_alpha_w_toc_mod(ps.f2.al, group="Group", metric="InvSimpson")
p_alpha <- ggpubr::ggarrange(p_a1, p_a2, p_a3, nrow = 1)
p_alpha
estimate_richness(ps.f2) %>%
cbind2(sample_data(ps.f2) %>%
data.frame() %>%
select(Type)) %>%
relocate(Type) %>%
group_by(Type) %>%
select(c("Observed", "Shannon", "InvSimpson")) %>%
summarise_all(mean) %>%
DT::datatable(caption = "Alpha diversity metrics")0.5 Beta
clr transformation from compositions library
https://www.rdocumentation.org/packages/compositions/versions/2.0-6/topics/clr
ps.f2.clr <- ps.f2
clr.f2.otu <- as.data.frame(t(ps.f2@otu_table)) %>%
compositions::clr() %>%
t()
ps.f2.clr@otu_table <- otu_table(clr.f2.otu, taxa_are_rows = FALSE)
beta_all <- beta_custom_norm_NMDS_elli_w(ps.f2.clr, Group = "Group", Color = "Region") +
# labs(title = "A. all samples") +
labs(title = "NMDS - Bray-Curtis") +
scale_colour_viridis_d(option = "magma",
aesthetics = "color",
begin = 0,
end = 0.8)## Run 0 stress 0.1541643
## Run 1 stress 0.1234039
## ... New best solution
## ... Procrustes: rmse 0.08175852 max resid 0.3218974
## Run 2 stress 0.1486047
## Run 3 stress 0.2177005
## Run 4 stress 0.2140094
## Run 5 stress 0.1234039
## ... Procrustes: rmse 0.00001325593 max resid 0.00009014157
## ... Similar to previous best
## Run 6 stress 0.1817778
## Run 7 stress 0.1840804
## Run 8 stress 0.1234039
## ... New best solution
## ... Procrustes: rmse 0.000006835381 max resid 0.00004115244
## ... Similar to previous best
## Run 9 stress 0.1746337
## Run 10 stress 0.1234039
## ... Procrustes: rmse 0.00001005475 max resid 0.00006153985
## ... Similar to previous best
## Run 11 stress 0.1234039
## ... Procrustes: rmse 0.00001306177 max resid 0.00008105919
## ... Similar to previous best
## Run 12 stress 0.1234039
## ... Procrustes: rmse 0.00002374299 max resid 0.0001604382
## ... Similar to previous best
## Run 13 stress 0.1975259
## Run 14 stress 0.1234039
## ... Procrustes: rmse 0.0000109071 max resid 0.00005281759
## ... Similar to previous best
## Run 15 stress 0.1234039
## ... New best solution
## ... Procrustes: rmse 0.000003859064 max resid 0.00002262156
## ... Similar to previous best
## Run 16 stress 0.1234039
## ... Procrustes: rmse 0.00002133015 max resid 0.0001403891
## ... Similar to previous best
## Run 17 stress 0.1841554
## Run 18 stress 0.1736377
## Run 19 stress 0.1922535
## Run 20 stress 0.2033537
## *** Best solution repeated 2 timesbeta_all
0.5.1 heatmaps
phylum
amp_f2 <- phyloseq_to_ampvis2(ps.f2)
amp_heatmap(amp_f2,
tax_show = 13,
group_by = "Type",
facet_by = "Region",
tax_aggregate = "Phylum",
tax_add = "Kingdom",
normalise=TRUE,
plot_values_size = 5.5,
showRemainingTaxa = TRUE,
color_vector = colorspace::lighten(viridis::magma(n = 3, begin = 0.35, direction = -1), 0.2)) +
theme(axis.text.x = element_text(size=14, angle=45, hjust=0.4),
axis.text.y = element_text(size=14)) 
Genus
amp_heatmap(amp_f2,
tax_show = 40,
group_by = "Group",
facet_by = "Region",
tax_aggregate = "Genus",
tax_add = "Phylum",
normalise=TRUE,
plot_values_size = 4.5,
showRemainingTaxa = TRUE,
color_vector = colorspace::lighten(viridis::magma(n = 3, begin = 0.35, direction = -1), 0.2)) +
theme(axis.text.x = element_text(size=14, angle=45, hjust=0.4),
axis.text.y = element_text(size=13)) 
0.5.2 agrochemestry
ps.f2@sam_data %>%
data.frame() %>%
# rownames_to_column("shit") %>%
# column_to_rownames("Name") %>%
select(which(sapply(., is.numeric)), Group) %>%
select(-c("lat", "long")) %>%
group_by(Group) %>%
summarise(across(everything(), mean),
.groups = 'drop') %>%
column_to_rownames("Group") %>%
scale() %>%
heatmaply::ggheatmap(color = viridis::magma(n = 3, begin = 0.1, direction = -1)) 
0.6 Pimpirev
A small stand-alone investigation of the ornitogenic and nonornitogenic samples from Antarctica
ancom_liv_asv <- ANCOMBC::ancombc2(data = ps.liv,
fix_formula = "Type",
rand_formula = NULL,
p_adj_method = "fdr",
pseudo = 0,
pseudo_sens = TRUE,
prv_cut = 0.10,
s0_perc = 0.05,
group = "Type",
struc_zero = TRUE,
neg_lb = TRUE,
alpha = 0.05,
n_cl = 20,
tax_level = NULL,
verbose = TRUE,
global = TRUE,
pairwise = TRUE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-2, max_iter = 20, verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
nrow = 2,
byrow = TRUE),
matrix(c(-1, 0, 1, -1),
nrow = 2,
byrow = TRUE)),
node = list(2, 2),
solver = "ECOS",
B = 100)
)amp_liv <- phyloseq_to_ampvis2(ps.liv)
amp_heatmap(amp_liv,
tax_show = 40,
group_by = "Type",
tax_aggregate = "Species",
tax_add = "Genus",
normalise=TRUE,
showRemainingTaxa = TRUE)
p.venn.cc <- amp_venn(amp_liv,
group_by = "Type",
normalise = TRUE,
cut_a = 0.000001,
cut_f=0.000001
)
p.venn.cc
Library normalization. The procedure is described in
http://www.bioconductor.org/packages/release/bioc/vignettes/ANCOMBC/inst/doc/ANCOMBC.html
chapter "Bias-corrected abundances"
An article with theory(the math is kind of sapplemented there):
Lin, H., Peddada, S.D. Analysis of compositions of microbiomes with bias correction. Nat Commun 11, 3514 (2020). https://doi.org/10.1038/s41467-020-17041-7
ps.anc.liv - phyloseq object with corrected reads( and log-transformation also)
use euclidean distance for plot
out.liv <- ANCOMBC::ancombc2(data=ps.liv,
fix_formula = "Type",
prv_cut = 0,
p_adj_method = "BH",
tax_level = NULL,
rand_formula= NULL,
pseudo = 0,
pseudo_sens = FALSE,
s0_perc = 0.05,
group = "Type",
struc_zero = FALSE,
neg_lb = TRUE,
alpha = 0.05,
n_cl = 20,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE
)
ps.anc.liv <- ps.liv
samp_frac <- out.liv$samp_frac
samp_frac[is.na(samp_frac)] <- 0
# Add pesudo-count (1) to avoid taking the log of 0
log_obs_abn = log(out.liv$feature_table + 1)
# Adjust the log observed abundances
log_corr_abn = t(t(log_obs_abn) - samp_frac)
otu_table(ps.anc.liv) <- otu_table(t(log_corr_abn), taxa_are_rows = FALSE)
beta_custom_norm_NMDS_elli_w_eu <- function(ps, seed = 7888, normtype="vst", Color="What", Group="Repeat"){
require(phyloseq)
require(ggplot2)
require(ggpubr)
library(ggforce)
ps@otu_table[ps@otu_table < 0] <- 0
ordination.b <- ordinate(ps, "NMDS", "euclidean")
mds <- as.data.frame(ordination.b$points)
p <- plot_ordination(ps,
ordination.b,
type="sample",
color = Color,
title="NMDS - euclidean",
# title=NULL,
axes = c(1,2) ) +
theme_bw() +
theme(text = element_text(size = 10)) +
geom_point(size = 3) +
annotate("text",
x=min(mds$MDS1) + abs(min(mds$MDS1))/4,
y=max(mds$MDS2),
label=paste0("Stress -- ", round(ordination.b$stress, 3))) +
geom_mark_ellipse(aes_string(group = Group, label = Group),
label.fontsize = 10,
label.buffer = unit(2, "mm"),
label.minwidth = unit(5, "mm"),
con.cap = unit(0.1, "mm"),
con.colour='gray') +
theme(legend.position = "none") +
scale_colour_viridis_d(option = "magma",
aesthetics = "color",
begin = 0,
end = 0.8)
return(p)
}
beta_custom_norm_NMDS_elli_w_eu(ps.anc.liv, Color = "Type", Group = "Type")## Run 0 stress 0.0380671
## Run 1 stress 0.1810879
## Run 2 stress 0.2055201
## Run 3 stress 0.0380671
## ... Procrustes: rmse 0.0000114568 max resid 0.00003143128
## ... Similar to previous best
## Run 4 stress 0.0380671
## ... Procrustes: rmse 0.00001192981 max resid 0.00003103653
## ... Similar to previous best
## Run 5 stress 0.03806709
## ... New best solution
## ... Procrustes: rmse 0.000001311323 max resid 0.000003218826
## ... Similar to previous best
## Run 6 stress 0.1921644
## Run 7 stress 0.0380671
## ... Procrustes: rmse 0.00002206572 max resid 0.00006009745
## ... Similar to previous best
## Run 8 stress 0.03806709
## ... New best solution
## ... Procrustes: rmse 0.000003322289 max resid 0.00000867924
## ... Similar to previous best
## Run 9 stress 0.0380671
## ... Procrustes: rmse 0.000008272935 max resid 0.00002251381
## ... Similar to previous best
## Run 10 stress 0.03806709
## ... Procrustes: rmse 0.000002921055 max resid 0.000007264082
## ... Similar to previous best
## Run 11 stress 0.0380671
## ... Procrustes: rmse 0.000008032725 max resid 0.00002208769
## ... Similar to previous best
## Run 12 stress 0.03806709
## ... Procrustes: rmse 0.000002114938 max resid 0.00000549362
## ... Similar to previous best
## Run 13 stress 0.20078
## Run 14 stress 0.2007796
## Run 15 stress 0.0380671
## ... Procrustes: rmse 0.000007948323 max resid 0.00002166962
## ... Similar to previous best
## Run 16 stress 0.2454791
## Run 17 stress 0.03806709
## ... Procrustes: rmse 0.000002244759 max resid 0.000004699939
## ... Similar to previous best
## Run 18 stress 0.2154322
## Run 19 stress 0.2357637
## Run 20 stress 0.0380671
## ... Procrustes: rmse 0.000002395101 max resid 0.000004875969
## ... Similar to previous best
## *** Best solution repeated 8 times
permanova
dist <- phyloseq::distance(ps.anc.liv, "euclidean")
metadata <- as(sample_data(ps.anc.liv@sam_data), "data.frame")
vegan::adonis2(dist ~ Type, data = metadata)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = dist ~ Type, data = metadata)
## Df SumOfSqs R2 F Pr(>F)
## Type 1 3351.2 0.26948 5.1644 0.002 **
## Residual 14 9084.7 0.73052
## Total 15 12435.9 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1DA visualization
taxa_family <- ps.anc.liv %>%
tax_table() %>%
as.data.frame() %>%
rownames_to_column("ID") %>%
add_column(abnd = taxa_sums(ps.anc.liv))
ancom_liv_asv$res %>%
filter(diff_TypePAOr == TRUE) %>%
select(c("taxon", "lfc_TypePAOr")) %>%
rename("ID" = "taxon") %>%
left_join(taxa_family, by="ID") %>%
mutate(Phylum = case_when(Phylum %in% "Pseudomonadota" ~ Class,
!Phylum %in% "Pseudomonadota" ~ Phylum) %>% as.factor()) %>%
mutate(`more in` = ifelse(lfc_TypePAOr < 0, "PAG", "PAOr") %>% as.factor()) %>%
group_by(Genus, `more in`, Phylum ) %>%
summarise(area = sum(abnd)) %>%
mutate(Genus = str_replace_na(Genus, replacement = ' ')) %>%
ggplot(aes(area = area,
fill = `more in`,
subgroup = Phylum,
label = Genus)) +
treemapify::geom_treemap() +
treemapify::geom_treemap_subgroup_border(color = "white") +
treemapify::geom_treemap_subgroup_text(place = "centre",
grow = T,
alpha = 0.5,
colour = "black",
fontface = "italic",
min.size = 0) +
treemapify::geom_treemap_text(colour = "white",
place = "topleft",
grow = T,
reflow = T,
layout = 'squarified') +
scale_colour_viridis_d(option = "magma",
aesthetics = "fill",
begin = 0.4,
end = 0.8) +
theme(legend.position="bottom") 
0.7 Core
Venn without any filtrations
code from https://github.com/a-zverev/16s-amplicon-processing
library(ggVennDiagram)
plot_vienn <- function(ps, group){
physeq <- prune_taxa(taxa_sums(ps) > 0, ps)
groups <- levels(sample_data(physeq)[[group]] %>% as.factor())
data <- merge_samples(physeq, group) %>%
psmelt() %>%
group_by(Sample, OTU) %>%
summarise(ASVs_abund = list(paste(OTU, 1:sum(Abundance))), Abund = sum(Abundance), .groups='keep') %>%
filter(Abund > 0)
asvs <- data %>% select(Sample, OTU) %>% group_by(Sample) %>% summarise(ASVs = list(OTU)) %>% as.list()
d1 <- asvs[[2]]
names(d1) <- asvs[[1]]
d1
weighted.asvs <- data %>% select(Sample, ASVs_abund) %>% group_by(Sample) %>% summarise(ASVs = list(unlist(ASVs_abund)))
d2 <- weighted.asvs[[2]]
names(d2) <- weighted.asvs[[1]]
d2
list(ggVennDiagram(d1) + ggtitle("ASVs") + scale_fill_distiller(palette = "OrRd", trans = "reverse"),
ggVennDiagram(d2) + ggtitle("Reads") + scale_fill_distiller(palette = "OrRd", trans = "reverse"))
}
plot_vienn(ps.cr, "Group")## [[1]]
##
## [[2]]
Three sampling region core phylotypes.
physeq <- ps.cr
l <- get_list_from_ps(physeq, "Group")
core <- Reduce(intersect, l)
ps.core <- prune_taxa(core, ps.cr)
length(core)## [1] 152Find structural zeros based on function with ancom-bc package
ancom_core <- ANCOMBC::ancombc2(data = ps.core,
fix_formula = "Group",
rand_formula = NULL,
p_adj_method = "fdr",
pseudo = 0,
pseudo_sens = TRUE,
prv_cut = 0.10,
s0_perc = 0.05,
group = "Group",
struc_zero = TRUE,
neg_lb = TRUE,
alpha = 0.05,
n_cl = 15,
tax_level = NULL,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-2, max_iter = 20, verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
nrow = 2,
byrow = TRUE),
matrix(c(-1, 0, 1, -1),
nrow = 2,
byrow = TRUE)),
node = list(2, 2),
solver = "ECOS",
B = 1)
)Only 32 phylotypes remain(from 152) for the core community after filtering based on the concept of structural zeros.
non_str_otu <- ancom_core$res$taxon
zeros <- ancom_core$zero_ind %>%
column_to_rownames("taxon") %>%
rowSums() > 0L
zeros.all <- zeros[zeros == FALSE]
zeros.all <- names(zeros.all)
zeros.all %>%
length()## [1] 32Bottom of the list - Lactobacillus apis, Frischella perrara, Snodgrassella alvi - bees gut endosymbionts, so maybe contamination from another dataset (try to find in documentation for our Illumina Miseq run something related to bees, honey of something - didn't find anything specific, but I think glaciers and the Antarctic are not regular vomiting places for bees).
At the top of the list - the "true core". Validate through https://www.gbif.org/ databases whether these phylotypes are really typical for glaciers. Eliminate phylotypes that look like a cross-contamination from different samples (like Seq1 (I think in real it was only in Mushketova place), Seq35 (Pipmperev reads looks like error) etc.) I think in future I should make threshold in structural zeros a little bit more strict (but I should double check this).
So I'm going to select only 4 phylotypes sequentially: Polaromonas sp.(Seq41), Cryobacterium sp.(Seq82), Rhodoferax sp.(Seq17) and Hymenobacter frigidus(Seq81). This is quite subjective, but I'm of the opinion that regular aplicon sequencing is like a leaky sieve. And I have some articles that prove that -
https://doi.org/10.1128/msystems.00186-19
https://doi.org/10.1128%2FmBio.00598-21
heatmap and species names for validation
ps.core.ancm <- prune_taxa(zeros.all, ps.core)
amp.core.ancm <- phyloseq_to_ampvis2(ps.core.ancm)
amp.core.ancm$tax %>%
select(c(Genus, Species))## Genus Species
## Seq1 Phormidesmis_ANT.L52.6 priestleyi
## Seq3 Tychonema_CCAP_1459-11B
## Seq4 Phormidium_CYN64
## Seq6 Granulicella
## Seq7 Tychonema_CCAP_1459-11B
## Seq8 Parafrigoribacterium amurskyense
## Seq12 Acidiphilium
## Seq13
## Seq16 Polaromonas
## Seq17 Rhodoferax
## Seq21 Blastocatella
## Seq24 Clostridium_sensu_stricto_9
## Seq28 Pseudomonas
## Seq35 Ferruginibacter
## Seq41 Polaromonas
## Seq43 Cutibacterium
## Seq46
## Seq48 Sphingomonas jaspsi
## Seq49 Sphingomonas
## Seq67 AAP99
## Seq74 Gemmatimonas
## Seq81 Hymenobacter frigidus
## Seq82 Cryobacterium
## Seq124 Pajaroellobacter
## Seq166 Hymenobacter
## Seq251 Pajaroellobacter
## Seq309 Massilia eurypsychrophila
## Seq548 Lactobacillus
## Seq864 Snodgrassella alvi
## Seq877 Frischella perrara
## Seq885
## Seq1123 Lactobacillus apisamp_heatmap(amp.core.ancm,
tax_show = 155,
tax_aggregate = "OTU",
tax_add = "Genus",
normalise=FALSE,
plot_values_size = 2,
facet_by = "Group",
showRemainingTaxa = TRUE)
0.8 Phylums DA
ancom_cr_phylum <- ANCOMBC::ancombc2(data = ps.cr,
tax_level = "Phylum",
fix_formula = "Group",
rand_formula = NULL,
p_adj_method = "fdr",
pseudo = 0,
pseudo_sens = TRUE,
prv_cut = 0.10,
s0_perc = 0.05,
group = "Group",
struc_zero = TRUE,
neg_lb = TRUE,
alpha = 0.05,
n_cl = 20,
verbose = TRUE,
global = TRUE,
pairwise = TRUE,
dunnet = TRUE,
trend = FALSE,
iter_control = list(tol = 1e-2, max_iter = 20, verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
nrow = 2,
byrow = TRUE),
matrix(c(-1, 0, 1, -1),
nrow = 2,
byrow = TRUE)),
node = list(2, 2),
solver = "ECOS",
B = 100)
)ancom_cr_phylum$res_pair %>%
select("taxon"|starts_with("p_")) %>%
rename_all(~ stringr::str_replace(., regex("p_Group", ignore_case = TRUE), "")) %>%
rename_all(~ stringr::str_replace(., regex("Group", ignore_case = TRUE), "")) %>%
rename_if( !stringr::str_detect(names(.), "_") & !stringr::str_detect(names(.), "taxon"), ~ paste0(., "_Garabashi")) %>%
reshape2::melt() %>%
rstatix::p_format(value, new.col = TRUE, digits = 2, accuracy = 1e-02) %>%
mutate(value.format = ifelse(value.format == "1", "", value.format)) %>%
separate(variable,into=c('Var1', 'Var2'), sep = "_") %>%
mutate_if(is.character, as.factor) %>%
ggplot(aes(Var1, Var2, fill = value)) +
geom_tile() +
scale_fill_gradientn(colours = c("#A5317EFF", "#FCFDBFFF", "#FCFDBFFF" ),
values = scales::rescale(c(0, 0.05, 1))) +
geom_text(aes(label = value.format)) +
facet_wrap(~taxon) +
theme_bw() +
theme(legend.position = "none",
axis.title.x=element_blank(),
axis.title.y=element_blank())
0.9 WGCNA
I'm trying to do variance stabilization with DESeq2 and clr from compositions before using the ancom-bc based normalization. Some of the intermediate plots are shown below.
vst
A caption
clr
0.9.1 ancom-bc
Filtration of minor phylotypes
ps.cr.f <- phyloseq::filter_taxa(ps.cr, function(x) sum(x > 10) > (0.1*length(x)), TRUE)
pruned <- setdiff(phyloseq::taxa_names(ps.cr.f), c("Seq396", "Seq321", "Seq66"))
ps.cr.f <- phyloseq::prune_taxa(pruned, ps.cr)out.f <- ANCOMBC::ancombc2(data=ps.cr.f,
fix_formula = "Group",
prv_cut = 0,
p_adj_method = "BH",
tax_level = NULL,
rand_formula= NULL,
pseudo = 0,
pseudo_sens = FALSE,
s0_perc = 0.05,
group = "Group",
struc_zero = FALSE,
neg_lb = TRUE,
alpha = 0.05,
n_cl = 20,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE
)ps.anc.f <- ps.cr.f
samp_frac <- out.f$samp_frac
samp_frac[is.na(samp_frac)] <- 0
# Add pesudo-count (1) to avoid taking the log of 0
log_obs_abn = log(out.f$feature_table + 1)
# Adjust the log observed abundances
log_corr_abn = t(t(log_obs_abn) - samp_frac)
otu_table(ps.anc.f) <- otu_table(t(log_corr_abn), taxa_are_rows = FALSE)
beta_custom_norm_NMDS_elli_w(ps.anc.f, Color = "Group", Group = "Group")## Wisconsin double standardization
## Run 0 stress 0.1012824
## Run 1 stress 0.1016154
## ... Procrustes: rmse 0.04439897 max resid 0.1066093
## Run 2 stress 0.1016154
## ... Procrustes: rmse 0.04440395 max resid 0.1066519
## Run 3 stress 0.1012824
## ... New best solution
## ... Procrustes: rmse 0.000006665935 max resid 0.00001330675
## ... Similar to previous best
## Run 4 stress 0.1493442
## Run 5 stress 0.1495521
## Run 6 stress 0.1701576
## Run 7 stress 0.1802249
## Run 8 stress 0.1675928
## Run 9 stress 0.1701001
## Run 10 stress 0.1674537
## Run 11 stress 0.1012824
## ... Procrustes: rmse 0.000005410664 max resid 0.00002359805
## ... Similar to previous best
## Run 12 stress 0.1012825
## ... Procrustes: rmse 0.00005489278 max resid 0.0003545365
## ... Similar to previous best
## Run 13 stress 0.1012824
## ... Procrustes: rmse 0.000002814839 max resid 0.00001620274
## ... Similar to previous best
## Run 14 stress 0.1012824
## ... Procrustes: rmse 0.000009900771 max resid 0.00005448954
## ... Similar to previous best
## Run 15 stress 0.1012824
## ... Procrustes: rmse 0.000003532457 max resid 0.00001419618
## ... Similar to previous best
## Run 16 stress 0.1012824
## ... Procrustes: rmse 0.000005205455 max resid 0.00003132901
## ... Similar to previous best
## Run 17 stress 0.1016154
## ... Procrustes: rmse 0.04440475 max resid 0.1066483
## Run 18 stress 0.1763084
## Run 19 stress 0.1016154
## ... Procrustes: rmse 0.04440135 max resid 0.106622
## Run 20 stress 0.1012824
## ... Procrustes: rmse 0.000006228778 max resid 0.00002047491
## ... Similar to previous best
## *** Best solution repeated 8 times
ps.vst.mod <- ps.anc.f
ps.vst.mod@sam_data <- ps.vst.mod@sam_data %>%
data.frame() %>%
mutate(Repeat = paste0(Group, "_", sapply(str_split(row.names(.), '\\.', 3), function(x) x[[3]]))
) %>%
sample_data()
data3 <- ps.vst.mod@otu_table@.Data %>%
as.data.frame()
rownames(data3) <- as.character(ps.vst.mod@sam_data$Repeat)
powers <- c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
unregister()
sft3 <- pickSoftThreshold(data3, powerVector = powers, verbose = 5, networkType = "signed hybrid")## pickSoftThreshold: will use block size 716.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 716 of 716
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1.0 0.689000 2.3400 0.654 123.00 123.00 183.0
## 2 1.5 0.849000 1.6900 0.806 96.60 93.40 142.0
## 3 2.0 0.810000 1.2100 0.840 78.90 75.80 130.0
## 4 2.5 0.616000 0.7670 0.938 66.40 62.50 122.0
## 5 3.0 0.231000 0.3690 0.802 57.10 54.60 115.0
## 6 3.5 0.000695 -0.0161 0.583 50.00 48.10 109.0
## 7 4.0 0.164000 -0.2380 0.601 44.30 42.40 104.0
## 8 4.5 0.396000 -0.3960 0.640 39.80 37.20 99.3
## 9 5.0 0.551000 -0.4900 0.711 36.00 32.70 95.1
## 10 5.5 0.676000 -0.5880 0.777 32.80 29.10 91.4
## 11 6.0 0.728000 -0.6390 0.792 30.10 26.40 87.9
## 12 6.5 0.752000 -0.7040 0.790 27.80 23.50 84.8
## 13 7.0 0.745000 -0.7460 0.769 25.80 21.10 81.9
## 14 7.5 0.760000 -0.7940 0.760 24.00 19.10 79.2
## 15 8.0 0.840000 -0.7920 0.874 22.50 17.30 76.7
## 16 8.5 0.869000 -0.7960 0.900 21.10 15.80 74.4
## 17 9.0 0.894000 -0.8210 0.923 19.90 14.40 72.2
## 18 9.5 0.909000 -0.8490 0.926 18.70 13.20 70.1
## 19 10.0 0.896000 -0.8580 0.900 17.70 12.00 68.2
## 20 11.0 0.788000 -0.9160 0.756 16.00 10.30 64.7
## 21 12.0 0.866000 -0.9110 0.860 14.50 8.80 61.6
## 22 13.0 0.918000 -0.9160 0.936 13.20 7.55 58.7
## 23 14.0 0.939000 -0.9270 0.963 12.10 6.57 56.2
## 24 15.0 0.949000 -0.9490 0.980 11.20 5.75 53.8
## 25 16.0 0.952000 -0.9600 0.981 10.40 5.10 51.6
## 26 17.0 0.955000 -0.9720 0.984 9.65 4.38 49.6
## 27 18.0 0.911000 -0.9920 0.935 9.00 3.82 47.7
## 28 19.0 0.915000 -1.0100 0.940 8.42 3.44 46.0
## 29 20.0 0.901000 -1.0300 0.920 7.89 3.19 44.4plot(sft3$fitIndices[,1], -sign(sft3$fitIndices[,3])*sft3$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft3$fitIndices[,1], -sign(sft3$fitIndices[,3])*sft3$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")
detachAllPackages()
library(WGCNA)
net3 <- WGCNA::blockwiseModules(data3,
power=9.5,
TOMType="signed",
networkType="signed hybrid",
nThreads=15)
library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)
mergedColors2 <- net3$colors
plotDendroAndColors(
net3$dendrograms[[1]],
mergedColors2[net3$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE,
hang = 0.03,
addGuide = TRUE,
guideHang = 0.05)
ids <- ps.anc.f@sam_data %>%
data.frame() %>%
rownames_to_column("ID") %>%
# filter(Region == "Yakutsk") %>%
pull(ID)
agra_data <- ps.anc.f@sam_data %>%
data.frame() %>%
select_if(., is.double) %>%
select( - c("long", "lat")) %>%
relocate(pH, P2O5) %>%
relocate(Ni, .after = N.NO3)
# mergedColors2 %>% unique()
nOTUs <- ncol(data3)
nSamples <- nrow(data3)
# Recalculate MEs with color labels
MEs0 <- moduleEigengenes(data3, mergedColors2)$eigengenes
MEs <- orderMEs(MEs0)
names(MEs) <- substring(names(MEs), 3)
names(MEs) <- c("black", "blue","green", "pink", "red", "turq.", "brown", "yellow", "grey")
moduleTraitCor <- cor(MEs, agra_data, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
# PLOT
sizeGrWindow(10,6)
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", format.pval(moduleTraitPvalue, digits = 2, eps = 0.001, nsmall = 3), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)
par(mar = c(4, 10, 2, 2))
# Display the correlation values within a heatmap
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(agra_data),
yLabels = names(MEs),
ySymbols = names(MEs),
font.lab.y = 2,
font.lab.x = 2,
# yLabelsPosition = "right",
cex.lab.y = 0.85,
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.8,
zlim = c(-1,1),
plotLegend = FALSE,
main = paste("Currelation of the WGCNA clusters with soil characteristics"))0.10 CCA
0.10.1 groups
Select "groups" from WGCNA clusters
PH - arctic - "brown", "yellow"
ME - caucasian - "black", "blue", "green", "pink", "red"
CN - antarctic - "turquoise"
"gray" cluster dropped
library(ggrepel)
ps.rel <- phyloseq::transform_sample_counts(ps.cr, function(x) x / sum(x) * 100)
clust.ph <- net3$colors[net3$colors %in% c("brown", "yellow")]
clust.me <- net3$colors[net3$colors %in% c("black", "blue", "green", "pink", "red")]
clust.cn <- net3$colors[net3$colors %in% c("turquoise")]
ps.ph <- prune_taxa(names(clust.ph), ps.anc.f)
ps.me <- prune_taxa(names(clust.me), ps.anc.f)
ps.cn <- prune_taxa(names(clust.cn), ps.anc.f)
tx.cn <- ps.cn %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Antarctic", length(taxa_names(ps.cn)))) %>%
add_column(abnd = taxa_sums(ps.rel)[taxa_names(ps.cn)])
tx.me <- ps.me %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Arctic", length(taxa_names(ps.me)))) %>%
add_column(abnd = taxa_sums(ps.rel)[taxa_names(ps.me)])
tx.ph <- ps.ph %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Caucasus", length(taxa_names(ps.ph)))) %>%
add_column(abnd = taxa_sums(ps.rel)[taxa_names(ps.ph)])
tx <- rbind(tx.cn, tx.me, tx.ph)
physeq <- ps.anc.f
veganifyOTU <- function(physeq){
require(phyloseq)
if(taxa_are_rows(physeq)){physeq <- t(physeq)}
return(as(otu_table(physeq), "matrix"))
}
otus.ps.vegan <- veganifyOTU(physeq)
rownames(otus.ps.vegan) <- physeq@sam_data$Group
metadata <- as(sample_data(physeq), "data.frame")
cca_w_varstab_asv <- vegan::cca(otus.ps.vegan ~ pH + N.NH4 + TOC + P2O5 + Cu + Pb + Zn + Cd + Ni + N.NO3, data=metadata)
tx.ids <- tx %>%
rownames_to_column("ID") %>%
mutate(Score = "sites")
cca.meta <- cca_w_varstab_asv$CCA$biplot %>%
as.data.frame() %>%
select(c(CCA1, CCA2)) %>%
rownames_to_column("ID") %>%
mutate(Score = "biplot")
fdat_amazing <- as.data.frame(cca_w_varstab_asv$CCA$v) %>%
select(c(CCA1, CCA2)) %>%
rownames_to_column("ID") %>%
right_join(tx.ids, by = join_by(ID)) %>%
mutate(label_size = 0) %>%
bind_rows(cca.meta) %>%
mutate(Group = as.factor(Group))
cca_clust <- ggplot(fdat_amazing) +
geom_point(data = fdat_amazing %>% dplyr::filter(Score == "sites"),
mapping = aes(x=CCA1,
y=CCA2,
color=Group,
size=abnd),
alpha=0.5) +
geom_segment(data = fdat_amazing %>% dplyr::filter(Score == "biplot"),
aes(x = 0,
xend = CCA1,
y = 0,
yend = CCA2),
alpha=0.8,
color = "red",
arrow = arrow(angle = 3)) +
geom_text_repel(data = fdat_amazing %>% dplyr::filter(Score == "biplot"),
aes(x=CCA1,
y=CCA2,
label= ID),
size=5) +
xlab(paste0(round(cca_w_varstab_asv$CCA$eig[1], 3)*100, "% CCA1")) +
ylab(paste0(round(cca_w_varstab_asv$CCA$eig[2], 3)*100, "% CCA2")) +
grids(linetype = "dashed") +
geom_vline(xintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
geom_hline(yintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
theme(legend.position = "none",
panel.background = element_rect(fill = "white", colour = "grey50")) +
scale_colour_viridis_d(option = "magma",
aesthetics = "color",
begin = 0.2,
end = 0.8) +
labs(title = "B. ASVs - WGCNA clusters")
cca_clust
0.10.2 samples
physeq <- ps.anc.f
otus.ps.vegan <- veganifyOTU(physeq)
rownames(otus.ps.vegan) <- physeq@sam_data$Group
metadata <- as(sample_data(physeq), "data.frame")
cca <- vegan::cca(otus.ps.vegan ~ pH + N.NH4 + TOC + P2O5 + Cu + Pb + Zn + Cd + Ni + N.NO3 + K2O, data=metadata)
biplot <- as.data.frame(cca$CCA$biplot)
wa <- as.data.frame(cca$CCA$wa)
biplot <- rownames_to_column(biplot, "Label") %>%
add_column(Score = rep("biplot", length(rownames(biplot)))) %>%
mutate(
Group = rep("biplot", length(rownames(biplot))),
Type = rep("biplot", length(rownames(biplot))),
Region = rep("biplot", length(rownames(biplot)))
)
wa <- rownames_to_column(wa, "Label") %>%
add_column(Score = rep("sites", length(rownames(wa)))) %>% mutate(
Group = metadata$Type,
Type = metadata$Group,
Region = metadata$Region
)
fdat_amazing <- rbind(biplot, wa)
fdat_amazing <- fdat_amazing %>%
mutate(label_size = ifelse(Score == "biplot", 4.5, 3)) %>%
mutate_if(is.character, as.factor)
cca_sampl <- ggplot(fdat_amazing) +
geom_vline(xintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
geom_hline(yintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
geom_point(data = fdat_amazing %>%
dplyr::filter(Score == "sites"), mapping = aes(x=CCA1, y=CCA2, colour = factor(Score))) +
geom_segment(data = fdat_amazing %>%
dplyr::filter(Score == "biplot"),
aes(x = 0, xend = CCA1, y = 0, yend = CCA2),
size = 0.6,
alpha=0.8,
color = "red",
arrow = arrow(angle = 3)) +
ggforce::geom_mark_ellipse(data = fdat_amazing %>%
dplyr::filter(Score == "sites"),
aes(x=CCA1, y=CCA2, group = Type, label = Type, col = 'grey'),
label.fontsize = 12,
label.buffer = unit(2, "mm"),
label.minwidth = unit(5, "mm"),
con.cap = unit(0.1, "mm"),
con.colour='salmon',
label.colour='salmon') +
ggrepel::geom_text_repel(data = fdat_amazing %>%
dplyr::filter(Score == "biplot"), aes(x=CCA1, y=CCA2, label= Label), size=5) +
xlab(paste0(round(cca$CCA$eig[1], 3)*100, "% CCA1")) +
ylab(paste0(round(cca$CCA$eig[2], 3)*100, "% CCA2")) +
grids(linetype = "dashed") +
theme(legend.position = "none",
panel.background = element_rect(fill = "white", colour = "grey50")) +
labs(title = "A. Sites")
cca_sampl
p.cca <- ggpubr::ggarrange(cca_sampl, cca_clust, legend = FALSE)
ggpubr::annotate_figure(p.cca,
top = text_grob("CCA - cryoconites",
color = "black",
face = "bold",
size = 12,
hjust = 2.5),
fig.lab.pos = "top.left")
0.10.3 clusters
library(ggrepel)
tx.all <- ps.anc.f %>%
tax_table() %>%
as.data.frame() %>%
mutate(Group = net3$colors) %>%
add_column(abnd = taxa_sums(ps.rel)[taxa_names(ps.anc.f)])
# %>%
# filter(Group == "")
otus.ps.vegan <- veganifyOTU(ps.anc.f)
rownames(otus.ps.vegan) <- ps.anc.f@sam_data$Group
metadata <- as(sample_data(ps.anc.f), "data.frame")
cca_w_varstab_asv <- vegan::cca(otus.ps.vegan ~ pH + N.NH4 + TOC + P2O5 + Cu + Pb + Zn + Cd + Ni + N.NO3, data=metadata)
tx.ids <- tx.all %>%
rownames_to_column("ID") %>%
mutate(Score = "sites")
cca.meta <- cca_w_varstab_asv$CCA$biplot %>%
as.data.frame() %>%
select(c(CCA1, CCA2)) %>%
rownames_to_column("ID") %>%
mutate(Score = "biplot")
fdat_amazing <- as.data.frame(cca_w_varstab_asv$CCA$v) %>%
select(c(CCA1, CCA2)) %>%
rownames_to_column("ID") %>%
right_join(tx.ids, by = join_by(ID)) %>%
mutate(label_size = 0) %>%
bind_rows(cca.meta)
ggplot(fdat_amazing) +
geom_point(data = fdat_amazing %>% dplyr::filter(Score == "sites"),
mapping = aes(x=CCA1,
y=CCA2,
color=Group,
size=abnd),
alpha=0.5) +
geom_segment(data = fdat_amazing %>% dplyr::filter(Score == "biplot"),
aes(x = 0,
xend = CCA1,
y = 0,
yend = CCA2),
alpha=0.8,
color = "red",
arrow = arrow(angle = 3)) +
geom_text_repel(data = fdat_amazing %>% dplyr::filter(Score == "biplot"),
aes(x=CCA1,
y=CCA2,
label= ID),
size=6) +
xlab(paste0(round(cca$CCA$eig[1], 3)*100, "% CCA1")) +
ylab(paste0(round(cca$CCA$eig[2], 3)*100, "% CCA2")) +
grids(linetype = "dashed") +
geom_vline(xintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
geom_hline(yintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
theme(legend.position = "top",
panel.background = element_rect(fill = "white", colour = "grey50")) +
scale_colour_identity()
phylotypes(ASVs level) - without any specific filtration/normalisation at all
only top 50 ASVs based on relative abundance
library(ggrepel)
physeq <- ps.cr
ps.varstab <- physeq
veganifyOTU <- function(physeq){
require(phyloseq)
if(taxa_are_rows(physeq)){physeq <- t(physeq)}
return(as(otu_table(physeq), "matrix"))
}
otus.ps.vegan <- veganifyOTU(physeq)
metadata <- as(sample_data(physeq), "data.frame")
cca_w_varstab_asv <- vegan::cca(otus.ps.vegan ~ pH + N.NH4 + TOC + P2O5 + Cu + Pb + Zn + Cd + Ni + N.NO3, data=metadata)
wa <- as.data.frame(cca_w_varstab_asv$CCA$v) %>%
rownames_to_column("ID")
taxa.pruned <- as.data.frame(ps.varstab@tax_table@.Data) %>%
rownames_to_column("ID")
taxa.pruned$taxa <- ifelse(is.na(taxa.pruned$Genus),
ifelse(is.na(taxa.pruned$Family),
ifelse(is.na(taxa.pruned$Order),
ifelse(is.na(taxa.pruned$Class), taxa.pruned$Phylum, taxa.pruned$Class) , taxa.pruned$Order), taxa.pruned$Family), taxa.pruned$Genus)
taxa.pruned[taxa.pruned == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
taxa.pruned[taxa.pruned == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Pararhizobium"
taxa.pruned$taxa2 <- ifelse(is.na(taxa.pruned$Species),
with(taxa.pruned, paste0(taxa)),
with(taxa.pruned, paste0(taxa, " ", Species )))
taxa.pruned$phylum <- ifelse(taxa.pruned$Phylum == "Proteobacteria", with(taxa.pruned, paste0(taxa.pruned$Class)), with(taxa.pruned, paste0(taxa.pruned$Phylum)))
# taxa.pruned$Label <- paste0(taxa.pruned$ID, "_", taxa.pruned$taxa2)
taxa.pruned$Label <- taxa.pruned$taxa2
wa <- full_join(wa, taxa.pruned, by="ID") %>%
mutate(Score = "sites")
#For the top 50 most abundant taxa, skip if use des res
sw <- summarize_all(as.data.frame(ps.varstab@otu_table), sum)
head_asv <- as_tibble(t(sw), rownames = "ID") %>%
arrange(desc(V1)) %>% top_n(n = 50) %>%
pull(ID)
wa <- filter(wa, ID %in% head_asv)
biplot <- as.data.frame(cca_w_varstab_asv$CCA$biplot)
biplot <- rownames_to_column(biplot, "Label") %>%
add_column(Score = rep("biplot", length(rownames(biplot))))
fdat_amazing <- plyr::rbind.fill(biplot, wa) %>%
mutate(Label = as.factor(Label))
fdat_amazing <- fdat_amazing %>%
mutate(label_size = ifelse(Score == "biplot", 8, 5))
p.cca.species <- ggplot(fdat_amazing %>% filter(Score %in% c("sites","biplot"))) +
geom_point(data = fdat_amazing %>%
dplyr::filter(Score == "sites"), mapping = aes(x=CCA1, y=CCA2)) +
geom_segment(data = fdat_amazing %>%
dplyr::filter(Score == "biplot"), aes(x = 0, xend = CCA1, y = 0, yend = CCA2), alpha=0.8, color = "red", arrow = arrow(angle = 3)) +
geom_text_repel(aes(x=CCA1, y=CCA2, label = Label, colour = phylum),
force=3,
force_pull=4,
max.iter = 100000,
max.time = 3,
max.overlaps = 20,
size=fdat_amazing$label_size) +
xlab(paste0(round(cca$CCA$eig[1], 3)*100, "% CCA1")) +
ylab(paste0(round(cca$CCA$eig[2], 3)*100, "% CCA2")) +
grids(linetype = "dashed") +
geom_vline(xintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
geom_hline(yintercept = 0, size = 0.75, color = "#737373", alpha=0.5) +
theme(legend.position = "top", panel.background = element_rect(fill = "white", colour = "grey50"))
p.cca.species
0.11 treeplot
first - ASV level
secons - genus level
tx.cn <- ps.cn %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Antarctic", length(taxa_names(ps.cn)))) %>%
add_column(abnd = taxa_sums(ps.cn)[taxa_names(ps.cn)])
tx.me <- ps.me %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Arctic", length(taxa_names(ps.me)))) %>%
add_column(abnd = taxa_sums(ps.me)[taxa_names(ps.me)])
tx.ph <- ps.ph %>%
tax_table() %>%
as.data.frame() %>%
add_column(Group = rep("Caucasus", length(taxa_names(ps.ph)))) %>%
add_column(abnd = taxa_sums(ps.ph)[taxa_names(ps.ph)])
tx <- rbind(tx.cn, tx.me, tx.ph)
tx %>%
rownames_to_column("ID") %>%
mutate(Phylum = case_when(Phylum %in% "Pseudomonadota" ~ Class,
!Phylum %in% "Pseudomonadota" ~ Phylum) %>% as.factor()) %>%
select(c("ID", "Genus", "Phylum", "Group", "abnd")) %>%
group_by(ID, Genus, Phylum, Group ) %>%
mutate(Genus = str_replace_na(Genus, replacement = ' ') %>% as.factor()) %>%
summarise(area = sum(abnd)) %>%
mutate_if(is.character, as.factor) %>%
drop_na() %>%
mutate(Group = forcats::fct_relevel(Group, c("Antarctic","Arctic", "Caucasus"))) %>%
ggplot(aes(area = area,
subgroup = Phylum,
label = Genus,
fill = Group)) +
treemapify::geom_treemap() +
treemapify::geom_treemap_subgroup_border(color = "white") +
treemapify::geom_treemap_subgroup_text(place = "centre",
grow = T,
alpha = 0.7,
colour = "black",
fontface = "italic",
min.size = 2) +
treemapify::geom_treemap_text(colour = "white",
place = "topleft",
grow = T,
reflow = T,
layout = 'squarified',
min.size = 4) +
scale_colour_viridis_d(option = "magma",
aesthetics = "fill",
begin = 0.4,
end = 0.8) +
theme(legend.position="bottom") 
tx %>%
mutate(Phylum = case_when(Phylum %in% "Pseudomonadota" ~ Class,
!Phylum %in% "Pseudomonadota" ~ Phylum) %>% as.factor()) %>%
select(c("Genus", "Phylum", "Group", "abnd")) %>%
group_by(Genus, Phylum, Group ) %>%
mutate(Genus = str_replace_na(Genus, replacement = ' ') %>% as.factor()) %>%
summarise(area = sum(abnd)) %>%
mutate_if(is.character, as.factor) %>%
drop_na() %>%
mutate(Group = forcats::fct_relevel(Group, c("Antarctic","Arctic", "Caucasus"))) %>%
ggplot(aes(area = area,
subgroup = Phylum,
label = Genus,
fill = Group)) +
treemapify::geom_treemap() +
treemapify::geom_treemap_subgroup_border(color = "white") +
treemapify::geom_treemap_subgroup_text(place = "centre",
grow = T,
alpha = 0.7,
colour = "black",
fontface = "italic",
min.size = 2) +
treemapify::geom_treemap_text(colour = "white",
place = "topleft",
grow = T,
reflow = T,
layout = 'squarified',
min.size = 4) +
scale_colour_viridis_d(option = "magma",
aesthetics = "fill",
begin = 0.4,
end = 0.8) +
theme(legend.position="bottom") 
same, but in table format
tx %>%
rownames_to_column("ID") %>%
mutate(Phylum = case_when(Phylum %in% "Pseudomonadota" ~ Class,
!Phylum %in% "Pseudomonadota" ~ Phylum) %>% as.factor()) %>%
select(c("ID", "Genus", "Species", "Phylum", "Group", "abnd")) %>%
group_by(ID, Genus, Species, Phylum, Group ) %>%
mutate(Genus = str_replace_na(Genus, replacement = ' ') %>% as.factor()) %>%
summarise(area = sum(abnd)) %>%
mutate_if(is.character, as.factor) %>%
drop_na() %>%
mutate(area = round(area, 2)) %>%
pivot_wider(names_from = Group, values_from = area) %>%
DT::datatable(caption = "Groups based on WGCNA clusters")0.12 Packages info
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Moscow
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.3 KneeArrower_1.0.0 ggtree_3.6.2
## [4] heatmaply_1.4.2 viridis_0.6.2 viridisLite_0.4.1
## [7] plotly_4.10.1 ampvis2_2.7.28 ggpubr_0.6.0
## [10] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
## [13] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
## [16] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2
## [19] tidyverse_2.0.0 phyloseq_1.42.0 WGCNA_1.72-1
## [22] fastcluster_1.2.3 dynamicTreeCut_1.63-1
##
## loaded via a namespace (and not attached):
## [1] IRanges_2.32.0 gld_2.6.6
## [3] nnet_7.3-18 ggfittext_0.9.1
## [5] DT_0.27 Biostrings_2.66.0
## [7] TH.data_1.1-1 vctrs_0.6.2
## [9] energy_1.7-11 digest_0.6.31
## [11] png_0.1-8 proxy_0.4-27
## [13] Exact_3.2 registry_0.5-1
## [15] deldir_1.0-6 permute_0.9-7
## [17] MASS_7.3-58.2 reshape2_1.4.4
## [19] foreach_1.5.2 BiocGenerics_0.44.0
## [21] withr_2.5.0 ggfun_0.0.9
## [23] xfun_0.37 ellipsis_0.3.2
## [25] survival_3.5-3 doRNG_1.8.6
## [27] memoise_2.0.1 ggbeeswarm_0.7.1
## [29] emmeans_1.8.4-1 gmp_0.7-1
## [31] systemfonts_1.0.4 tidytree_0.4.2
## [33] zoo_1.8-11 DEoptimR_1.0-11
## [35] Formula_1.2-5 KEGGREST_1.38.0
## [37] httr_1.4.6 rstatix_0.7.2
## [39] rhdf5filters_1.10.0 rhdf5_2.42.0
## [41] rstudioapi_0.14 ggVennDiagram_1.2.2
## [43] units_0.8-1 generics_0.1.3
## [45] base64enc_0.1-3 S4Vectors_0.36.2
## [47] zlibbioc_1.44.0 ScaledMatrix_1.6.0
## [49] polyclip_1.10-4 ca_0.71.1
## [51] GenomeInfoDbData_1.2.9 xtable_1.8-4
## [53] ade4_1.7-22 doParallel_1.0.17
## [55] evaluate_0.20 preprocessCore_1.60.2
## [57] hms_1.1.3 GenomicRanges_1.50.2
## [59] irlba_2.3.5.1 colorspace_2.1-0
## [61] readxl_1.4.2 magrittr_2.0.3
## [63] lattice_0.20-45 robustbase_0.95-0
## [65] DECIPHER_2.26.0 scuttle_1.8.4
## [67] cowplot_1.1.1 matrixStats_0.63.0
## [69] class_7.3-21 Hmisc_4.8-0
## [71] pillar_1.9.0 nlme_3.1-162
## [73] iterators_1.0.14 decontam_1.18.0
## [75] compiler_4.3.0 beachmat_2.14.0
## [77] stringi_1.7.12 biomformat_1.26.0
## [79] treemapify_2.5.5 DescTools_0.99.48
## [81] TSP_1.2-2 sf_1.0-9
## [83] minqa_1.2.5 SummarizedExperiment_1.28.0
## [85] dendextend_1.16.0 plyr_1.8.8
## [87] crayon_1.5.2 abind_1.4-5
## [89] scater_1.26.1 gridGraphics_0.5-1
## [91] bit_4.0.5 mia_1.6.0
## [93] rootSolve_1.8.2.3 sandwich_3.0-2
## [95] codetools_0.2-19 multcomp_1.4-22
## [97] BiocSingular_1.14.0 crosstalk_1.2.0
## [99] bslib_0.4.2 e1071_1.7-13
## [101] lmom_2.9 multtest_2.54.0
## [103] MultiAssayExperiment_1.24.0 splines_4.3.0
## [105] Rcpp_1.0.10 sparseMatrixStats_1.10.0
## [107] cellranger_1.1.0 interp_1.1-3
## [109] knitr_1.42 egg_0.4.5
## [111] blob_1.2.3 utf8_1.2.3
## [113] lme4_1.1-31 checkmate_2.1.0
## [115] DelayedMatrixStats_1.20.0 Rdpack_2.4
## [117] expm_0.999-7 ggplotify_0.1.0
## [119] gsl_2.1-8 ggsignif_0.6.4
## [121] estimability_1.4.1 Matrix_1.5-3
## [123] tzdb_0.3.0 tweenr_2.0.2
## [125] pkgconfig_2.0.3 tools_4.3.0
## [127] cachem_1.0.7 rbibutils_2.2.13
## [129] RSQLite_2.3.0 DBI_1.1.3
## [131] numDeriv_2016.8-1.1 signal_0.7-7
## [133] impute_1.72.3 fastmap_1.1.1
## [135] rmarkdown_2.20 scales_1.2.1
## [137] grid_4.3.0 broom_1.0.4
## [139] sass_0.4.5 patchwork_1.1.2
## [141] coda_0.19-4 carData_3.0-5
## [143] rpart_4.1.19 farver_2.1.1
## [145] mgcv_1.8-42 yaml_2.3.7
## [147] latticeExtra_0.6-30 MatrixGenerics_1.10.0
## [149] foreign_0.8-84 bayesm_3.1-5
## [151] cli_3.6.0 stats4_4.3.0
## [153] webshot_0.5.4 lifecycle_1.0.3
## [155] Biobase_2.58.0 mvtnorm_1.1-3
## [157] backports_1.4.1 BiocParallel_1.32.5
## [159] timechange_0.2.0 gtable_0.3.1
## [161] ANCOMBC_2.0.2 parallel_4.3.0
## [163] ape_5.7 CVXR_1.0-11
## [165] jsonlite_1.8.4 seriation_1.4.1
## [167] bitops_1.0-7 bit64_4.0.5
## [169] assertthat_0.2.1 yulab.utils_0.0.6
## [171] vegan_2.6-4 BiocNeighbors_1.16.0
## [173] TreeSummarizedExperiment_2.6.0 jquerylib_0.1.4
## [175] highr_0.10 lazyeval_0.2.2
## [177] htmltools_0.5.4 GO.db_3.16.0
## [179] glue_1.6.2 XVector_0.38.0
## [181] RCurl_1.98-1.10 treeio_1.22.0
## [183] classInt_0.4-9 jpeg_0.1-10
## [185] gridExtra_2.3 boot_1.3-28.1
## [187] igraph_1.4.1 R6_2.5.1
## [189] SingleCellExperiment_1.20.0 Rmpfr_0.9-1
## [191] labeling_0.4.2 cluster_2.1.4
## [193] rngtools_1.5.2 Rhdf5lib_1.20.0
## [195] aplot_0.1.9 GenomeInfoDb_1.34.9
## [197] nloptr_2.0.3 compositions_2.0-6
## [199] DirichletMultinomial_1.40.0 DelayedArray_0.24.0
## [201] tidyselect_1.2.0 vipor_0.4.5
## [203] htmlTable_2.4.1 tensorA_0.36.2
## [205] ggforce_0.4.1 car_3.1-1
## [207] AnnotationDbi_1.60.0 KernSmooth_2.23-20
## [209] rsvd_1.0.5 munsell_0.5.0
## [211] data.table_1.14.8 htmlwidgets_1.6.1
## [213] RColorBrewer_1.1-3 rlang_1.1.1
## [215] lmerTest_3.1-3 fansi_1.0.4
## [217] beeswarm_0.4.0 RVenn_1.1.0