In krisrs1128/treelapse: Visualization for Tree-Structured Collections of Counts and Series

Often, it becomes easier to interpret the results of different confirmatory analysis when the results can be studied interactively. Here, we take an example hierarchical testing analysis from [@callahan2016bioconductor] and study the output interactively using the parallel coordinates version of timebox trees.

library("knitr")
library("reshape2")
library("plyr")
library("dplyr")
library("phyloseq")
library("DESeq2")
library("structSSI")
library("multtest")
library("treelapse")
opts_chunk$set(fig.width = 8, fig.height = 8, fig.align = "center")
data(mouse_diet)
ps_dds <- phyloseq_to_deseq2(mouse_diet,  ~ age_binned + family_relationship)

# geometric mean, set to zero when all coordinates are zero
geo_mean_protected <- function(x) {
  if (all(x == 0)) {
    return (0)
  }
  exp(mean(log(x[x != 0])))
}

geoMeans <- apply(counts(ps_dds), 1, geo_mean_protected)
ps_dds <- estimateSizeFactors(ps_dds, geoMeans = geoMeans)
ps_dds <- estimateDispersions(ps_dds)
abund <- getVarianceStabilizedData(ps_dds)

short_names <- substr(rownames(abund), 1, 5)%>%
  make.names(unique = TRUE)
rownames(abund) <- short_names

el <- phy_tree(mouse_diet)$edge
el0 <- el
el0 <- el0[nrow(el):1, ]
el_names <- c(short_names, seq_len(phy_tree(mouse_diet)$Nnode))
el[, 1] <- el_names[el0[, 1]]
el[, 2] <- el_names[as.numeric(el0[, 2])]
unadj_p <- treePValues(el, abund, sample_data(mouse_diet)$age_binned)

hfdr_res <- hFDR.adjust(unadj_p, el, .75)
summary(hfdr_res)

adj_p <- [e-mail]
adj_p$id <- rownames(adj_p)
adj_p$adjp[is.na(adj_p$adjp)] <- 1
adj_p <- cbind(adj_p, mt.rawp2adjp(adj_p$unadjp))

values <- adj_p %>%
  select(starts_with("adjp"), id) %>%
  dplyr::rename(structSSI = adjp) %>%
  melt(id.vars = "id",) %>%
  dplyr::rename(unit = id, time = variable) %>%
  mutate(
    value = - log(value),
    time = gsub("adjp.", "", time)
  )
values$time <- values$time %>%
  revalue(
    replace = c(
      "Bonferroni" = "Bonf.",
      "Hochberg" = "Hoch.",
      "SidakSS" = "SidSS",
      "SidakSD" = "SidSD"
    )
  )

colnames(el) <- c("parent", "child")
timebox_tree(
  values,
  el,
  height = 350,
  width = 450,
  display_opts = list(
    "axis_font_size" = 11,
    "n_ticks_y" = 2
  )
)

tax <- tax_table(mouse_diet) %>%
  data.frame()
tax$seq <- short_names

[e-mail]$seq <- rownames([e-mail])
tax %>%
  left_join([e-mail]) %>%
  arrange(adjp) %>% head(10)

interest_taxa <- c(
  "GCAAG.176",
  "GCAAG.52",
  "GCAAG.253",
  "GCAAG.81"
)

tax %>%
  filter(seq %in% interest_taxa)

interest_taxa <- c(
  "GCAAG.139",
  "GCAAG.212",
  "GCAAG.164",
  "GCAAG.47"
)
tax %>%
  filter(seq %in% interest_taxa)

krisrs1128/treelapse documentation built on Jan. 3, 2020, 6:29 p.m.