articletwo.Rmd

---
title: "Hypertension and microbiome"
author: "Joonatan Palmu"
date: "`r format(Sys.time(), '%d.%m.%Y')`"
output: 
  html_document:
    number_sections: true
---
	
```{r setup, echo=FALSE}
knitr::opts_chunk$set(include = TRUE, echo = TRUE, message = FALSE, results='asis',
                      cache=FALSE, warning=FALSE)
knitr::opts_chunk$set(cache.path = 'cache/', output.dir="cache/",
                      file.path = 'cache/', fig.path = 'cache/')

options(knitr.kable.NA = "")

dir.create("cache/", showWarnings = FALSE)
dir.create("rds/", showWarnings = FALSE)
dir.create("session/", showWarnings = FALSE)

now <- format(Sys.time(), '%Y%m%d-%H%M%S')
```

```{r Command line arguments, echo = FALSE, eval = FALSE}
file <- paste0("session/", sort(list.files("session"), decreasing = TRUE)[1])
message("Loading variables from file ", file)
load(file)
```

# Importing libraries:

<details><summary>Libraries</summary>

```{r libraries, cache = FALSE}
library(dplyr)
library(tibble)
library(phyloseq)
library(nortest)
library(microbiome)
library(knitr)
library(tidyr)
library(vegan)
library(reshape)
library(parallel)
library(officer)
library(flextable)
library(xtable)
library(rvg)
library(tableone)
library(scales)
library(ggplot2)
library(egg)
library(gridExtra)
library(png)
library(ggpubr)
library(broom)
library(ggfortify)
library(RColorBrewer)
library(gvlma)
library(purrr)
library(gtable)
library(car)
library(M3C)
library(emmeans)
library(DESeq2)
library(ggeffects)
library(pander)
library(hexbin)
library(readr)
library(KEGGREST)
```

</details>

# Session info

<details><summary>Session</summary>

```{r Session info}
pander(sessionInfo(), compact = TRUE)
```

</details>

# Import helper functions

<details><summary>Functions</summary>

```{r Functions}
sourcefiles <- c("articletwo-officer.R",
                 "articletwo-import.R",
                 "articletwo-rrbiome.R",
                 "articletwo-ggplot.R")
```

```{r import files, echo = FALSE}
for (f in sourcefiles) {
    source(f)
}
```

```{r embbed files, echo = FALSE}
xfun::embed_files(c("rrnmr.Rmd", sourcefiles))
```

</details>


# Loading data
	
Loading descriptions for clinical data

```{r variables, warning = FALSE}
names.dset <- getdescriptions()
```

Loading phyloseq object

```{r data}
pseq.species <- import_filter_data("data/phfinrisk_species_all_drop50k_2018-12-21.RDs")
pseq.genus <- import_filter_data("data/phfinrisk_genus_all_drop50k_2018-11-16.RDs")
pseq.genus.coretaxa <- coretaxa(pseq.genus, detection = 0.1/100, prevalence = 1/100)
pseq.genus.core.deseq <-  prune_taxa(pseq.genus.coretaxa, pseq.genus)
```

At species level meta has dimensions (`r dim(meta(pseq.species))`) and
there ntaxa is `r ntaxa(pseq.species)`. At genus level meta has
dimensions (`r dim(meta(pseq.genus))`) and there ntaxa is 
`r ntaxa(pseq.genus)`. Core size is `r length(pseq.genus.coretaxa)`.

# Variables

```{r my variables}
var.BP <- c("MAP", "SYSTM", "DIASM", "PULSEPRESSURE", "HYPERTENSION")
var.CL.min <- c("BL_AGE", "SEX")
var.CL <- c("BL_AGE", "SEX", "BMI", "CURR_SMOKE", "Q57X", "PREVAL_DIAB",
            "BL_USE_RX_C03","BL_USE_RX_C07", "BL_USE_RX_C08", "BL_USE_RX_C09")
```

# Heavy models

## Bray curtis distance matrix

```{r matrix calculation}
if (!file.exists("rds/bray.dist.m.species.rds")) {
    bray.dist.m.species <- calculate.beta.matrix(pseq.species)
    saveRDS(bray.dist.m.species, file = "rds/bray.dist.m.species.rds")
} else {
    bray.dist.m.species  <- readRDS("rds/bray.dist.m.species.rds")
}
```

