For publication

jjmpal · Mar 26, 2020 · 09ef836 · 09ef836
1 parent e3d9981
commit 09ef836
Show file tree

Hide file tree

Showing 5 changed files with 138 additions and 240 deletions.
diff --git a/articleone-functions.R b/articleone-functions.R
@@ -24,12 +24,20 @@ mygrep <- function(..., word, ignorecase = TRUE, complement = FALSE) {
 }
 
 
-pub.mzrt.naive <- function(mzid, formatter = "%.4f / %.2f") {
+pub.mzrt.naive <- function(mzid, formatter = "%.4f / %.4f") {
     sprintf(formatter,
             as.numeric(gsub("_.*", "", gsub("mzid_", "", mzid))),
             as.numeric(gsub(".*_", "", mzid)))
 }
 
+mzrt.standardize <- function(mzid, formatter = "mzid_%.6f_%.4f") {
+    lapply(mzid, function(id) {
+        masscharge <- strsplit(id, "_")[[1]][1] %>% as.numeric
+        retentiontime <- strsplit(id, "_")[[1]][2] %>% as.numeric
+        sprintf(formatter, masscharge, retentiontime)
+    }) %>% unlist
+}
+
 pub.mzrt <- function(mzid,
                      metabolitenames,
                      mark = "",
@@ -269,8 +277,10 @@ stepwise.bonfp <- function(full.model,
 
 #' Formats p values
 pub.p <- function(p) {
-  p <- as.numeric(p)
-  ifelse(p < 0.01, ifelse(p<0.001, "p<0.001", sprintf("%.3f", p)), sprintf("%.2f", p))
+    p <- as.numeric(p)
+    case_when(p<0.001 ~ "<0.001",
+              p < 0.01 ~ sprintf("%.3f", p),
+              TRUE ~ sprintf("%.2f", p))
 }
 
 ctolist <- function(c) {
@@ -334,12 +344,12 @@ comparecohorts <- function(fr02, fhs, mapping, eicosanoids) {
 pub.lmrank <- function(...,
                        by="term",
                        mark = "",
-                       models = c2l("MAP", "SBP", "DBP", "PP", "HT"),
+                       models = c2l("SBP", "DBP", "MAP", "PP", "HT"),
                        arrangebyp = FALSE) {
   linrunall <- dplyr::bind_rows(...) %>%
     dplyr::filter(grepl("mzid", term), qval < 0.05) %>%
     dplyr::mutate("betase" = sprintf("%.2f±%.2f", estimate, std.error),
-                  p = pub.p(p.value)) %>%
+                  p = sprintf("%0.0e", p.value)) %>%
     dplyr::select(response, term, betase, p)
 
   linrunarray <- lapply(models,
@@ -351,3 +361,15 @@ pub.lmrank <- function(...,
     Reduce(function(dtf1, dtf2) dplyr::full_join(dtf1, dtf2, by=by, suffix=c(".1", ".2")), .) %>%
     dplyr::arrange(term)
 }
+
+bionames <- function(x) {
+    name <- gsub("Novel EIC_", "Novel-", x) %>%
+        gsub("EIC_", "Putative-", .) %>%
+        gsub("Eicosanoid_", "", .) %>%
+        gsub("FFA", "", .) %>%
+        gsub(" *\\[M-H\\]", "", .) %>%
+        gsub("[ _]", "", .) %>%
+        gsub(" ", "-", .) %>%
+        gsub(";", "; ", .) %>%
+        gsub("alpha", "α", .)
+}
diff --git a/articleone-plots.R b/articleone-plots.R
@@ -1,12 +1,3 @@
-mycolors <- function(i) {
-  colors <- list("red" = "#FF7575",
-                 "blue" = "#4B81B6",
-                 "gray" = "#bebebe",
-                 "white" = "#FFFAFA",
-                 "purewhite" = "#FFFFFF")
-  get(i, colors)
-}
-
 mzrt2mz <- function(mzid) {
   mzid %>%
     gsub("mzid_", "", .) %>%
@@ -18,11 +9,12 @@ plot.manhattanplot <- function(linres,
                                bonfp = NULL,
                                response = "SBP",
                                nlabels = 10,
+                               pointsize = 1.0,
                                xlab = "Eicosanoids Rank Ordered by Mass to Charge Ratio",
                                ylab = "Negative Log P Value",
-                               color_scheme = c("positive" = mycolors("red"),
-                                                "negative" = mycolors("blue"),
-                                                "insignificant" = mycolors("gray"))) {
+                               color_scheme = c("positive" = "red",
+                                                "negative" = "blue",
+                                                "insignificant" = "gray")) {
     dat <- linres %>%
     dplyr::filter(response == "SBP", grepl('^mzid', term)) %>%
     dplyr::arrange(term) %>%
@@ -39,7 +31,7 @@ plot.manhattanplot <- function(linres,
   ggplot2::ggplot(dat, ggplot2::aes(x = mz, y = neg_log10_pvalue)) +
     ggplot2::ylab(ylab) +
     ggplot2::xlab(xlab) +
-    ggplot2::geom_point(ggplot2::aes(color = coloring), size = 1.5) +
+    ggplot2::geom_point(ggplot2::aes(color = coloring), size = pointsize) +
     ggrepel::geom_text_repel(data = label_dat,
                              ggplot2::aes(label = mzrt),
                              segment.size = 0.1,
@@ -59,17 +51,17 @@ plot.manhattanplot <- function(linres,
 }
 
 spearmancorrelation  <- function(data, mzids) {
-  dat <- data$data.full %>% dplyr::select(mzids)
-  colnames(dat) <- colnames(dat) %>% pub.mzrt(., data$metabolite.names)
-  cor <- cor(dat, method = 'spearman')
+    dat <- data$data.full %>% dplyr::select(mzids)
+    colnames(dat) <- colnames(dat) %>% pub.mzrt.naive
+    cor <- cor(dat, method = 'spearman')
 }
 
 replicationforestplot <- function(forest.fr02, forest.fhs, file) {
-    png(width = 1600, height = 700, res = 150, file = file)
+    png(width = 1600, height = 780, res = 150, file = file)
     forestplot::forestplot(
-                    labeltext = cbind(c("Eicosanoid", "", forest.fhs$name),
-                                      c("FINRISK\n β (95% CI)", "", forest.fr02$mean_ci),
-                                      c("FHS\n β (95% CI)", "", forest.fhs$mean_ci)),
+                    labeltext = cbind(c("Eicosanoid", NA, forest.fhs$name),
+                                      c("FINRISK\n β (95% CI)", NA, forest.fr02$mean_ci),
+                                      c("FHS\n β (95% CI)", NA, forest.fhs$mean_ci)),
                     mean = cbind(c(NA, NA, forest.fr02$estimate),
                                  c(NA, NA, forest.fhs$estimate)),
                     lower = cbind(c(NA, NA, forest.fr02$conf.low),
@@ -78,7 +70,7 @@ replicationforestplot <- function(forest.fr02, forest.fhs, file) {
                                   c(NA, NA, forest.fhs$conf.high)),
                     legend = c("FINRISK", "FHS"),
                     align = c("l", "l", "l"),
-                    graph.pos = 4,
+                    graph.pos = 3,
                     title = "",
                     xlog = FALSE,
                     xlab = " β (95%-CI)",
@@ -88,11 +80,11 @@ replicationforestplot <- function(forest.fr02, forest.fhs, file) {
                                      title = gpar(cex = 1.2)),
                     xticks = seq(-1, 4),
                     clip =c(-1, 4),
-                    col = fpColors(box=c("blue", "darkred")),
+                    col = fpColors(box=c("#0433ff", "#ff2500")),
                     zero = 0, 
-                    lineheight = unit(14, "mm"), 
-                    boxsize = 0.2, 
-                    colgap = unit(8, "mm"),
+                    lineheight = unit(16, "mm"),
+                    boxsize = 0.15,
+                    colgap = unit(4, "mm"),
                     lwd.ci = 1)
     dev.off()
 }
@@ -115,7 +107,7 @@ plot.riskmodel <- function(riskmodel.fhs,
                                           y = estimate,
                                           ymin = conf.low,
                                           ymax = conf.high,
-                                          col = cohort)) +
+                                          col = factor(cohort, c("FINRISK", "FHS")))) +
         facet_wrap(~adjustment,
                    strip.position = "top", 
                    labeller = as_labeller(c("unadjusted" = "Unadjusted",

diff --git a/articleone-replication.R b/articleone-replication.R
@@ -3,9 +3,9 @@ compare.names <- function() {
             "mzid_339.179900_1.7205", "11-dehydro-2,3-dinor-TXB2",
             "mzid_279.196600_3.7247", "12-HHTrE",
             "mzid_295.227900_4.8902", "295.2279/4.89 (putative eicosanoid)",
-            "mzid_319.228000_5.6733", "5,6-EET",
-            "mzid_331.267900_6.5274", "Adrenic Acid",
-            "mzid_265.181000_3.5705", "Tetranor-12(R)-HETE")
+            "mzid_319.228000_5.6733", "319.2280/5.67 (unknown)",
+            "mzid_331.267900_6.5274", "Adrenic acid",
+            "mzid_265.181000_3.5705", "265.1810/3.57 (putative eicosanoid)")
 }
 
 compare.fhs <- function() {
@@ -23,7 +23,7 @@ compare.fhs <- function() {
         arrange(name)
 }
 
-riskmodel.fhs <- function() {
+riskmodel.results.fhs <- function() {
     tribble(
         ~adjustment, ~response, ~term, ~estimate, ~std.error, ~statistic, ~p.value, ~conf.low, ~conf.high,
         "adjusted", "HT", "fwdbonf_riskclass2", 1.267338843, 0.117173977, 2.021944707, 0.043182064, 1.007414265, 1.59493293,
@@ -44,6 +44,7 @@ riskmodel.fhs <- function() {
         "unadjusted", "SBP", "fwdbonf_riskpersd", 2.721582913, 0.318103303, 8.555657498, 1.88E-17, 2.098111895, 3.345053931)
 }
 
+
 compare.fr02 <- function(lmrank, eicosanoids) {
     lmrank %>%
         filter(term %in% eicosanoids, response == "SBP") %>%