CoMet.Rmd

---
title: "Comparative Metataxonomy: example"
author: "Roberto Siani"
date: "23/11/2021"
---

PROJECT NOTES:
insert short intro to project here

# SETUP

Set up working environment: packages and default graphic themes.

```{r, cache = TRUE, echo = FALSE, include = FALSE}

# pacman to install and load libraries

if (!require("pacman")) {
  install.packages(
    "pacman",
    verbose = F
  )
}

# and BiocManager for Bioconductor libraries

if (!require("BiocManager")) {
  install.packages(
    "BiocManager",
    verbose = F
  )
}

# install/load GitHub libraries

pacman::p_load_gh(
  "jbisanz/qiime2R",
  "benjjneb/decontam",
  "mikemc/speedyseq",
  "adw96/breakaway",
  "adw96/DivNet"
)

# install/load CRAN/BioConductor libraries

pacman::p_load(
  vegan,
  ape,
  microbiome,
  hrbrthemes,
  phangorn,
  DECIPHER,
  MicrobiotaProcess,
  hues,
  ggtree,
  ggridges,
  ANCOMBC
)

# load sensible themes and palettes (from a separate script!)

source("~/Desktop/library/scripts/helpeR.R")

# define functions

select = dplyr::select
transform = microbiome::transform

```

# PRE-PROCESS

## SELECT PATHS

```{r}

metadata_table = "map.txt"

taxonomic_assignment = "taxonomy-dada2.qza"

counts_table = "table-dada2.qza"

control_samples = "rob22"

```


## PRE-PROCESS TO PHYLOSEQ OBJECT

```{r}

# import metadata

metaTab =
  read_tsv(metadata_table) %>%
  mutate_if(is_character, as.factor)

# import and clean taxonomy

raw_taxTab =
  qiime2R::read_qza(taxonomic_assignment)$data %>%
  separate(Taxon,
           sep = ";",
           c("Domain",
             "Phylum",
             "Class",
             "Order",
             "Family",
             "Genus",
             "Species"),
           fill = "right",
           extra = "drop") %>%
  map_df(~gsub(
    pattern = "D_[0-9]__",
    replacement = "",
    .x)) %>%
  map_df(~ gsub(
    pattern = "metagenome|uncultured|unidentified|Unknown| |  ",
    replacement = NA,
    .x
  )) %>%
  mutate_if(is_character, str_trim) %>%
  mutate(
    Domain = ifelse(is.na(Domain),
                    "U. Domain",
                    Domain
    ),
    Phylum = coalesce(
      Phylum,
      ifelse(grepl("^U.", Domain),
        Domain,
        paste("U.", Domain)
      )
    ),
    Class = coalesce(
      Class,
      ifelse(grepl("^U.", Phylum),
        Phylum,
        paste("U.", Phylum)
      )
    ),
    Order = coalesce(
      Order,
      ifelse(grepl("^U.", Class),
        Class,
        paste("U.", Class)
      )
    ),
    Family = coalesce(
      Family,
      ifelse(grepl("^U.", Order),
        Order,
        paste("U.", Order)
      )
    ),
    Genus = coalesce(
      Genus,
      ifelse(grepl("^U.", Family),
        Family,
        paste("U.", Family)
      )
    ),
    Species = coalesce(
      Species,
      ifelse(grepl("^U.", Genus),
        Genus,
        paste("U.", Genus)
      )
    )
  ) %>%
  column_to_rownames("Feature.ID") %>%
  filter(Domain %in% "Bacteria" &
    !Order %in% "Chloroplast" &
    !Family %in% "Mitochondria") %>% 
  select(-Confidence)

# import counts table

raw_abuTab =
  qiime2R::read_qza(counts_table)$data %>%
  as.data.frame()

# reads per sample

summary((colSums(raw_abuTab)))

# clean counts table

abuTab =
  raw_abuTab %>%
  filter(!rownames(.) %in%
    isContaminant(
      as.matrix(.),
      neg = c(colnames(.) %in%
        control_samples),
      threshold = .1,
      normalize = T,
      detailed = F
    )) %>%
  select_if(!names(.) %in% control_samples) %>%
  filter(rowSums(.) > 0 &
    rownames(.) %in% rownames(raw_taxTab))

# clean reads per sample

summary((colSums(abuTab)))

# filtered reads

colSums(raw_abuTab[,-14]) - c(colSums(abuTab))

# final cleaning of taxonomy table

taxTab =
  raw_taxTab %>%
  filter(rownames(.) %in%
    rownames(abuTab)) %>%
  as.matrix()

# multi-alignment to predict phylogenetic tree
# this step is long and not always necessary

treeObj =
  phyDat(
    as(
      AlignSeqs(
        DNAStringSet(
          x = rownames(abuTab) %>%
            set_names(),
          use.names = T
        ),
        anchor = NA,
        verbose = F
      ),
      "matrix"
    ),
    type = "DNA"
  ) %>%
  pml(dist.ml(.) %>%
    NJ(),
  data = .
  ) %>%
  optim.pml(.)

# import into phyloseq-object

meCom =
  phyloseq(
    otu_table(abuTab,
      taxa_are_rows = T),
    tax_table(taxTab),
    sample_data(metaTab %>%
      column_to_rownames("#SampleID")))

# create directory for results

dir.create("OUT")

# save phyloseq-class object

save(meCom,
  file = paste("meCom",
    format(Sys.time(), "_%d%m%y_%H%M"),
    ".RData",
    sep = ""
  )
)

```

