articleone.Rmd

---
title: "Eicosanoids mediators of hypertension"
author: "Joonatan Palmu"
date: "`r format(Sys.time(), '%d.%m.%Y')`"
---

```{r options, include=FALSE}
knitr::opts_chunk$set(include = TRUE, echo = TRUE, message = FALSE,
                      results='asis', cache=FALSE, warning=FALSE)
options(knitr.kable.NA = "")
now <- format(Sys.time(), '%Y%m%d-%H%M%S')
```

```{r read last session, include = FALSE, eval = FALSE}
message("Reading file ", paste0("rds/", sort(list.files("rds"), decreasing = TRUE)[1]))
load(paste0("rds/", sort(list.files("rds"), decreasing = TRUE)[1]))
```

```{r Remove temporary data, include = FALSE}
unlink("cache", recursive=TRUE)
dir.create("cache", showWarnings = FALSE)
dir.create("report", showWarnings = FALSE)
```

# Libraries

<details>
  <summary>Open/Close</summary>

```{r libraries, results = 'hide'}
library(knitr)
library(tidyr)
library(broom)
library(dplyr)
library(parallel)
library(xtable)
library(glmnet)
library(ggplot2)
library(scales)
library(magrittr)
library(forestplot)
library(survival)
library(png)
library(data.table)
library(tableone)
library(flextable)
library(ggrepel)
library(corrplot)
library(purrr)
library(gplots)
library(ggcorrplot)
library(Rtsne)
library(RColorBrewer)
library(gridExtra)
library(dendextend)
library(ggpubr)
library(readr)
library(heatmaply)
library(GeneNet)
library(Rgraphviz)
```


</details>

<details><summary>Session</summary>

```{r Session info}
sessionInfo()
```

</details>

# Functions

<details><summary>Functions</summary>

```{r Functions}
sourcefiles <- c("articleone-importer.R",
                 "articleone-functions.R",
                 "articleone-replication.R",
                 "articleone-riskscore.R",
                 "articleone-plots.R",
                 "articleone-tables.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>

# Data import

```{r import data}
dset <- importdata(replacenawithmin = TRUE, normalize = TRUE)
```

# Characteristics

```{r Characteristics}
characteristics(dset$data.full) %>% kable
```

# Linear models

## FR02

```{r Linear models}
linear.lmrank <- calculateglm(dset = dset$data.full, 
                                           loops = dset$metabolite.mzids,
                                           responses = c2l("MAP", "SBP", "DBP", "PP", "HT"),
                                           covariates = c("AGE", "sex", "BMI", "curr_smk", "curr_diab",
                                                          "HRX", "plate")) %>%
    mutate(mzrt = pub.mzrt.naive(term),
           qval = bonf.adjust(p.value, n=5))
```

Total number of `r dset$metabolite.mzids %>% length` eicosanoids
available and number of eicosanoids with significant associations is
`r pull.signf(linear.lmrank) %>% length`.

Number of associations per BP indice

```{r print number of associations per BP indice}
lapply(c2l("MAP", "SBP", "DBP", "PP", "HT"), function(x)
    linear.lmrank %>% filter(response == x)) %>%
    purrr::map_df(., ~as.data.frame(.x)) %>%
    group_by(response) %>%
    summarize(n = n(),
              significant = sum(qval < 0.05),
              proportion = significant/n,
              positive = sum(qval < 0.05 & estimate > 0),
              negative = sum(qval < 0.05 & estimate < 0)) %>%
    kable
```


Effect size of significant associations in FR02

```{r print linear results}
linear.lmrank %>%
    filter(qval < 0.05) %>%
    select(term, response, estimate) %>%
    mutate(estimate = round(estimate, 3)) %>%
    spread(response, estimate) %>% 
    kable
```

Number of significant associations to map

```{r number of associations for map}
filter.model(linear.lmrank, plimit = 0.05, filterterm = "SBP") %>%
    dplyr::mutate(sign = ifelse(estimate >0, 1, -1)) %>%
    group_by(sign) %>%
    summarize(n = n()) %>%
    mutate(prop = round(n/sum(n)*100, 1)) %>%
    kable
```

## Regularisation

```{r Stepwise calculation, results = 'hide'}
regularisation <- regularisationmodel(dset$data.full,
                                      pull.signf(linear.lmrank),
                                      bonfp = 0.05/length(dset$metabolite.mzids))
```


## Heatmap of metabolites

Correlation matrix

```{r heatmap for supplement}
my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 300)