## Beta diversity

 ```{r adonis calculation}
 if (!file.exists("rds/adonis.species.rds")) {
     adonis.species <- calculate.betadiversity(pseq = pseq.species,
                                               matrix = bray.dist.m.species,
                                               npermutations = 999,
                                               vars = list("max" = var.CL,
                                                            "min" = var.CL.min))
     saveRDS(adonis.species, file = "rds/adonis.species.rds")
 } else {
     adonis.species  <- readRDS("rds/adonis.species.rds")
 }
```

## PCoA

```{r pcoa calculate}
if (!file.exists("rds/pcoa.ordinate.rds")) {
    pcoa.abundances <- microbiome::transform(pseq.species, 'compositional')
    pcoa.ordinate <- ordinate(pcoa.abundances, method="PCoA", distance="bray")
    saveRDS(pcoa.ordinate, file = "rds/pcoa.ordinate.rds")
} else {
    pcoa.ordinate <- readRDS("rds/pcoa.ordinate.rds")
}
```

## DeSeq2

Associations for common genera

```{r deseq2}
if (!file.exists("rds/dds.rds")) {
    dds <- lapply(c2l(var.BP), function(x, pseq) {
        dds.data <- phyloseq_to_deseq2(pseq, deseq.formula(x, var.CL))
        DESeq(dds.data,
              test="Wald",
              fitType="parametric",
              parallel = TRUE,
              BPPARAM=MulticoreParam(16))
    }, pseq = pseq.genus.core.deseq)
    saveRDS(dds, file = "rds/dds.rds")
} else {
    dds <- readRDS(file = "rds/dds.rds")
}
```

The affect of third covariate for hypertension

```{r deseq2 third covariate}    
if (!file.exists("rds/dds3.rds")) {
    dds3 <- lapply(deseq.list(var.CL, var.CL.min), function(x, pseq) {
        dds.data <- phyloseq_to_deseq2(pseq, deseq.formula(x, "HYPERTENSION"))
        DESeq(dds.data,
              test="Wald",
              fitType="parametric",
              parallel = TRUE,
              BPPARAM=MulticoreParam(16))
    }, pseq = pseq.genus.core.deseq)
    saveRDS(dds3, file = "rds/dds3.rds")
}  else {
    dds3 <- readRDS(file = "rds/dds3.rds")
}
```

Associations for Lactobacili

```{r deseq2 lactobacilli}
if (!file.exists("rds/ddsLactobacilli.rds")) {
    ddsLactobacilli <- mclapply(c2l(var.BP), function(x) {
        myDESeq(pseq = pseq.species,
                vars = c(x, var.CL),
                coreterm = "^Lactobacillus.*Bacteria")
    }, mc.cores = length(var.BP))
    saveRDS(ddsLactobacilli, file = "rds/ddsLactobacilli.rds")
} else {
    ddsLactobacilli <- readRDS(file = "rds/ddsLactobacilli.rds")
}
```

Subgroup analyses for Lactobacilli

```{r heavy subgroup analyses for lactobasilli-BP}
if (!file.exists("rds/ddsSubLactobacilli.rds")) {
    ddsSubLactobacilli <- lapply(myDESeq.subsets(), function(subset) {
        mclapply(c2l(var.BP), function(bp) {
            myDESeq(pseq = pseq.species,
                    vars = bp %union% var.CL %difference% subset$exclude,
                    coreterm = "^Lactobacillus.*Bacteria",
                    FUN = subset$fun)
        }, mc.cores = length(var.BP))
    })
    saveRDS(ddsSubLactobacilli, file = "rds/ddsSubLactobacilli.rds")
} else {
    ddsSubLactobacilli <- readRDS(file = "rds/ddsSubLactobacilli.rds")
}

#ddsSubLactobacilli <-  readRDS(file = "rds/ddsSub-temp.rds") %>%
#    purrr::list_modify(drug = readRDS(file = "rds/ddsSubLactobacilli.rds")$drug)
```


# Characteristics

```{r Characteristics}
characteristicsTable(meta(pseq.species)) %>% kable
```

<details><summary>Salt subset</summary>

```{r Characteristics subset, echo = FALSE}
pseqsubset(pseq.genus, coretaxa = pseq.genus.coretaxa, saltsubset = TRUE) %>%
    meta %>%
    characteristicsTable %>%
    kable
```

</details>

<details><summary>Chi-square tests and t-tests</summary>