# INPUT

```{r, echo = FALSE, include = FALSE}

# load your meCom dataset

load("meCom_221121_1128.RData")

# check your summary

summarize_phyloseq(meCom)

# check if you have low depth samples to discard!

colSums(meCom@otu_table) %>% 
  as.data.frame() %>% 
  rownames_to_column("sample_names") %>% 
  left_join(meta(meCom) %>% 
              rownames_to_column("sample_names")) %>% 
  ggplot() +
  geom_bar(stat = "identity",
           position  = "dodge",
           aes(x = `.`,
               y = sample_names,
               fill = Bacteria)) +
  scale_fill_iwanthue()

# it looks like sample R20 should be discarded

meCom = 
  meCom %>% 
  subset_samples(Bacteria == "LjR87" |
                   Bacteria == "WT1" |
                   Bacteria == "LjR176") %>% 
  prune_taxa(taxa_sums(.) > 0, .)

# control again

summarize_phyloseq(meCom)

# define a color palette for data viz

palette_list =
  list(
    Bacteria = c("#006D77", "#E9C46A", "#C1666B"),
    Description = c("#ff506e", "#69005f")
  )

map(
  palette_list,
  ~ swatches::show_palette(.x)
)

# check that all levels are approximately equal

table(meta(meCom)$Bacteria, 
      meta(meCom)$Description)
```

# OUTPUT 

## MICROBIOME OVERVIEW

### General stats

