Skip to content
Snippets Groups Projects
Select Git revision
  • main default protected
1 result
Blame Permalink
.Rhistory 20.05 KiB 
top = set$top_n,
ignore_na = set$ignore_na)
set$topnASVs <- names(sort(taxa_sums(ps), decreasing = TRUE))[1:set$topASVs]
set$ps.topnASVs <- prune_taxa(set$topnASVs, ps.trans)
if (set$unify_colors | exists("highlight", envir = set) | set$fuse_ASVs) {
set$toptax <- union(phyloseq::tax_table(set$ps.topnTax)[,set$taxlvl],
phyloseq::tax_table(set$ps.topnASVs)[,set$taxlvl])
set$toptax <- sort(set$toptax)
set$taxlvlPalette <- viridis(length(set$toptax))
names(set$taxlvlPalette) <- set$toptax
if (exists("highlight", envir = set)) {
# It is possible to change the highlight color here by substituting
# 'sub_viridis$reds[4]' with a hexcode-string, e.g. '#ff7f7f"'
set$taxlvlPalette[set$highlight] <- sub_viridis$reds[4]
}
set$taxlvlPalette <- set$taxlvlPalette[sort(names(set$taxlvlPalette))]
set$my_scale_fill <- scale_fill_manual(values = set$taxlvlPalette,
na.value = "grey")
}else{
set$my_scale_fill <- my_scale_fill
}
# Plot
if (set$unify_colors | exists("highlight", envir = set) | set$fuse_ASVs) {
set$my_scale_fill <- scale_fill_manual(
values = set$taxlvlPalette[
sort(unique(phyloseq::tax_table(set$ps.topnTax)[,set$taxlvl]))],
na.value = "grey")
}
plots$topn_tax <- plot_bar(set$ps.topnTax,
x = set$x_axis_value,
fill = set$taxlvl,
title = paste0("Top ", set$top_n, " ", set$taxlvl)) +
facet_grid(paste0("~", set$panel_by), scales = "free", space = "free") +
set$my_scale_fill +
ylab(set$y_axis_label) +
xlab("Sample") +
theme(strip.background = element_blank(), strip.text = element_text(size = 16),
axis.title=element_text(size=16), legend.text = element_text(size=14))
if (set$fuse_ASVs) {
plots$topn_tax <- plots$topn_tax + geom_bar(
aes_string(color = set$taxlvl, fill = set$taxlvl),
stat = "identity", position = "stack") +
scale_color_manual(values = set$taxlvlPalette, na.value = NA)
}
if (set$unify_colors | exists("highlight", envir = set) | set$fuse_ASVs) {
set$my_scale_fill <- scale_fill_manual(
values = set$taxlvlPalette[
sort(unique(phyloseq::tax_table(set$ps.topnASVs)[,set$taxlvl]))],
na.value = "grey")
}
plots$topn_ASVs <- plot_bar(set$ps.topnASVs,
x = set$x_axis_value,
fill = set$taxlvl,
title = paste0("Top", set$topASVs, "_ASVs")) +
facet_grid(paste0("~", set$panel_by), scales = "free_x", space = "free") +
set$my_scale_fill +
ylab(set$y_axis_label) +
xlab("Sample") +
theme(strip.background = element_blank())
# save
save_plot(plots$topn_tax, plot_name = paste0("Top", set$top_n, "_", set$taxlvl),
filetype = tmp$out)
save_plot(plots$topn_ASVs, plot_name = paste0("Top", set$topASVs, "_ASVs"),
filetype = tmp$out)
# Clean up plot parameters
rm(list = ls(set), envir = set)
# Print to standard out
plots$topn_tax
plots$topn_ASVs
# CHANGE ME to the desired sample categories on the x-axis.
# Accepted values are the column headers in the descriptor file.
set$x_axis_value = "ndvi"
# CHANGE ME to the taxonomic level of interest (color coding). Accepted values
# are the column headers in your descriptor file.
set$taxlvl = "Genus"
# CHANGE ME to change the number of Top n taxa to be plotted at
# taxlvl.
set$top_n = 10
# Can be changed to include (FALSE) or exclude (TRUE) NA values in the barplot
set$ignore_na = FALSE
# CHANGE ME to an entry at the chosen taxonomic level you want to highlight.
# Comment out to not highlight anything.
# set$highlight =
##### Optional settings (sensible defaults) #####
# Can be changed to change the width (in cm) of the saved plot.
set$wp = 20
# Can be changed to change the height (in cm) of the saved plot.
set$hp = 13
# Can be changed to change the resolution (in dpi) of the saved plot.
set$res = 300
# Can be changed to change the y-axis label
set$y_axis_label = "Relative abundance"
# Can be changed to change the x-axis label
set$x_axis_label = "Sample"
############### NO NEED FOR CHANGES BELOW ###############
# Estimate Alpha-diversity (Shannon)
set$alpha_div_ps_trans <- estimate_richness(ps.trans, measures = "Shannon") %>%
as_tibble(rownames = "Sample")
# Make physeq objects of top n taxa and top n ASVs
set$ps.topnTax <- cuphyr::abundant_tax_physeq(ps.trans, lvl = set$taxlvl,
top = set$top_n,
ignore_na = set$ignore_na)
# Plot
set$my_scale_fill <- my_scale_fill
set$topntax_tbl <- psmelt(set$ps.topnTax) %>%
as_tibble() %>%
left_join(set$alpha_div_ps_trans, by = "Sample") %>%
select(Genus, Alias, ndvi, ndvi_july, Abundance, Shannon) %>%
filter(Abundance > 0) %>%
group_by(Genus, Alias, ndvi, ndvi_july, Shannon) %>%
summarise(Abundance = sum(Abundance)) %>%
arrange(ndvi) %>%
mutate(ndvi_rank = c(1:length(ndvi)))
plots$topn_tax_custom <- ggplot(set$topntax_tbl, aes(x = fct_reorder(Alias, ndvi),
y = Abundance,
fill = Genus)) +
#title = paste0("Top ", set$top_n, " ", set$taxlvl))) +
geom_col(color = "black") +
set$my_scale_fill +
ylab(set$y_axis_label) +
xlab(set$x_axis_label) +
theme(strip.background = element_blank(),
#strip.text = element_text(size = 12),
#axis.title=element_text(size=12),
#legend.text = element_text(size=12),
legend.position = "bottom")
plots$ndvi_dot_plot <- ggplot(set$topntax_tbl, aes(fct_reorder(Alias, ndvi),
y = ndvi)) +
geom_point() +
theme(strip.background = element_blank(),
# strip.text = element_text(size = 12),
# axis.title=element_text(size=12),
# legend.text = element_text(size=12),
axis.title.x=element_blank()) +
ylab("NDVI")
plots$shannon_dot_plot <- ggplot(set$topntax_tbl,
aes(fct_reorder(Alias, ndvi),
y = Shannon)) +
geom_point() +
theme(strip.background = element_blank(),
# strip.text = element_text(size = 12),
# axis.title=element_text(size=12),
# legend.text = element_text(size=12),
axis.title.x=element_blank()) +
ylab("Shannon")
plots$combo_topn_custom <- ggarrange(plots$ndvi_dot_plot,
plots$shannon_dot_plot,
plots$topn_tax_custom,
nrow = 3,
heights = c(1, 1, 3),
align = "v")
plots$combo_topn_custom
save_plot(plots$combo_topn_custom, plot_name = paste0("Customized_NDVI_Shannon_plot"),
filetype = tmp$out)
set$my_formula <- y ~ x
set$plot_title <- "Nematodes"
topntax_data <- set$topntax_tbl %>%
mutate(Taxa = 'nematodes') %>%
ungroup() %>%
select(Alias, ndvi, Shannon, Taxa) %>%
distinct()
write.csv(topntax_data, file = "../topntax_all_taxa/topntax_data_nematodes.csv")
# Excel file with field data
morphological_id <- file.path("2020-patch-larvik-a-counts.csv")
# Reduce abundance per sample and genus data to genera of interest
genera_of_interest_mol_morph <- c("Globodera", "Meloidogyne", "Pratylenchus")
genus_abundance_tbl_per_sample_mm <- genus_abundance_tbl_per_sample %>%
filter(Genus %in% genera_of_interest_mol_morph)
# read in and parse morphological count data
morph_data <- read_delim(morphological_id) %>%
mutate(genus = str_replace(genus, "cyst", "Globodera"),
genus = str_replace(genus, "mel", "Meloidogyne"),
genus = str_replace(genus, "prat", "Pratylenchus")) %>%
dplyr::rename(Sample = rute, Genus = genus)
# Combine morphological and metabarcoding data
set$data_mol_morph <- genus_abundance_tbl_per_sample_mm %>%
left_join(morph_data)
### Plotting
# New facet label names for Genus variable
genus_names <- list(
'Globodera'=expression(paste(italic("Globodera "))), #"spp.")),
'Meloidogyne'=expression(paste(italic("Meloidogyne "))), #"spp.")),
'Pratylenchus'=expression(paste(italic("Pratylenchus ")))#, "spp."))
)
genus_labeller <- function(variable,value){
return(genus_names[value])
}
signif.labs <- c("pval > 0.05", "pval <= 0.05")
names(signif.labs) <- c("p value > 0.05", "p value <= 0.05")
# Subset data for modeling
set$glob_data <- subset(set$data_mol_morph, Genus == "Globodera")
set$mel_data <- subset(set$data_mol_morph, Genus == "Meloidogyne")
set$prat_data <- subset(set$data_mol_morph, Genus == "Pratylenchus")
# Quasi-binomial model for Globodera
model_glob = glm(
Genus_Sample_Abundance ~ total_count,
data = set$glob_data,
family = quasibinomial
)
# Quasi-binomial model for Meloidogyne
model_mel = glm(
Genus_Sample_Abundance ~ total_count,
data = set$mel_data,
family = quasibinomial
)
# Quasi-binomial model for Pratylenchus
model_prat = glm(
Genus_Sample_Abundance ~ total_count,
data = set$prat_data,
family = quasibinomial
)
pval_glm <- function(glm){
coef(summary(glm))[,4][2]
}
set$data_mol_morph <- set$data_mol_morph %>%
mutate(pval = ifelse(Genus == "Globodera", pval_glm(model_glob), NA)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus", pval_glm(model_prat), pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne", pval_glm(model_mel), pval)) %>%
mutate(signif = ifelse(pval <= 0.05, "pval <= 0.05", "pval > 0.05"))
# Other plotting variables
# Can be changed to change the width (in cm) of the saved plot.
set$wp = 20
# Can be changed to change the height (in cm) of the saved plot.
set$hp = 13
# Can be changed to change the resolution (in dpi) of the saved plot.
set$res = 300
set$my_formula <- y ~ x
set$plot_title <- paste("Molecular analysis vs",
"morpological analysis", sep = "\n")
# Plotting quasibinomial model
plots$mol_morph_quasi <- ggplot(set$data_mol_morph,
aes(x=total_count,
y = Genus_Sample_Abundance)) +
geom_point(alpha = 0.7) +
#ggtitle(set$plot_title) +
theme(plot.title = element_text(hjust = 0.5, size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
#axis.line = element_line(colour = "dark grey"),
strip.background = element_blank(),
strip.text = element_text(face = "italic", size = 16),
axis.text = element_text(),
text = element_text()) +
xlab("Nematodes / 250 ml soil") +
ylab("Relative abundance") +
geom_smooth(method = "glm", formula = set$my_formula, se = F,
method.args = list(family = quasibinomial),
aes(color = signif, linetype = signif)) +
scale_color_manual(values=c('black', 'lightgrey'),
name= 'Significance',
labels = names(signif.labs)) +
scale_linetype_discrete(name= 'Significance',
labels = names(signif.labs)) +
facet_wrap(~Genus, scales = "free",
#labeller = labeller(Genus = genus_labs)) +
labeller = genus_labeller)
#Print and save plot
plots$mol_morph_quasi
save_plot(plots$mol_morph_quasi, plot_name = paste0("Molecular_and_Morphological_Plot"),
filetype = ".pdf")
# Generate table overview of morphological counts for manuscript
morph_table_manuscript <- pivot_wider(set$data_mol_morph %>%
select(-Genus_Sample_Abundance,
-ndvi),
names_from = Genus,
values_from = total_count) %>%
arrange(Alias)
write.xlsx(morph_table_manuscript, file = "t1_morph_table_manuscript.xlsx")
# Create properly formatted tibble with columns Sample, ndvi_01 (ndvi translated to (0, 1) interval) and one column for each genus
# containing the sample abundances for that genus.
ldf <- data.frame(genus_abundance_tbl_per_sample %>% pivot_wider(id_cols = c('Sample', 'ndvi'), names_from = 'Genus', values_from = 'Genus_Sample_Abundance'))
ldf_genus_data <- data.frame(ldf) %>% select(!c('Sample', 'ndvi'))
colnames(ldf_genus_data) <- gsub(' ', '.', colnames(ldf_genus_data))
ldf <- cbind.data.frame(
Sample = ldf$Sample,
ndvi_01 = (ldf$ndvi + 1) / 2.0
)
ldf <- tibble(cbind(ldf, ldf_genus_data))
n_samples_by_genus <- data.frame(ldf_genus_data > 0) %>% mutate_if(is.logical, as.numeric) %>% colSums() %>% sort(decreasing = TRUE)
keep_n <- 100 # Maximum number of genera to include in the analysis
top_n_occurence_genuses <- names(n_samples_by_genus[1:keep_n])
top_n_occurence_genuses <- top_n_occurence_genuses[!is.na(top_n_occurence_genuses)]
l_genus_ldf <- ldf %>% select(all_of(top_n_occurence_genuses))
l_genus_ldf_transposed <- data.frame(t(l_genus_ldf))
l_meta_ldf <- ldf %>% select('ndvi_01')
l_model <- linda(
l_genus_ldf_transposed,
l_meta_ldf,
formula = '~ ndvi_01',
feature.dat.type = 'proportion',
is.winsor = FALSE,
alpha = 0.05
)
# Print model info
l_model
# Show effect size and significance plots
linda.plot(
l_model,
variables.plot = c('ndvi_01'),
alpha = 0.05,
lfc.cut = 1,
legend = TRUE
)
# Write supplementary table for manuscript
l_model_df <- as.data.frame(l_model$output)
write.xlsx(l_model_df, file = "supplementary_table_ndvi_regression_nematodes.xlsx", rowNames = TRUE, colnames = TRUE)
# Prat
set$prat_data$ndvi_temp <- (set$prat_data$ndvi+1)/2
model_prat_count_ndvi = glm(
ndvi_temp ~ total_count,
data = set$prat_data,
family = quasibinomial
)
summary(model_prat_count_ndvi)
# Glob
set$glob_data$ndvi_temp <- (set$glob_data$ndvi+1)/2
model_glob_count_ndvi = glm(
ndvi_temp ~ total_count,
data = set$glob_data,
family = quasibinomial
)
summary(model_glob_count_ndvi)
# Mel
set$mel_data$ndvi_temp <- (set$mel_data$ndvi+1)/2
model_mel_count_ndvi = glm(
ndvi_temp ~ total_count,
data = set$mel_data,
family = quasibinomial
)
summary(model_mel_count_ndvi)
data_mol_morph_long_temp <- pivot_longer(set$data_mol_morph,
cols = c(Genus_Sample_Abundance,
total_count),
names_to = "type",
values_to = "value")
data_mol_morph_long_temp$ndvi <- (data_mol_morph_long_temp$ndvi+1)/2
pval_glm <- function(glm){
coef(summary(glm))[,4][2]
}
set$data_mol_morph_long <- data_mol_morph_long_temp %>%
mutate(pval = ifelse(Genus == "Globodera", pval_glm(model_glob), NA)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus", pval_glm(model_prat), pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne", pval_glm(model_mel), pval)) %>%
mutate(signif = ifelse(pval <= 0.05, "pval <= 0.05", "pval > 0.05"))
# Count data vs NDVI
plots$morph_count_ndvi <- ggplot(set$data_mol_morph_long,
aes(x = value,
y = ndvi)) +
geom_point(alpha = 0.7) +
#ggtitle(set$plot_title) +
theme(plot.title = element_text(hjust = 0.5, size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "dark grey"),
strip.background = element_blank(),
strip.text.x = element_text(face = "italic")#,
#axis.text = element_text(size = 12),
#text = element_text(size = 14)
) +
#xlab("Nematodes / 250 ml soil") +
#ylab("NDVI") +
geom_smooth(method = "glm", formula = set$my_formula, se = F, color = "dark grey",
method.args = list(family = quasibinomial)) +
facet_wrap(~Genus, scales = "free_x")
plots$morph_count_ndvi
# Print and save plot
plots$morph_count_ndvi
save_plot(plots$morph_count_ndvi, plot_name = paste0("morph_count_ndvi"),
filetype = tmp$out)
set$data_mol_morph_long <- data_mol_morph_long_temp %>%
mutate(pval = ifelse(Genus == "Globodera", pval_glm(model_glob_count_ndvi), NA)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus", pval_glm(model_prat_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne", pval_glm(model_mel_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Globodera" & type == "Genus_Sample_Abundance", 0.5574120890, pval)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus" & type == "Genus_Sample_Abundance", 0.5177548997, pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne" & type == "Genus_Sample_Abundance", 0.0001920515, pval)) %>%
mutate(signif = ifelse(pval <= 0.05, "pval <= 0.05", "pval > 0.05"))
# Morph data
set$data_morph <- set$data_mol_morph_long[set$data_mol_morph_long$type == 'total_count',]
set$data_morph <- set$data_morph %>%
rename("value" = "Number_of_nematodes")
set$data_mol_morph_long
set$data_mol_morph_long <- data_mol_morph_long_temp %>%
mutate(pval = ifelse(Genus == "Globodera", pval_glm(model_glob_count_ndvi), NA)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus", pval_glm(model_prat_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne", pval_glm(model_mel_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Globodera" & type == "Genus_Sample_Abundance", 0.5574120890, pval)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus" & type == "Genus_Sample_Abundance", 0.5177548997, pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne" & type == "Genus_Sample_Abundance", 0.0001920515, pval)) %>%
mutate(signif = ifelse(pval <= 0.05, "pval <= 0.05", "pval > 0.05"))
# Morph data
set$data_morph <- set$data_mol_morph_long[set$data_mol_morph_long$type == 'total_count',]
set$data_morph <- set$data_morph %>%
dplyr::rename("Number_of_nematodes" = value)
# Relative abundance data
set$data_relabu <- set$data_mol_morph_long[set$data_mol_morph_long$type == 'Genus_Sample_Abundance',]
set$data_relabu <- set$data_relabu %>%
rename("value" = "Relative_abundance")
set$data_mol_morph_long <- data_mol_morph_long_temp %>%
mutate(pval = ifelse(Genus == "Globodera", pval_glm(model_glob_count_ndvi), NA)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus", pval_glm(model_prat_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne", pval_glm(model_mel_count_ndvi), pval)) %>%
mutate(pval = ifelse(Genus == "Globodera" & type == "Genus_Sample_Abundance", 0.5574120890, pval)) %>%
mutate(pval = ifelse(Genus == "Pratylenchus" & type == "Genus_Sample_Abundance", 0.5177548997, pval)) %>%
mutate(pval = ifelse(Genus == "Meloidogyne" & type == "Genus_Sample_Abundance", 0.0001920515, pval)) %>%
mutate(signif = ifelse(pval <= 0.05, "pval <= 0.05", "pval > 0.05"))
# Morph data
set$data_morph <- set$data_mol_morph_long[set$data_mol_morph_long$type == 'total_count',]
set$data_morph <- set$data_morph %>%
dplyr::rename("Number_of_nematodes" = value)
# Relative abundance data
set$data_relabu <- set$data_mol_morph_long[set$data_mol_morph_long$type == 'Genus_Sample_Abundance',]
set$data_relabu <- set$data_relabu %>%
dplyr::rename("Relative_abundance" = value)
signif.labs <- c("pval > 0.05", "pval <= 0.05")
names(signif.labs) <- c("p value > 0.05", "p value <= 0.05")
breaks_morph <- function(x) { if (max(x) < 200) seq (0, 100, 100) else seq(0, 600, 50)}
breaks_mb <- function(x) { if (max(x) < 0.1) seq(0, 0.002, 3) else seq(0, 3, 1) }
breaks_fun_morph <- function(x) {
if (max(x) > 500) {
seq(0, 800, 200)
} else if (max(x) > 150) {
seq(0, 200, length.out = 5)
} else {
seq(0, 100, length.out = 5)
}}
# MORPHOLOGIAL AND NDVI PLOT
# Relative abundance and NDVI plotting
# Axis ticks MB
breaks_fun_mb <- function(x) {
if (max(x) > 0.8) {
seq(0, 0.9, length.out = 4)
} else if (max(x) > 0.2) {
seq(0, 0.6, 0.2)
} else {
seq(0, 0.003, length.out = 4)
}}
p1 <- ggplot(set$data_relabu,
aes(x = Relative_abundance,
y = ndvi)) +
geom_point(alpha = 0.7, color = "black") +
ggtitle("a)") +
theme(plot.title = element_text(hjust = 0, size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
strip.background = element_blank(),
axis.title = element_text(size = 9),
strip.text.x = element_text(face = "italic"),
plot.caption = element_text(size = 7, hjust = 0, face = "italic",
margin = margin(t = 0.2, unit = "cm")),
plot.caption.position = "plot",
plot.subtitle = element_text(size = 10, hjust = 0, vjust = 0.5,
margin = margin(b = 0.8, unit = "cm"))) +
xlab("Relative abundance") +
ylab("NDVI") +
geom_smooth(method = "glm", formula = set$my_formula, se = FALSE,
method.args = list(family = quasibinomial),
aes(color = signif, linetype = signif)) +
scale_color_manual(values=c('black', 'lightgrey'),
name= 'Significance',
labels = names(signif.labs)) +
scale_linetype_discrete(name= 'Significance',
labels = names(signif.labs)) +
facet_wrap(~Genus, scales = "free_x") +
scale_x_continuous(breaks = breaks_fun_mb, limits = c(0, NA))
p1
# Plotting NDVI vs metabarcoding
p2 <- ggplot(set$data_morph,
aes(x = Number_of_nematodes,
y = ndvi)) +
geom_point(alpha = 0.7, color = "black") +
ggtitle("b)") +
theme(plot.title = element_text(hjust = 0, size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
strip.background = element_blank(),
axis.title = element_text(size = 9),
panel.spacing=unit(2, "lines"),
strip.text.x = element_text(face = "italic"),
plot.caption = element_text(size = 7, hjust = 0, face = "italic",
margin = margin(t = 0.2, unit = "cm")),
plot.caption.position = "plot",
plot.subtitle = element_text(size = 10, hjust = 0, vjust = 0.5,
margin = margin(b = 0.8, unit = "cm"))) +
xlab("Nematodes / 250 ml soil") +
ylab("NDVI") +
geom_smooth(method = "glm", formula = set$my_formula, se = FALSE,
method.args = list(family = quasibinomial),
aes(color = signif, linetype = signif)) +
scale_color_manual(values=c('black', 'lightgrey'),
name= 'Significance',
labels = names(signif.labs)) +
scale_linetype_discrete(name= 'Significance',
labels = names(signif.labs)) +
facet_wrap(~Genus, scales = "free_x") +
scale_x_continuous(breaks = breaks_fun_morph, limits = c(0, NA))
p2
# Combining the two plots
plots$combo_morph_metab_ndvi <- plot_grid(p1, p2,
ncol = 1,
align = "hv")
plots$combo_morph_metab_ndvi
# Print and save plot
save_plot(plots$combo_morph_metab_ndvi, plot_name = paste0("morph_relabu_ndvi"),
filetype = tmp$out)
found <- c(40, 24, 16, 37, 10, 27, 74, 100, 100, 78, 32, 28)
not_found <- c(24,12, 7, 7, 24, 52)
wilcox.test(found, not_found)