```{r chi squared}
lapply(c2l(var.CL, var.BP), function(x) {
    meta(pseq.species) %>%
        mutate(nagroup = ifelse(is.na(dUNA), 1, 0)) %$%
        { if (x %in% c("BL_AGE", "BMI", "MAP", "SYSTM", "DIASM", "PULSEPRESSURE"))
              a <<- t.test(as.formula(sprintf("%s ~ nagroup", x))) 
          else
              table(nagroup, .[[x]]) %>%
                  chisq.test() } %>%
        tidy}) %>%
    map_df(~.x, .id = "term") %>%
    mutate (p.value = pub.p(p.value)) %>%
    select(term, p.value)
```

```{r chi squared multilevel}
lapply(c2l(1, 2, 3), function(x) 
    meta(pseq.species) %>%
    mutate(nagroup = ifelse(is.na(dUNA), 1, 0),
           Q57X = ifelse(Q57X == x, 1, 0)) %$%
    table(nagroup, Q57X) %>%
    chisq.test() %>%
    tidy) %>%
    map_df(~.x, .id = "term") %>%
    mutate (p.value = pub.p(p.value)) %>%
    select(term, p.value)
```

</details>

# Alpha and beta diversity plots

```{r alphabeta definitions}
diversity <- diversities(pseq = pseq.species,
                           vars = list("max" = var.CL, "min" = var.CL.min),
                           betadiversity = adonis.species,
                           names.dset = names.dset)
```

```{r talbe alpha div min, echo = FALSE}
diversities.tidy(diversity) %>%
    kable
```

Diversity plot

```{r save grob, echo = FALSE}
g.diversity <- plot.diversities(diversity)
ggsave(file = "cache/alpha-beta.png", plot = g.diversity, height = 6, width = 9)
```

<img src="cache/alpha-beta.png" />

# Principal coordinate analysis

```{r pcoa plot, echo = FALSE}
pcoa.plot <- pcoaplot(pcoa.ordinate)
ggsave(file = "cache/pcoa-species.png", plot=pcoa.plot, height=3.0, width=6, units = "in", dpi = 300)
```

<img src="cache/pcoa-species.png" />

# Associations between genera and blood pressure indices
	
```{r deseq continuous}
dset.deseq.signf <- deseqresults(dds, names.dset)
g.deseq <- deseqheatmap(dset.deseq.signf)
```

Significant associations found between blood pressure indices and N:o `r nrow(dset.deseq.signf)`
genera and N:o `r dset.deseq.signf %>% pull(Feature) %>% unique %>% length` features. N:o 
`r dset.deseq.signf %>% filter(log2FoldChange > 0) %>% nrow` positive associations found.

## N:o associations for hypertension per included 3rd covariate

Base model had 
`r sum(results(dds3$BL_AGE, name = "HYPERTENSION_1_vs_0")$padj < 0.05, na.rm = TRUE)`
significant associations and inclusion of BMI dropped it to
`r sum(results(dds3$BMI, name = "HYPERTENSION_1_vs_0")$padj < 0.05, na.rm = TRUE)`.

```{r results compare deseq2}
deseqresults(dds3, names.dset, hypertensionvars = names(dds3)) %>%
    select(Name, Feature, log2FoldChange) %>%
    spread(Name, log2FoldChange) %>%
    kable
```

# Associations between Lactobacilli and blood pressure indices

```{r lactobacillus species results}
dset.deseq.signf.na <- deseqresults(ddsLactobacilli, names.dset)
g.deseq.na <- deseqheatmap(dset.deseq.signf.na)
```

Significant associations found between blood pressure indices and N:o `r nrow(dset.deseq.signf.na)`
genera and N:o `r dset.deseq.signf.na %>% pull(Feature) %>% unique %>% length` features. N:o 
`r dset.deseq.signf.na %>% filter(log2FoldChange > 0) %>% nrow` positive associations found.

```{r deseq dual plot, echo = FALSE}
g.deseqgenus <- gtable_frame(ggplotGrob(g.deseq + theme(legend.position = "none")),
                  height = unit(length(unique((dset.deseq.signf$Feature))), 'mm'),
                  width = unit(25, 'mm'))

g.deseqspecies <- gtable_frame(ggplotGrob(g.deseq.na + theme(legend.position = c(0.7, 1.45))),
                  height = unit(length(unique((dds.results$Feature))), 'mm'),
                                width = unit(25, 'mm'))

g.deseq <- arrangeGrob(g.deseqgenus, g.deseqspecies, ncol = 2)

ggsave(file = "cache/deseq-full.png", plot = g.deseq, width = 14, height = 25, dpi = 300, unit = "cm")
```

