In mvazquezs/chmitools: Package to manage & obtain results of the CHMI Study

knitr::opts_chunk$set(comment = NA,
  fig.path = '01_figures_aim4_q01_q04_univar/',
  fig.width = 20,
  fig.height = 12,
  results = 'markup',
  tidy = FALSE,
  message = FALSE,
  warning = FALSE,
  echo = FALSE)

alpha <- 0.05

r packages

library(devtools)
library(prettydoc)

# avoid 'objects' in workspace
rm(list = ls())

setwd(Sys.getenv('HOME'))

if(Sys.info()[["nodename"]] == 'mbp-de-csc.csc.es'){

load_all('~/Documents/git/chmitools')

} else if(Sys.info()[["nodename"]] != 'mbp-de-csc.csc.es'){

load_all('~/git/chmitools') 

}

load data

Main bullet points in load data section. chmi.phen() load data by group and aim. filter data by t_igg, antigen_ag and t_point_ag/t2_point_ag columns.

# load data
dat <- chmi.phen(
  study_chmi = 'antibodies',
  aim_ab_chmi = 'aim_4') %>%
  mutate_at(vars(c('original_id', 'n_study_ag')), as.factor)

dat_f <- chmi.phen(
  study_chmi = 'antibodies',
  aim_ab_chmi = 'aim_4',
  fold_change = TRUE) %>%
  mutate_at(vars(c('original_id', 'n_study_ag')), as.factor)

# list to perform different functions
  n_igg <- sort(unique(dat$t_igg))
  l_t2point_ag_1 <- c('C-1', 'D7', 'D11_D13')
  l_t2point_ag_2 <- c('D7')
  l_group <- c('t_igg', 'antigen_ag', 't2_point_ag')
  l_varsy <- c('pcr_pos_pid', 'tbs_pos_qpcr')
  l_varsx <- c('log10_mfi', 'log10_rmfi')
  l_dbl_or <- c('or', 'ci_lower', 'ci_upper')
  l_pval <- c('raw_pval', 'adjust_pval')
  l_showed <- c(l_group, 'vars_y', 'vars_x', l_dbl_or, l_pval)

# subsets for questions
dat_1 <- dat %>%
  filter(t2_point_ag %in% l_t2point_ag_1) %>%
  filter(mal_exposure != 'placebo') %>%
  filter(complete.cases(log10_mfi)) %>%
  droplevels() ## Q1 & Q3

dat_2 <- dat_f %>%
    filter(t2_point_ag %in% l_t2point_ag_2) %>%
  filter(mal_exposure != 'placebo') %>%
    filter(!is.na(log10_rmfi)) %>%
    droplevels() ## Q2 & Q4

tables: odd ratios & p_values in glm models

Odd ratios, lower & upper confidence interval, and adjusted p_values. chmi.stat.glm_or() obtain all parameters for all groups. 1st subset (dat_1): variable (pcr_pos_pid), t2_point_ag ('C-1', 'D7', 'D11_D13'). 2nd subset (dat_2): variable (pcr_pos_pid), t2_point_ag ('D7 / C-1'). 3rd subset (dat_1): variable (tbs_pos_qpcr), t2_point_ag ('C-1', 'D7'). 4th subset (dat_2): variable (tbs_pos_qpcr), t2_point_ag ('D7 / C-1'). used p.adjust() function to adjust p-values by Benjamini & Hochberg method.

# glm models for `l_varsy` & `log10_mfi`
  tab_1 <- chmi.stat.glm_or(
        phen = dat_1,
        l_group = l_group,
        l_varsy = l_varsy,
        l_varsx = l_varsx[1],
        arg_multivar = 'univariate',
        arg_info_model = NULL,
      arg_info_beta = NULL,
        format_pval = TRUE,
        padj_method = 'BH') %>%
    mutate(varsy_t2_point_ag = factor(paste0(vars_y, '_', t2_point_ag))) %>%
    filter(varsy_t2_point_ag != 'tbs_pos_qpcr_D11_D13') %>%
    select(-varsy_t2_point_ag)

# glm models for `l_varsy` & `log10_rmfi`
  tab_2 <- chmi.stat.glm_or(
        phen = dat_2,
        l_group = l_group,
        l_varsy = l_varsy,
        l_varsx = l_varsx[2],
        arg_multivar = 'univariate',
      arg_info_model = NULL,
      arg_info_beta = NULL,
        format_pval = TRUE,
        padj_method = 'BH')


# join `tab_1` & `tab_2`
    tab <- bind_rows(tab_1, tab_2) %>%
        mutate_at(vars(c('t2_point_ag')),  ~ fct_relevel(., c('C-1', 'D7', 'D11_D13'))) %>%
        mutate_at(vars(c('vars_x', 'adjust_text', 'adjust_signif')), as.factor)

tables for q01 & q03

# 'tab_1' datatable
  datatable(tab_1 %>% select(l_showed),
    class = 'cell-border stripe',
    filter = 'top',
    extensions = c('Buttons', 'FixedColumns'),
      options = list(dom = 'Bfrtip',
        buttons = c('csv', 'copy'),
        pageLength = 10,
        scrollX = TRUE,
        text = 'Download')) %>%
    formatRound(l_dbl_or, 3) %>%
    formatRound(l_pval, 5)

 ```

## tables for `q02` & `q04`

```r

# 'tab_2' datatable
  datatable(tab_2 %>% select(l_showed),
    class = 'cell-border stripe',
    filter = 'top',
    extensions = c('Buttons', 'FixedColumns'),
      options = list(dom = 'Bfrtip',
        buttons = c('csv', 'copy'),
        pageLength = 10,
        scrollX = TRUE,
        text = 'Download')) %>%
    formatRound(l_dbl_or, 3) %>%
    formatRound(l_pval, 5)

plots for odd_ratio, raw_pval & adjust_pval

# heatmaps plots
h_odd <- chmi.stat_plot.or_heatmap(
    phen = tab,
    arg_plot = 'or')

h_raw <- chmi.stat_plot.or_heatmap(
    phen = tab,
    arg_plot = 'raw_pval')

h_adj <- chmi.stat_plot.or_heatmap(
    phen = tab,
    arg_plot = 'adjust_pval')

# pointrange plots
p_odd <- chmi.stat_plot.pointrange_or(
    phen = tab,
    arg_igg = n_igg)

heatmaps for odd ratios

# heatmap
print(h_odd)

heatmap for raw p_values

# loop
print(h_raw)

heatmap for adjusted p_values

# loop
print(h_adj)

point_range for odd ratios

# point_range
print(p_odd)

save in r_data folder

# save

if(Sys.info()[["nodename"]] == 'mbp-de-csc.csc.es'){

save.image(file = '~/Documents/git/chmitools/projects/04_aim4/r_data/01_aim4_q1_q4_univar.RData')

} else if(Sys.info()[["nodename"]] != 'mbp-de-csc.csc.es'){

save.image(file = '~/git/chmitools/projects/04_aim4/r_data/01_aim4_q1_q4_univar.RData')

}

mvazquezs/chmitools documentation built on May 1, 2020, 2:06 a.m.