```{r}

# for heavily left-skewed dataset, you might want to filter out some rare taxa

prevFilt =
  meCom %>%
  metagMisc::phyloseq_filter_prevalence(
    prev.trh = 3 / nsamples(meCom),
    abund.trh = 100,
    threshold_condition = "OR",
    abund.type = "total"
  )

ntaxa(meCom)

ntaxa(prevFilt)

# just some operations to get an overview of your metacommunity

sum_df =
  prevFilt %>%
  transform("clr") %>%
  speedyseq::psmelt() %>%
  group_by_at(vars(-Abundance, -BarcodeSequence, -LinkerPrimerSequence, -Bacteria, -Description, -Pot, -Sample)) %>%
  summarise(
    CLR_mean = mean(Abundance),
    CLR_SD = sd(Abundance)
  ) %>%
  ungroup() %>%
  mutate(
    TRA = taxa_sums(meCom %>%
      transform("compositional")),
    Prevalence =
      prevalence(meCom, count = T),
    Phylum = as.factor(Phylum)
  )

# scatter plot of abundances against standard deviation

plot_1 =
  ggplot(
    sum_df,
    aes(
      x = CLR_mean,
      y = CLR_SD,
      color = Phylum
    )
  ) +
  geom_point(aes(
    alpha = Prevalence,
    size = TRA
  ),
  size = 3,
  shape = 19
  ) +
  theme(legend.position = "none") +
  scale_color_iwanthue(
    hmin = 0,
    hmax = 360,
    cmin = 20,
    cmax = 50,
    lmin = 50,
    lmax = 100
  ) +
  facet_wrap(~Phylum) +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    strip.text.x = element_text(size = 12)
  ) +
  labs(
    x = "Abundance",
    y = "Dispersion"
  )

# density ridges of prevalence per phylum

plot_2 =
  ggplot(
    sum_df,
    aes(Prevalence,
      y = Phylum,
      fill = Phylum,
      color = Phylum
    )
  ) +
  geom_density_ridges(
    alpha = 1 / 4,
    size = 0.5
  ) +
  theme(legend.position = "none") +
  scale_color_iwanthue(
    hmin = 0,
    hmax = 360,
    cmin = 20,
    cmax = 50,
    lmin = 50,
    lmax = 100,
    aesthetics = c("fill", "color")
  ) +
  xlim(c(0, nsamples(meCom)))

cairo_pdf("OUT/overlook.pdf",
  width = 16,
  height = 16
)
plot_1 + plot_2 + plot_layout(widths = c(2, 1))
dev.off()

# tree, because we all like trees

plot_3 =
  ggtree(prevFilt %>%
    transform("compositional"),
  layout = "fan",
  open.angle = 10,
  size = 0.1
  ) +
  geom_rootedge(0.1,
    size = 0.1
  ) +
  geom_tippoint(
    aes(
      color = Phylum,
      size = Abundance,
      alpha = Abundance
    ),
    shape = 19
  ) +
  scale_color_iwanthue(
    hmin = 0,
    hmax = 360,
    cmin = 20,
    cmax = 50,
    lmin = 50,
    lmax = 100,
    aesthetics = c("fill", "color")
  )

cairo_pdf("OUT/tree.pdf",
  width = 16,
  height = 16
)
plot_3
dev.off()

# in case you want to continue with the filtered dataset, do this!

meCom =
  prevFilt


summarize_phyloseq(meCom)

```

### Relative Abundance

```{r}

# function for relative abundance

plotRelativeAbundance = function(level, thresh, var1, var2) {
  meCom %>%
    transform("compositional") %>%
    speedyseq::psmelt() %>%
    group_by(
      {{ var1 }},
      {{ var2 }},
      {{ level }}
    ) %>%
    summarise(Abundance = sum(Abundance)) %>%
    mutate(taxa = ifelse(
      Abundance >= quantile(Abundance, thresh),
      {{ level }},
      paste("below q", thresh, sep = "")
    )) %>%
    ggplot() +
    geom_bar(
      stat = "identity",
      position = "fill",
      aes(
        x = {{ var1 }},
        y = Abundance,
        fill = taxa
      )
    ) +
    scale_fill_iwanthue(
      hmin = 0,
      hmax = 360,
      cmin = 20,
      cmax = 50,
      lmin = 25,
      lmax = 100,
      random = T
    ) +
    facet_wrap(vars({{ var2 }}),
      scales = "free_x"
    ) +
    theme(
      legend.position = "right",
      axis.title.x = element_blank(),
      legend.key = element_rect(
        size = 0,
        linetype = 0
      ),
      axis.title.y = element_blank()
    )
}

# plot relative abundance per Phylum, grouping Phyla below the X quantile and displayed by var1 and var2 ()

plot_4 =
  plotRelativeAbundance(Family, 0.75, Bacteria, Description)

# save all in a pdf

cairo_pdf("OUT/relativeAbundance.pdf",
  width = 16,
  height = 16,
  onefile = T
)
plot_4
dev.off()

```


## DIVERSITY

### Richness and Alpha Diversity

Here we switch from plug-in to estimated measures of diversity (Shannon, Simpson)