<img src="cache/deseq-full.png" />


# Association between Lactobacillus genus and urinary sodium

```{r significant results for salt subset}
(dds.salt <- myDESeq(pseq = pseq.genus,
                    coretaxa = pseq.genus.coretaxa,
                    vars = c("dUNA", var.CL),
                    saltsubset = TRUE)) %>%
    results(., name = "dUNA", tidy = TRUE) %>%
    dplyr::mutate(qval = p.adjust(pvalue, method="BH"),
                  lfc_se = sprintf("%.4f±%.4f", log2FoldChange, lfcSE)) %>%
    select(row, lfc_se, qval) %>%
    filter(row == "Lactobacillus (Bacteria)") %>%
    kable
```

Graphical representation for findings

```{r salt subset char}
pseq.genus.salt <- pseqsubset(pseq.genus, coretaxa = pseq.genus.coretaxa, saltsubset = TRUE)
df.genus.salt <- saltdset(pseq.genus.salt, dds.salt)
```

```{r genus salt levels}
df.genus.salt %>%
    group_by(group) %>%
    summarize(mean = round(mean(dUNA), 1),
              sd = round(sd(dUNA), 1)) %>%
    kable
```

Mean urinary salt `r df.genus.salt$dUNA %>% mean` and sd `r df.genus.salt$dUNA %>% sd`.

# Association between Lactobacilli and urinary sodium

Analysis between *Lactobacilli* and urinary salt. 

```{r species salt}
( dds.salt.species <- myDESeq(pseq = pseq.species,
                        vars = c("dUNA", var.CL),
                        coreterm = "^Lactobacillus.*Bacteria",
                        saltsubset = TRUE) ) %>%
    results(., name = "dUNA", tidy = TRUE) %>%
    dplyr::mutate(qval = p.adjust(pvalue, method="BH"),
                  row = renametaxa(row)) %>%
    filter(qval < 0.05) %>%
    mutate(qval_exp = sprintf("%.4e", qval)) %>%
    kable
```

```{r test}
dds.salt.species %>%
      results(., name = "dUNA", tidy = TRUE) %>%
    dplyr::mutate(qval = p.adjust(pvalue, method="BH"),
                  row = renametaxa(row)) %>%
                                        #    filter(qval < 0.05) %>%
    arrange(qval) %>%
    mutate(qval_exp = sprintf("%.4e", qval)) %>%
    kable
```

```{r species salt}
dds.salt.species %>%
    results(., name = "SEX_1_vs_0", tidy = TRUE) %>%
    dplyr::mutate(qval = p.adjust(pvalue, method="BH"),
                  row = renametaxa(row)) %>%
    filter(qval < 0.05) %>%
    mutate(qval_exp = sprintf("%.4e", qval)) %>%
    kable
```

The box plot for *Lactobacilli* groups and urinary salt to demonstrate the association.

```{r species plot}
( df.species.salt <- deseq.abundances(dds.salt.species) %>%
    full_join(deseq.abundances(dds.salt), by = "sampleid") %>%
    left_join(meta(pseq.species) %>% tibble::rownames_to_column("sampleid"), by = "sampleid") %>%
    mutate_at(vars(contains('Lactobacillus')), .funs = list(group = ~factor(dplyr::ntile(., 4)))) %>%
    select(dUNA,
           Lactobacillus_paracasei.Bacteria_group,
           Lactobacillus_salivarius.Bacteria_group,
           Lactobacillus.Bacteria_group) ) %>%
    gather(taxa, abundance, -dUNA) %>%
    saltboxplot %>%
    ggsave(file = "cache/boxplot-lacto-salt.png", plot = ., height = 3, width = 3.5, dpi = 300, unit = "in")
```

<img src = "cache/boxplot-lacto-salt.png" />

Mean and SD for urinary sodium quartiles

```{r species salt levels paracasei, echo = FALSE}
lapply(c2l(df.species.salt %>% select(contains("Lactobacillus")) %>% colnames),
       function(x) df.species.salt %>%
                   dplyr::rename(group = !!x) %>%
                   group_by(group) %>%
                   summarize(mean = round(mean(dUNA), 2),
                             sd = round(sd(dUNA), 2)) ) %>%
    map_df(~as.data.frame(.x), .id = "taxa")  %>%
    mutate(taxa = renametaxa(taxa)) %>%
    kable
```