eicosanoidcorrelation <- spearmancorrelation(dset, pull.signf(linear.lmrank, filterterm = "SBP"))

mark.position <- mygrep(regularisation$fwdbonf$term, word = "mzid") %>% pub.mzrt.naive
markings <- eicosanoidcorrelation
markings[] <- NA
diag(markings)[mark.position] <- "*"

png(height = 1560, width = 1200, pointsize = 14, res = 300, file = "cache/heatmap.png")
gplots::heatmap.2(eicosanoidcorrelation,
                  cellnote = markings,
                  notecex = 2,
                  notecol = "black",
                  distfun = dist,
                  hclustfun = hclust,
                  trace = "none",
                  col = my_palette,
                  dendrogram = "column",
                  density.info = "none",
                  labCol = NA,
                  labRow = NA,
                  asp = 1,
                  lmat = rbind(c(0, 3), c(2, 1), c(0, 4)),
                  lhei = c(2, 10, 1),
                  lwid = c(0.001, 9),
                  margins = c(0, 0),
                  key.par = list(mar=c(2.0, 2.8, 0.5, 2.8), cex=0.5),
                  keysize=0.5,
                  key.title = NA,
                  key.xlab = "")
dev.off()
```

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

## Manhattan plot

For FR02 bonferroni limit for significance is `r 0.05/length(dset$metabolite.mzids)`.

```{r manhattan descriptions}
manhattanplot <- plot.manhattanplot(filter.model(linear.lmrank),
                                    bonfp = 0.05/length(dset$metabolite.mzids),
                                    nlabels = 0)
ggsave(file = "cache/manhattanplot.png", plot = manhattanplot, height = 6, width = 7, dpi = 300)
```

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

# Multivariate models

Number of eicosanoids in three regularisation models

```{r number of eicosanoids in regularisation}
purrr::map_df(regularisation, ~as.data.frame(.x), .id="id") %>%
    dplyr::filter(grepl("mzid", term)) %>%
    group_by(id) %>%
    summarize(n = n()) %>%
    spread(id, n) %>%
    kable
```

```{r regularisation eicosanoids}
compare.eicosanoids <- mygrep(regularisation$fwdbonf$term, word = "mzid")
```

Regularisation results

```{r regularisation results}
purrr::map_df(regularisation, ~as.data.frame(.x), .id="id") %>%
    dplyr::filter(grepl("mzid", term)) %>%
    select(term, id, estimate) %>%
    spread(id, estimate) %>%
    mutate_if(is.numeric, round, 3) %>%
    kable
```

Results for Bonferroni forward selection

```{r regression bonf kable}
regularisation$fwdbonf %>%
    mutate(response = "SBP") %>%
    compare.fr02(compare.eicosanoids) %>%
    select(mzrt, name, contains("mean_ci"), contains("p.value")) %>%
    mutate_at(vars(contains("p.value")), .funs = format,  digits=3, scientific=TRUE) %>%
    arrange(mzrt) %>%
    kable
```

# Replication

```{r results for the six eicosanoids}
forest.fr02 <- compare.fr02(linear.lmrank, compare.eicosanoids)
forest.fhs <- compare.fhs()
```

```{r Comparing six-eicosanoid results}
full_join(forest.fr02, forest.fhs, by = "term", suffix = c(".fr02", ".fhs")) %>%
    rename(mzrt = mzrt.fr02, name = name.fr02) %>%
    select(mzrt, name, contains("mean_ci"), contains("p.value")) %>%
    mutate_at(vars(contains("p.value")), .funs = format,  digits=3, scientific=TRUE) %>%
    arrange(mzrt) %>%
    kable
```

```{r Plot the forest plot image}
forest.reg <- compare.fr02(regularisation$fwdbonf %>% mutate(response = "SBP"), compare.eicosanoids)
replicationforestplot(forest.fr02, forest.fhs, file = "cache/forestmap-comparison.png")
```

<img src = "cache/forestmap-comparison.png" />

# Eicosanoid profiles

## FINRISK

```{r Risk models}
rset <- getriskset(dset$data.full, regularisation)
riskmodel.fr02 <- lapply(list("unadjusted" = c(),
                              "adjusted" = c("AGE", "sex", "BMI", "curr_smk", "curr_diab")), function(adj)
                                  calculateglm(dset = rset,
                                               loops = c("fwdbonf_riskpersd", "fwdbonf_riskclass"),
                                               covariates = adj,
                                               responses = c("SBP", "HT"),
                                               filter = "risk")) %>%
    purrr::map_df(~as.data.frame(.x), .id="adjustment")