```{r}
data("GlobalPatterns")

glo_patt = 
  subset_samples(GlobalPatterns, SampleType == "Freshwater" | SampleType == "Soil")

plot_1 = sample_richness(glo_patt) %>% 
  summary() %>% 
  as_tibble() %>% 
  left_join(meta(glo_patt) %>% 
              rownames_to_column("sample_names")) %>% 
  ggplot(
    aes(
      x = SampleType,
      y = estimate,
      lower = lower,
      upper = upper,
      fill = SampleType
    )
  ) +
  geom_boxplot(size = 0.5) +
  scale_fill_manual(values = palette_list$Description) +
  ylab("observed richness") +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_blank()
  ) +
  geom_signif(
    comparisons = list(c("Freshwater", "Soil")),
    map_signif_level = F,
    tip_length = 0,
    test = "wilcox.test",
    step_increase = 0.05
  ) 
plot_2 =
  breakaway(glo_patt) %>% 
  summary() %>% 
  as_tibble() %>% 
  left_join(meta(glo_patt) %>% 
              rownames_to_column("sample_names")) %>% 
  ggplot(
    aes(
      x = SampleType,
      y = estimate,
      fill = SampleType
    )
  ) +
  geom_boxplot(size = 0.5) +
  scale_fill_manual(values = palette_list$Description) +
  ylab("breakaway estimate") +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_blank()
  )  +
  geom_signif(
    comparisons = list(c("Freshwater", "Soil")),
    map_signif_level = F,
    tip_length = 0,
    test = "wilcox.test",
    step_increase = 0.05
  ) 

## estimates

res_breakaway =
  breakaway(meCom) %>% 
  summary() %>% 
  as_tibble() %>% 
  left_join(meta(meCom) %>% 
              rownames_to_column("sample_names"))

a =
  ggplot(
    data = res_breakaway,
    aes(
      x = Bacteria,
      y = estimate,
      lower = lower,
      upper = upper,
      fill = Bacteria
    )
  ) +
  geom_boxplot(size = 0.5) +
  scale_fill_manual(values = palette_list$Bacteria) +
  ylab("breakaway estimate") +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_blank()
  ) +
  geom_signif(
    comparisons = list(
      c("WT1", "LjR176"),
      c("LjR87", "LjR176"),
      c("WT1", "LjR87")
    ),
    map_signif_level = F,
    tip_length = 0,
    test = "wilcox.test",
    step_increase = 0.05
  ) 

# your design matrix

des_mat =
  model.matrix(~ Description + Bacteria, data = meta(meCom))

# calculate diversity estimates for your samples

set.seed(666)

res_divnet =
  meCom %>%
  divnet(
    ncores = parallel::detectCores(),
    tuning = "fast")

# group results

alpha_diversity_estimates =
  map(
    c("shannon", "simpson") %>%
      set_names(),
    ~ res_divnet %>%
      pluck(.x) %>%
      summary() %>%
      left_join(meta(meCom) %>%
        rownames_to_column("sample_names"),
      by = "sample_names"
      )
  )

# shannon estimates boxplot

b = 
  ggplot(
    data = alpha_diversity_estimates$shannon,
    aes(
      x = Bacteria,
      y = estimate,
      lower = lower,
      upper = upper,
      fill = Bacteria
    )
  ) +
  geom_boxplot(size = 0.5) +
  scale_fill_manual(values = palette_list$Bacteria) +
  ylab("shannon estimate") +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_blank()
  ) +
  geom_signif(
    comparisons = list(
      c("WT1", "LjR176"),
      c("LjR87", "LjR176"),
      c("WT1", "LjR87")
    ),
    map_signif_level = F,
    tip_length = 0,
    test = "wilcox.test",
    step_increase = 0.05
  )

# simpson estimates boxplot

c = 
  ggplot(
    data = alpha_diversity_estimates$simpson,
    aes(
      x = Bacteria,
      y = estimate,
      lower = lower,
      upper = upper,
      fill = Bacteria
    )
  ) +
  geom_boxplot(size = 0.5) +
  scale_fill_manual(values = palette_list$Bacteria) +
  ylab("simpson estimate") +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_blank()
  ) +
  geom_signif(
    comparisons = list(
      c("WT1", "LjR176"),
      c("LjR87", "LjR176"),
      c("WT1", "LjR87")
    ),
    map_signif_level = F,
    tip_length = 0,
    test = "wilcox.test",
    step_increase = 0.05
  ) 

# all together

plot_5 =
  a + b + c +
  plot_layout(guides = "collect")

cairo_pdf("OUT/alphaDiversity.pdf",
  width = 16,
  height = 16,
  onefile = T
)
plot_5
dev.off()

## hypothesis testing

res_simpson =
  testDiversity(res_divnet,
    h0 = "simpson"
  )

#  p-value for global test: 2.6132296326864e-10

res_shannon =
  testDiversity(res_divnet,
    h0 = "shannon"
  )

# p-value for global test: 0 

```