# Subgroup analyses

## Sex

```{r deseq dual plot, echo = FALSE}
g.ddssub.male <- ddsSubLactobacilli$male %>%
    deseqresults(., names.dset) %>%
    deseqheatmap(breaks = seq(-2, 2)) %>%
    { . + theme(legend.position = "none",
                plot.title = element_text(size = 12, face = "plain")) } %>%
    myggtitle(., title= "A. Men")

g.ddssub.female <- ddsSubLactobacilli$female %>%
    deseqresults(., names.dset) %>%
    deseqheatmap(breaks = seq(-2, 2)) %>%
    { g.ddssub.female.raw <<- . } %>%
    { . + theme(legend.position = "none",
                 plot.title = element_text(size = 12, face = "plain")) } %>%
    myggtitle(., title= "B. Women")

g.ddssub.legend <- get_legend(g.ddssub.female.raw) %>% as_ggplot

g.deseq.subsample.sex <- arrangeGrob(g.ddssub.male,
                                     g.ddssub.female,
                                     g.ddssub.legend ,
                                     ncol = 3,
                                     widths = c(2, 2, 1))
ggsave(file = "cache/deseq-subsample-sex.png",
       plot = g.deseq.subsample.sex,
       width = 17,
       height = 11,
       dpi = 300,
       unit = "cm")
```

## Antihypertensive

```{r deseq dual plot, echo = FALSE}
g.ddssub.nodrug <- ddsSubLactobacilli$nodrug %>%
    purrr::list_modify(HYPERTENSION = NULL) %>%
    deseqresults(., names.dset) %>%
    deseqheatmap(breaks = seq(-2, 2)) %>%
    { . + theme(legend.position = "none",
                 plot.title = element_text(size = 12, face = "plain")) } %>%
    myggtitle(., title= "A. No antihypertensive medication use")

g.ddssub.drug <- ddsSubLactobacilli$drug %>%
    purrr::list_modify(HYPERTENSION = NULL) %>%
    deseqresults(., names.dset) %>%
    deseqheatmap(breaks = seq(-2, 2)) %>%
    { . + theme(legend.position = "none",
                plot.title = element_text(size = 12, face = "plain")) } %>%
    myggtitle(., title= "B. Antihypertensive medication use")

g.deseq.subsample.drug <- arrangeGrob(g.ddssub.nodrug,
                                      g.ddssub.drug,
                                      g.ddssub.legend,
                                     ncol = 3,
                                     widths = c(2, 2, 1))
ggsave(file = "cache/deseq-subsample-drug.png",
       plot = g.deseq.subsample.drug,
       width = 19,
       height = 14.2,
       dpi = 300,
       unit = "cm")
```


# Sodium subsample

```{r subset salt sex}
( ddsSubLactobacilliSodium <- mclapply(myDESeq.subsets(), mc.cores = 4, function(subset)
    myDESeq(pseq = pseq.species,
            vars = "dUNA" %union% var.CL %difference% subset$exclude,
            coreterm = "^Lactobacillus.*Bacteria",
            saltsubset = TRUE,
            FUN = subset$fun)) ) %>%
    myDESeq.tidy() %>%
    map_df(~.x, .id = "subgroup") %>%
    kable
```

```{r subset salt sex}
ddsSubLactobacilliSodium %>%
    myDESeq.tidy()) %>%
    map_df(~.x, .id = "subgroup") %>%
    kable
```

## Lactobacilli prevalences

Prevalence

```{r prevalence 0.01%}
lapply(pseq.species %>% taxa %>% mygrep(word = "^Lactobacillus_", .) %>% c2l, function(x) 
    pseq_prevalence(pseq.species, taxa = x, limit = 0.01/100)) %>%
    tibble::enframe() %>% 
    unnest(cols = value) %>%
    arrange(desc(value)) %>%
    kable


females <- pseq.species %>%
    meta %>%
    tibble::rownames_to_column("sampleid") %>%
    filter(SEX == 1) %>%
    pull(sampleid)

males <- pseq.species %>%
    meta %>%
    tibble::rownames_to_column("sampleid") %>%
    filter(SEX == 0) %>%
    pull(sampleid)

taxa(pseq.species) %>% mygrep(word = "Salivarius")

pseq_prevalence(prune_samples(males, pseq.species), limit = 0.01/100, taxa = "Lactobacillus_paracasei (Bacteria)")
pseq_prevalence(prune_samples(females, pseq.species), limit = 0.01/100, taxa = "Lactobacillus_paracasei (Bacteria)")


```