riskmodel.fhs <- riskmodel.results.fhs()
```

P for trend

```{r risk models p for trend}
calculateglm(dset = rset %>% mutate(fwdbonf_riskclass = as.numeric(fwdbonf_riskclass)),
             loops = "fwdbonf_riskclass",
             covariates = c("AGE", "sex", "BMI", "curr_smk", "curr_diab"),
             responses = c("SBP", "HT"),
             filter = "risk") %>%
    mutate(p.value = pub.p(p.value)) %>%
    select(response, term, p.value) %>%
    kable
```

Risk model results for fr02

```{r risk model results}
riskmodel.fr02 %>%
    select(adjustment, term, response, estimate, conf.low, conf.high, p.value) %>%
    kable
```

```{r risk model hypertensive per class}
rset %>% group_by(fwdbonf_riskclass) %>% summarize(htn = sum(HT == 1), n=n()) %>% kable
```

## FHS

```{r replication for FHS}
riskmodel.fhs %>%
    select(adjustment, term, response, estimate, conf.low, conf.high, p.value) %>%
    kable
```

## Comparing results

```{r Plot risk models}
plot.riskmodel(riskmodel.fhs,
               riskmodel.fr02,
               variable = "HT",
               ylab = "Odds for hypertension (95% CI)",
               file = "cache/riskforest-htn.png")
```

<img src = "cache/riskforest-htn.png" />

```{r Plot risk models SBP}
plot.riskmodel(riskmodel.fhs,
               riskmodel.fr02,
               variable = "SBP",
               ylab = "Effect size for systolic BP (95% CI)",
               file = "cache/riskforest-sbp.png")
```

<img src = "cache/riskforest-sbp.png" />

# Extra analyses in review process


## Sub-group analyses

Sub-group analyses by GFR, Aspirin, Asthma, and Age. 

```{r subgroups}
rset.extra <- 
    rset %>% mutate(gBMI = as.factor(ifelse(BMI < 30, 0, 1)),
                    gAGE = as.factor(ifelse(AGE < median(AGE), 0, 1)),
                    gM01sub = as.factor(PREVAL_ANTI_INFMEDS_M01sub),
                    asthma = as.factor(PREVAL_ASTHMA),
                    ASA = as.factor(ASA),
                    GFR = 0.993^round(AGE) *
                        case_when(sex == 1 & KREA_ENTS <= 62 ~ 144*(KREA_ENTS/61.9)^-0.329, 
                                  sex == 1 & KREA_ENTS >  62 ~ 144*(KREA_ENTS/61.9)^-1.209,
                                  sex == 0 & KREA_ENTS <= 80 ~ 141*(KREA_ENTS/79.6)^-0.411,
                                  sex == 0 & KREA_ENTS >  80 ~ 141*(KREA_ENTS/79.6)^-1.209),
                    gGFR = as.factor(ifelse(GFR < 90, 0, 1))) %>%
    filter(!is.na(gGFR), !is.na(ASA), !is.na(asthma))
```

```{r subgroups number of participants}
rset.extra %>%
    select(gBMI, gAGE, gM01sub, asthma, ASA, gGFR) %>%
    gather(key, value) %>%
    group_by(key) %>%
    summarize(N = sum(value == 1)) %>%
    kable
```

```{r subgroup models}
ret.extra <- lapply(X = c2l("gBMI", "gGFR", "ASA", "asthma", "gAGE", "gM01sub"),
                    FUN = function(df, extra, medianage) subgroupanalysis(df, extra, medianage),
                    df = rset.extra,
                    medianage = round(median(rset$AGE))) %>%
    purrr::map_df(~as.data.frame(.x), .id="subgroup")


g.extra <- plot.subgroupanalysis(ret.extra %>% filter(model != "interaction"))

ggsave(file = "cache/extra.png", plot = g.extra, height = 4, width = 5.2, dpi = 300)
```

```{r subgroup table}
ret.extra %>%
    select(subgroup, term, estimate, std.error, p.value) %>%
    kable
```

### Correlations

```{r extra correlations}
lapply(c2l("BMI", "GFR", "AGE"), function(x)
    cor.test(rset.extra$fwdbonf_riskpersd, rset.extra[[x]], method = "pearson") %>% tidy) %>%
    purrr::map_df(~as.data.frame(.x), .id="subgroup") %>%
    mutate(beta_ci = sprintf("%.2f (%.2f–%.2f) p=%.2E", estimate, conf.low, conf.high, p.value)) %>%
    kable