### Beta Diversity

```{r}

## testing differences in distance

beta_diversity_estimates =
  map(
    c("bray-curtis", "euclidean") %>%
      set_names(),
    ~ res_divnet %>%
      pluck(.x) %>%
      as.data.frame() %>%
      rownames_to_column("sample_names") %>% 
      left_join(meta(meCom) %>%
        rownames_to_column("sample_names"),
      by = "sample_names"
      )
  )

res_means =
  compare_means(
    value ~ Bacteria,
    data = 
      beta_diversity_estimates$`euclidean` %>%
      pivot_longer(2:(nsamples(meCom) + 1)),
    method = "wilcox.test",
    p.adjust.method = "BH"
  )

write_tsv(res_means, "OUT/aitchison_wilcox.tsv")

## PCA 

res_prcomp =
  prcomp(
    meCom %>% 
      abundances("clr") %>%
      t()
  )
biplot(res_prcomp)

plot_6 =
  res_prcomp$x %>%
  as.data.frame() %>%
  rownames_to_column("sample_names") %>%
  left_join(meta(meCom) %>%
    rownames_to_column("sample_names")) %>%
  ggplot(
    aes(
      x = PC1,
      y = PC2
    )
  ) +
  geom_point(
    aes(
      group = Bacteria,
      color = Bacteria,
      shape = Description
    ),
    size = 7
  ) +
  scale_shape_manual(values = c(1, 2)) +
  stat_ellipse(
    linetype = 0,
    geom = "polygon",
    type = "norm",
    alpha = 0.05,
    aes(fill = Bacteria)
  ) +
  scale_color_manual(
    values = palette_list$Bacteria,
    aesthetics = c("fill", "color")
  ) +
  labs(
    x = paste("Dim. 1,", round(res_prcomp$sdev[1]^2 / sum(res_prcomp$sdev^2), 2)),
    y = paste("Dim. 2,", round(res_prcomp$sdev[2]^2 / sum(res_prcomp$sdev^2), 2))
  ) +
  theme(legend.position = "right") +
  geom_point(
    data =
      res_prcomp$x %>%
        as.data.frame() %>%
        rownames_to_column("sample_names") %>%
        left_join(meta(meCom) %>%
          rownames_to_column("sample_names")) %>%
        group_by(Bacteria) %>%
        summarise(
          PC1 = mean(PC1),
          PC2 = mean(PC2)
        ),
    aes(
      x = PC1,
      y = PC2,
      color = Bacteria
    ),
    size = 5,
    shape = 19
  ) +
  geom_polygon(
    data =
      res_prcomp$x %>%
        as.data.frame() %>%
        rownames_to_column("sample_names") %>%
        left_join(meta(meCom) %>%
          rownames_to_column("sample_names")) %>%
        group_by(Bacteria) %>%
        summarise(
          PC1 = mean(PC1),
          PC2 = mean(PC2)
        ) %>%
        mutate(group = "temp"),
    aes(
      x = PC1,
      y = PC2,
      group = group
    ),
    alpha = 0,
    fill = "white",
    color = "#aaaaaa",
    size = 0.5,
    linetype = "dotted"
  )

cairo_pdf("OUT/PCA.pdf",
  width = 16,
  height = 16
)
plot_6
dev.off()

```

## DIFFERENTIAL ABUNDANCE

### ANCOM-BC