```{r prevalence 0.1%}
lapply(c("Lactobacillus_paracasei (Bacteria)", "Lactobacillus_salivarius (Bacteria)") %>% c2l, function(x) 
    pseq_prevalence(pseq.species, taxa = x, limit = 0.1/100)) %>%
    tibble::enframe() %>% 
    unnest(cols = value) %>%
    arrange(desc(value)) %>%
    kable
```


# Functional analysis

```{r import KEGG data}
dset.ko  <- readRDS("data/function_activity_copied_20200213.rds")$ko %>%
    as.data.frame %>%
    rownames_to_column(var = "term") %>%
    gather(sampleid, value, -term) %>% 
    spread(term, value) %>%
    mutate_at(vars(starts_with("K")), ~scale(log10(. + 1))) %>%
    left_join(meta(pseq.species) %>% rownames_to_column(var = "sampleid"), by = "sampleid")
```

```{r model KEGG}
ret.ko <- loop.lm(dset = dset.ko,
               response = "SYSTM",
               loops = dset.ko %>% colnames %>% mygrep(word = "^K[0-9]+"),
               covariates = var.CL) %>%
    loop.results()
```

```{r process KEGG results}
ret.ko %>%
    filter(qval < 0.05) %>%
    arrange(qval) %>%
    mutate(mean_ci = sprintf("%.2f (%.2f to %.2f)", estimate, conf.low, conf.high),
           qval = pub.p(qval)) %>%
    mutate(desc = mykeggget(term)) %>%
    select(term, desc, mean_ci, qval) %>%
    kable
```

<details><summary>Diagnostic plots</summary>

```{r KEGG qq plots}
diagnosticqqplot(dset.ko, ret.ko %>% filter(qval < 0.05) %>% pull(term)) %>%
    ggsave(file = "cache/kegg.png", plot = ., height = 40, width = 10, dpi = 300, unit = "in")
```

<img src = "cache/kegg.png" />

```{r KEGG scatter plots}
diagnosticscatterplot(dset.ko,
                      ret.ko %>% filter(qval < 0.05) %>% pull(term),
                      SYSTM) %>%
    ggsave(file = "cache/kegg-scatter.png", plot = ., height = 40, width = 10, dpi = 300, unit = "in")
```

<img src = "cache/kegg-scatter.png" />

</details>

```{r process KEGG results}
ret.ko %>%
    filter(qval < 0.05, estimate < 0) %>%
    mutate(term = paste0("n-", term),
           size = paste0("v-", 1/qval)) %>%
    select(term, size) %>%
    write_tsv(path = "cache/functree_negative.tsv", col_names = FALSE)
```

# Supplements

Listing our "core" genera

```{r write supplement core names}
core.names.incols <- pseq.genus.coretaxa %>%
    renametaxa %>%
    sort %>%
    list.partition %>%
    mycbind %>%
    as.data.frame

core.names.flextable <- flextable(core.names.incols) %>%
    flextable::fontsize(size = 10, part = "header") %>%
    flextable::fontsize(size = 10, part = "body") %>%
    flextable::width(j=1:3, width = 1.8) %>%
    flextable::align(align = "left", part = "header") %>%
    flextable::align(align = "left")

core.names.flextable
```

Listing Lactobacillus spp

```{r write supplement lactobacilli}
lactobacilli.names.incols <- pseq.species %>%
    taxa %>%
    mygrep(word = "^Lactobacillus.*Bacteria") %>%
    renametaxa %>%
    sort %>%
    list.partition(parts = 4) %>%
    mycbind %>%
    as.data.frame

lactobacilli.names.flextable <- flextable(lactobacilli.names.incols) %>%
    flextable::fontsize(size = 10, part = "header") %>%
    flextable::fontsize(size = 10, part = "body") %>%
    flextable::width(j=1:3, width = 1.8) %>%
    flextable::align(align = "left", part = "header") %>%
    flextable::align(align = "left")

lactobacilli.names.flextable 
```


Listing results for alpha diversity