```

```{r extra t-test}
lapply(c2l("ASA", "ASTHMA", "gM01sub"), function(x)
    t.test(fwdbonf_riskpersd ~ group, data = rset.extra %>% rename(group := !!x)) %>% tidy) %>%
    purrr::map_df(~as.data.frame(.x), .id="subgroup") %>%
    mutate(beta_ci = sprintf("%.2f (%.2f–%.2f) p=%.2E", estimate, conf.low, conf.high, p.value)) %>%
    kable
```

### Dependency network 

The dependency network for eicosanoids associated with systolic BP.

```{r genenet}
df.mzids <- filter.model(linear.lmrank, filterterm = "SBP", plimit = 0.05) %>% pull(term)

df.net <- dset$data.full %>% dplyr::select_(.dots = df.mzids)

net.pc <- GeneNet::ggm.estimate.pcor(df.net)
net.edges <- GeneNet::network.test.edges(net.pc, direct = TRUE, fdr = TRUE)
net.network <- extract.network(net.edges, method.ggm="prob",
                               method.dir="prob",
                               cutoff.ggm = (1-0.05/length(dset$metabolite.mzids)),
                               cutoff.dir = (1-0.05/length(dset$metabolite.mzids)))

net.gr <- network.make.graph(net.network, pub.mzrt.naive(colnames(df.net)), drop.singles=TRUE)

edi <- edge.info(net.gr)

```

```{r network attributes}
getnodelist <- function(list, color = "red") {
    nodeAttrs <- list()
    temp <- rep(color, length(list))
    names(temp) <- list
    nodeAttrs$fillcolor <- temp
    nodeAttrs
}

nodeAttrs <- getnodelist(pub.mzrt.naive(forest.reg$term))


globalAttrs = list()
#globalAttrs$edge = list(color = "black", lty = "solid", lwd = 1, arrowsize=1)
globalAttrs$node = list(fillcolor = "gray", shape = "ellipse")

edgeAttrs <- list()
edgeAttrs$dir <- edi$dir
edgeAttrs$lty <- ifelse(edi$weight < 0, "dotted", "solid") # negative correlation -> dotted
edgeAttrs$color <- ifelse(edi$dir == "none", "black", "red")
edgeAttrs$fontsize <- 20*edi$weight/edi$weight
edgeAttrs$label <- round(edi$weight, 2) # use partial correlation as edge labels
```

```{r plot network}
png(file = "myplot.png", width = 2000, height = 2000, pointsize = 12, units = "px")
plot(net.gr, attrs = globalAttrs, edgeAttrs = edgeAttrs, nodeAttrs = nodeAttrs, "fdp")
dev.off()
```

# Supplement

## Linear and logistic regression results in FR02

```{r linear results for supplement}
sup.linres <- pub.lmrank(linear.lmrank) %>%
    mutate(MZ = getmz(term),
           RT = getrt(term)) %>%
    select(-term) %>%
    select(MZ, RT, everything()) 

typology.tbls2 <- data.frame(
  col_keys = colnames(sup.linres),
  what = c("", "",
           rep("Systolic BP", 2), 
           rep("Diastolic BP", 2), 
           rep("MAP", 2), 
           rep("Pulse pressure", 2), 
           rep("Hypertension", 2)),
  measure = c("MZ", "RT", rep(c("β±SE", "p"), 5)),
  stringsAsFactors = FALSE) 

typologyformatter(data = sup.linres, font = 9, typology = typology.tbls2) %>%
        colformat_num(., col_keys = c("MZ", "RT"), digits = 4)

```

```{r other covariates}
sup.riskmodel <- calculateglm(dset = rset,
                              loops = "fwdbonf_riskpersd",
                              covariates = c("AGE", "sex", "BMI", "curr_smk", "curr_diab"),
                              responses = c("SBP", "HT"),
                              filter = ".") %>%
    mutate(beta_ci = sprintf("%.2f (%.2f–%.2f)_p=%.2E", estimate, conf.low, conf.high, p.value)) %>%
    select(response, term, beta_ci)
           
sup.riskmodel %>%
    spread(response, beta_ci) %>%
    separate(HT, c("A", "B"), sep = "_") %>%
    separate(SBP, c("C", "D"), sep = "_") %>%
    typologyformatter
```

```{r Write docx to file, results = FALSE, include = FALSE}
save.image(file = paste0("rds/session-", now, ".Rdata"))
```