```{r}

# create unique identifier for Genus

new_tax_table = 
  meCom@tax_table@.Data %>% 
  as.data.frame() %>% 
  mutate(Genus = vctrs::vec_as_names(Genus, repair = "unique") %>% 
           sub("...", "_", fixed = T, .))

# workaround for something I dindn't have the patience for

meCom_id = 
  phyloseq(
    otu_table(meCom@otu_table, taxa_are_rows = T),
    sample_data(meCom@sam_data),
    tax_table(new_tax_table %>% 
                as.matrix))

## M26 vs AR10-3-9

res_ancombc =
  ancombc(meCom_id,
  formula = "Description",
  p_adj = "fdr",
  zero_cut = 1,
  conserve = T
  )

res_ancombc_df =
  res_ancombc$res$beta %>%
  as.data.frame() %>%
  rownames_to_column("ASV") %>%
  left_join(
    res_ancombc$res$q_val %>%
      as.data.frame() %>%
      rownames_to_column("ASV"),
    by = "ASV",
    suffix = c("_beta", "_q")
  ) %>% 
  left_join(
    res_ancombc$res$W %>% 
      as.data.frame() %>% 
      rownames_to_column("ASV")) %>% 
  dplyr::rename(DescriptionInoculated_W = DescriptionInoculated) %>% 
  mutate(DescriptionInoculated_DA = as.factor(case_when(
    DescriptionInoculated_q > 0.01 & DescriptionInoculated_q <= 0.05 ~ "*",
    DescriptionInoculated_q > 0.001 & DescriptionInoculated_q <= 0.01 ~ "**",
    DescriptionInoculated_q > 0.0001 & DescriptionInoculated_q <= 0.001 ~ "***",
    DescriptionInoculated_q <= 0.0001 ~ "****",
    TRUE ~ "ns"
  ))) %>% 
  pivot_longer(-ASV,
    names_to = c("var", ".value"),
    names_sep = "_") %>% 
  mutate(var = factor(var, levels = c("DescriptionInoculated"), 
                      labels = c("<- Control|Treatment ->")))

res_ancombc_df$DA = relevel(res_ancombc_df$DA, ref = "ns")

relAbu_ancombc = 
  res_ancombc_df %>%
  filter(q <= 0.05) %>% 
  left_join(tax_table(meCom_id) %>% 
              as.data.frame() %>% 
              rownames_to_column("ASV"))

plot_7 =
  ggplot(
    res_ancombc_df
  ) +
  geom_point(
    aes(
      x = beta,
      y = W,
      color = DA,
      alpha = DA
    ),
    size = 2,
    shape = 19
  ) +
  scale_color_manual(values = michelangelo$grad3) +
  facet_wrap(~var) +
  theme(legend.position = "right") +
  scale_alpha_discrete(range = c(0.2, 1))

plot_8 = 
  relAbu_ancombc  %>% 
  ggplot() +
  geom_bar(aes(y = fct_reorder(Genus, W),
               x = W, 
               fill = Family),
           stat = "identity",
           position = "dodge") +
  scale_fill_iwanthue(hmin = 0,
                      hmax = 360,
                      cmin = 0,
                      cmax = 30,
                      lmin = 50,
                      lmax = 75) +
  theme(axis.text.y = element_text(size = 12),
        axis.title.y = element_blank(),
        legend.position = "right",
        strip.text = element_blank()) +
  guides(fill = guide_legend(ncol = 1)) +
  facet_wrap(~var) +
  geom_vline(xintercept = 0,
             size = 0.5,
             color = "#696969")

cairo_pdf("OUT/ANCOMBC.pdf",
  width = 16,
  height = 16,
  onefile = T
)
plot_7 / plot_8
dev.off()

```

## NetComi

```{r}

p_load(NetCoMi)

res_netCon =
  netConstruct(
    data = meCom,
    measure = "spring",
    cores = 6,
    measurePar = list(nlambda = 10,
                      rep.num = 10),
    dissFunc = "signed",
    verbose = 3,
    seed = 666
  )

res_netAna = 
  netAnalyze(res_netCon, 
             centrLCC = TRUE,
             clustMethod = "cluster_fast_greedy",
             hubPar = "eigenvector",
             weightDeg = FALSE,
             normDeg = FALSE)

summary(res_netAna)

colPal = iwanthue(n = 12, 0, 360, 25, 50, 50, 75)

plot(res_netAna,
     layout = "circle",
     labels = F,
     nodeSize = "eigenvector",
     nodeColor = "cluster",
     colorVec = colPal,
     nodeTransp = 1,
     borderCol = "#696969",
     borderWidth = 0.5,
     hubBorderWidth = 1.5,
     hubBorderCol = "#111111",
     edgeWidth = 0.5,
     posCol = c("powderblue"),
     negCol = c("indianred"),
     title1 = "Network on OTU level with SPRING associations",
     showTitle = T)
```