```{r supplement table for alpha diversity}
alphadiversity.table <- diversity %>% map_df(~as.data.frame(.x), .id = "model") %>%
    mutate(mean_ci = sprintf("%.2f (%.2f - %.2f)", alpha.effect, alpha.low, alpha.high),
           pubp = pub.p(alpha.p)) %>%
    myspread(list = c2l("mean_ci", "pubp"), term = "Name", key = "model") %>%
    arrange(match(Name, c("Systolic BP", "Diastolic BP", "Mean arterial pressure", "Pulse pressure", "Hypertension"))) %>%
        select(Name, min_mean_ci, min_pubp, max_mean_ci, max_pubp)

typology.tbls2 <- data.frame(
    col_keys = colnames(alphadiversity.table),
    what = c("", rep("Age- and sex adjusted model", 2), rep("Multivaraible adjusted model", 2)),
    measure = c("", rep(c("beta (95%-CI)", "p"), 2)),
    stringsAsFactors = FALSE)

alphadiversity.flextable <-
    typologyformatter(data = alphadiversity.table, font = 12, typology = typology.tbls2) %>%
    flextable::width(j=1, width = 1.9) %>%
    flextable::width(j=c(2:2,4:4), width = 1.8) %>%
    flextable::width(j=c(3:3,5:5), width = 0.7)
alphadiversity.flextable 
```


```{r supplement table for beta diversity}
betadiversity.table <- diversity %>% map_df(~as.data.frame(.x), .id = "model") %>%
    mutate(rsquared = sprintf("%.3f%%", beta.R2*100),
           pubp = pub.p(beta.p)) %>%
    myspread(list = c2l("rsquared", "pubp"), term = "Name", key = "model") %>%
    arrange(match(Name, c("Systolic BP", "Diastolic BP", "Mean arterial pressure", "Pulse pressure", "Hypertension"))) %>%
    select(Name, min_rsquared, min_pubp, max_rsquared, max_pubp)

typology.tbls3 <- data.frame(
    col_keys = colnames(betadiversity.table),
    what = c("", rep("Age- and sex adjusted model", 2), rep("Multivaraible adjusted model", 2)),
    measure = c("", rep(c("R2", "p"), 2)),
    stringsAsFactors = FALSE)

betadiversity.flextable <-
    typologyformatter(data = betadiversity.table, font = 12, typology = typology.tbls3) %>%
    flextable::width(j=1, width = 1.9) %>%
    flextable::width(j=c(2:5), width = 0.8)

betadiversity.flextable
```

```{r supplement table deseq}
deseq.table <- dset.deseq.signf %>%
    mutate(lfc_se = sprintf("%.2f±%.2f", log2FoldChange, lfcSE),
           p.value = pub.p(qval)) %>%
    myspread %>%
    select(Feature, starts_with("Systolic"), starts_with("Diastolic"), starts_with("Pulse pressure"),
           starts_with("MAP"), starts_with("Hypertension"))

typology.deseq <- data.frame(
    col_keys = colnames(deseq.table),
    what = c("", rep("Systolic BP", 2), rep("Diastolic BP", 2),
             rep("Pulse pressure", 2), rep("Mean arterial pressure", 2), rep("Hypertension", 2)),
    measure = c("", rep(c("Log2FC±SE", "p"), 5)),
    stringsAsFactors = FALSE)


typologyformatter(data = deseq.table, font = 9, typology = typology.deseq) %>%
    flextable::width(j = 1, width = 1.2) %>%
    flextable::width(j = seq(2, 11, 2), width = 0.75) %>%
    flextable::width(j = seq(3, 11, 2), width = 0.55)
```

```{r supplement table deseq na}
deseq.table.na <- dset.deseq.signf.na %>%
    mutate(lfc_se = sprintf("%.2f±%.2f", log2FoldChange, lfcSE),
           p.value = pub.p(qval)) %>%
    myspread %>%
    select(Feature, starts_with("Systolic"), starts_with("Diastolic"), starts_with("Pulse pressure"),
           starts_with("MAP"), starts_with("Hypertension"))

typologyformatter(data = deseq.table.na, font = 9, typology = typology.deseq) %>%
    flextable::width(j = 1, width = 1.2) %>%
    flextable::width(j = seq(2, 11, 2), width = 0.75) %>%
    flextable::width(j = seq(3, 11, 2), width = 0.55) 
```

```{r save session}
save.image(file = paste0("session/session-", now, ".Rdata"))
```