R/chmi.stat_plots.R

Defines functions chmi.stat_plot.or_heatmap chmi.stat_plot.median_line chmi.stat_plot.pval_boxfold chmi.stat_plot.pval_boxplot

#---------------------------------------------
# Statistic functions to plot in `CHMI` Study
#---------------------------------------------


### p_values & boxplot for 'log10_mfi'--------------------------------------------------------------


# '@export
chmi.stat_plot.pval_boxplot <- function(
  phen, pval_phen,
  l_igg,
  var_x, var_y,
  aes_color,
  aes_geom = c('double', 'variable'))
## `variable` have two or more levels
## `double` only has one level
{
# enquos
  var_x <- enquo(var_x)
  var_y <- enquo(var_y)
  aes_color <- enquo(aes_color)

# list
l_plots <- list()

	# loop
    for(i in seq_along(l_igg)) {
      # plot
      p_1 <- ggplot(
          data = phen %>% filter(t_igg == l_igg[i]),
            aes(x = !!var_x, y = !!var_y, color = !!aes_color)) +
        geom_boxplot(aes(x = factor(!!var_x), y = !!var_y), outlier.shape = NA) +
        geom_jitter(aes(x = factor(!!var_x), y = !!var_y), width = .2, height = 0, alpha = .25) +
        facet_wrap(~ antigen_ag) +
        scale_color_viridis(discrete = TRUE, begin = 0, end = .5, option = 'D') +
        labs(title = paste0(l_igg[i]), x = paste0(quo_name(var_x)), y = paste0(quo_name(var_y)))

				## argument `aes_geom`
				if(aes_geom == 'variable') {
        p_1 <- p_1 +
         geom_text(
          data = pval_phen %>% filter(t_igg == l_igg[i]),
            aes(x = !!var_x, y = max(phen[, quo_name(var_y)] + 1), label = adjust_text),
            family = 'Georgia',  colour = 'black', size = 4)

				} else {
      	p_1 <- p_1 +
       	geom_text(
        	data = pval_phen %>% filter(t_igg == l_igg[i]),
          	aes(x = 2, y = max(phen[, quo_name(var_y)] + 1), label = adjust_text),
          family = 'Georgia',  colour = 'black', size = 4)
				}
		# list
   	l_plots[[length(l_plots) + 1]] = p_1
	  }

# return
  return(l_plots)
}



### boxplot & p_values for 'mfi' in 'fold changes'----------------------------------------

# '@export
chmi.stat_plot.pval_boxfold <- function(
  phen, pval_phen,
  l_group,
  l_igg,
  var_y, aes_color)
{
# enquos
  var_y <- enquo(var_y)
  aes_color <- enquo(aes_color)

# list
  l_plots <- list()

  # loop
  for(i in seq_along(l_igg)) {
    # plot
    p_1 <- ggplot(
        data = phen %>% filter(t_igg == l_igg[i]),
          aes(x = t2_point_ag, y = !!var_y, color = !!aes_color)) +
      geom_boxplot(outlier.shape = NA) +
      geom_point(
        position = position_jitterdodge(jitter.width = .2, jitter.height = 0),
        alpha = .25) +
      facet_wrap(~ antigen_ag) +
      scale_color_viridis(discrete = TRUE, begin = 0, end = .5, option = 'D') +
      geom_text(data = pval_phen %>% filter(t_igg == l_igg[i]),
        aes(x = t2_point_ag, y = p_ymax, label = adjust_text),
        family = 'Georgia',  colour = 'black', size = 3) +
    labs(title = paste0(l_igg[i]), x = 't2_point', y = paste0(quo_name(var_y)))

  # list
   	l_plots[[length(l_plots) + 1]] = p_1
    }

# return
  return(l_plots)
}


### medians & geom_line for 'mfi'-----------------------------------------------------------

# '@export
chmi.stat_plot.median_line <- function(
  phen,
  l_igg,
  var_x, var_y,
  aes_color)
{
# enquos
  var_x <- enquo(var_x)
  var_y <- enquo(var_y)
  aes_color <- enquo(aes_color)

  # list
    l_plots <- list()
    # loop
    for(i in seq_along(l_igg)) {
      # plot
      p_1 <- ggplot(
          data = phen %>% filter(t_igg == l_igg[[i]]),
            aes(x = !!var_x, y = !!var_y, color = !!aes_color, group = original_id)) +
        geom_line(alpha = .2) +
        stat_summary(fun.y = 'median',
          aes(group = !!aes_color), geom = 'point', size = 2, alpha = .7) +
        facet_wrap(~ antigen_ag) +
        scale_color_viridis(discrete = T, begin = 0, end = .5, option = 'D')

    # list
      l_plots[[length(l_plots) + 1]] = p_1
      }

# return
  return(l_plots)
}



### or plots ---------------------------------------------------------------------------------------

# '@export
chmi.stat_plot.or_heatmap <- function(
	phen,
	l_varsy, l_varsx,
  arg_igg = c('IgG', 'IgG1', 'IgG2', 'IgG3', 'IgG4', 'IgM'),
	arg_plot = c('or', 'raw_pval', 'adjust_pval', 'or_raw', 'or_adj'),
	stat = c('ggpubr'))
{
# stat
  stat <- match.arg(stat, c('ggpubr'))

### variables to avoid & updates
l_vars <- c('null_deviance', 'df_null', 'log_lik', 'aic', 'bic', 'deviance', 'df_residual',
	'r_squared', 'adj_r_squared', 'sigma', 'mod_statistic', 'mod_pval', 'df', 'deviance', 'df_residual',
	'std_error', 't_value')


  # fix round variables
  phen[, c('or')] <- round(phen[, c('or')], 2)
  phen[, c('raw_pval', 'adjust_pval')] <- round(phen[, c('raw_pval', 'adjust_pval')], 5)

### argument to filter `l_varsy`
if(!missing(l_varsy)) {
# filter 'phen'
	phen <- phen %>%
		filter_at(vars(vars_y), all_vars(. %in% l_varsy)) %>%
    droplevels()
}

### argument to filter `l_varsx`
if(!missing(l_varsx)) {
# filter 'phen'
	phen <- phen %>%
		filter_at(vars(vars_x), all_vars(. %in% l_varsx)) %>%
    droplevels()
}

### argument to filter `igg`
if(!missing(arg_igg)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(t_igg), all_vars(. %in% arg_igg)) %>%
    droplevels()
}

# group_split
phen_1 <- phen %>%
  select(-one_of(l_vars)) %>%
  filter(complete.cases(.)) %>%
  droplevels() %>%
	group_by(vars_y, vars_x) %>%
	group_split()

  # list
  l_or <- list()
  l_raw <- list()
  l_adj <- list()
  l_orraw <- list()
  l_oradj <- list()

### make all plots
	# loop
  for(i in seq_along(phen_1)) {
	# plot 'or'
 		p_1 <- ggplot(
 				phen_1[[i]],
      	aes(x = t2_point_ag, y = antigen_ag, fill = or)) +
      geom_tile() +
      scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
        breaks = seq(
          min(phen_1[[i]] %>% select(or)),
          max(phen_1[[i]] %>% select(or)), by = 1),
        limits = c(
          min(phen_1[[i]] %>% select(or)),
          max(phen_1[[i]] %>% select(or)))) +
      facet_wrap( ~ t_igg, nrow = 1, scales = 'free_x') +
      labs(title = unique(paste0(
		   			droplevels(phen_1[[i]][['vars_y']]), ' ~ ',
		   			droplevels(phen_1[[i]][['vars_x']]))),
		   		x = 't_point', y = 'antigen', fill = 'odd_ratio')

   # plot 'raw_pval'
   	 p_2 <- ggplot(
 				phen_1[[i]],
      	aes(x = t2_point_ag, y = antigen_ag, fill = raw_pval)) +
      geom_tile() +
      scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
        breaks = seq(
          min(phen_1[[i]] %>% select(raw_pval)),
          max(phen_1[[i]] %>% select(raw_pval)), by = .5),
        limits = c(
          min(phen_1[[i]] %>% select(raw_pval)),
          max(phen_1[[i]] %>% select(raw_pval)))) +
      facet_wrap( ~ t_igg, nrow = 1, scales = 'free_x') +
      labs(title = unique(paste0(
		   			droplevels(phen_1[[i]][['vars_y']]), ' ~ ',
		   			droplevels(phen_1[[i]][['vars_x']]))),
		   		x = 't_point', y = 'antigen', fill = 'raw_pval')

	 # plot 'adj_pval'
   	 p_3 <- ggplot(
 				phen_1[[i]],
      	aes(x = t2_point_ag, y = antigen_ag, fill = adjust_pval)) +
      geom_tile() +
      scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
        breaks = seq(
          min(phen_1[[i]] %>% select(adjust_pval)),
          max(phen_1[[i]] %>% select(adjust_pval)), by = .5),
        limits = c(
          min(phen_1[[i]] %>% select(adjust_pval)),
          max(phen_1[[i]] %>% select(adjust_pval)))) +
      facet_wrap( ~ t_igg, nrow = 1, scales = 'free_x') +
      labs(title = unique(paste0(
		   			droplevels(phen_1[[i]][['vars_y']]), ' ~ ',
		   			droplevels(phen_1[[i]][['vars_x']]))),
		   		x = 't_point', y = 'antigen', fill = 'adj_pval')

  # combine list
    p_raw <- ggpubr::ggarrange(p_1, p_2,
    	ncol = 1, nrow = 2,
			common.legend = FALSE,
			legend = 'right',
			labels = c('a. odd_ratio', 'b. raw_pvalues'))

		p_adj <- ggpubr::ggarrange(p_1, p_3,
			ncol = 1, nrow = 2,
			common.legend = FALSE,
			legend = 'right',
			labels = c('a. odd_ratio', 'b. adjust_pvalues'))

	# list
  	l_or[[length(l_or) + 1]] = p_1
  	l_raw[[length(l_raw) + 1]] = p_2
  	l_adj[[length(l_adj) + 1]] = p_3
  	l_orraw[[length(l_orraw) + 1]] = p_raw
  	l_oradj[[length(l_oradj) + 1]] = p_adj
 }


### argument 'plots'
	if(arg_plot == 'or') {
 	# return
  return(l_or)

	} else if(arg_plot == 'raw_pval') {
  # return
  return(l_raw)

	} else if(arg_plot == 'adjust_pval') {
	# return
	return(l_adj)

	} else if (arg_plot == 'or_raw') {
	# return
	return(l_orraw)

	} else {
	# return
	return(l_oradj)
	}
}


### beta plots ---------------------------------------------------------------------------------------

# '@export
chmi.stat_plot.heatmap <- function(
	phen,
  l_group,
	l_varsy = c('NULL'),
  l_varsx = c('NULL'),
  arg_igg = c('IgG', 'IgG1', 'IgG2', 'IgG3', 'IgG4', 'IgM', 'NULL'),
	arg_plot = c('beta', 'or'))
{
### list of variables to choose
  l_vars_arg <- c(l_group, 'vars_y', 'vars_x', arg_plot)

### update 'data'
  phen_1 <- phen %>%
    select(one_of(l_vars_arg)) %>%
    filter(complete.cases(.)) %>%
    droplevels()


### argument to filter `l_varsy`, `l_vars_x`
  # list
  l_vars_1 <- c(l_varsy[2], l_varsx[1], arg_igg)
  l_vars_2 <- c(l_varsy[1], l_varsx[1])
	l_vars_3 <- c(l_varsy[1])
	l_vars_4 <- NULL


  ### conditions
  if(length(l_vars_1) != 0) {
   phen_1 <- phen_1 %>% filter_at(vars(t_igg, vars_y, vars_x), all_vars(. %in% l_vars_1))

	} else {
	 phen_1

  }

### group_split
  phen_1 <- phen_1
	  group_by_at(l_names[c(2:3)]) %>%
	  group_split()

  # list
  l_plots <- list()

### make all plots
	# loop
  for(i in seq_along(phen_1)) {
	# plot 'or'
 		p_1 <- ggplot(
 				phen_1[[1]],
      	aes(x = t2_point_ag, y = antigen_ag, fill = beta)) +
      geom_tile() +
      scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
        breaks = seq(
          min(phen_1[[1]] %>% select(beta)),
          max(phen_1[[1]] %>% select(beta)), by = 1),
        limits = c(
          min(phen_1[[1]] %>% select(beta)),
          max(phen_1[[1]] %>% select(beta)))) +
      facet_wrap( ~ t_igg, nrow = 1, scales = 'free_x') +
      labs(title = unique(paste0(
		   			droplevels(phen_1[[1]][['vars_y']]), ' ~ ',
		   			droplevels(phen_1[[1]][['vars_x']]))),
		   		x = 't_point', y = 'antigen', fill = paste0(arg_plot))

  # list
  	l_plots[[length(l_plots) + 1]] = p_1
  }

### return
  return(l_plots)
}



### geom_point plots (for OR)-------------------------------------------------------------------------------

# '@export
chmi.stat_plot.pointrange_or <- function(
	phen,
	l_varsy, l_varsx,
	arg_igg = c('IgG', 'IgG1', 'IgG2', 'IgG3', 'IgG4', 'IgM'),
	stat = c('ggpubr'))
{
# stat
  stat <- match.arg(stat, c('ggpubr'))

# variables to avoid
	l_vars <- c('null_deviance', 'df_null', 'log_lik', 'aic', 'bic', 'deviance', 'df_residual',
		'r_squared', 'adj_r_squared', 'sigma', 'mod_statistic', 'mod_pval', 'df', 'deviance', 'df_residual',
		'std_error', 't_value')

	# fix round variables
	phen[, c('or', 'ci_lower', 'ci_upper')] <- round(phen[, c('or', 'ci_lower', 'ci_upper')], 2)
	phen[, c('raw_pval', 'adjust_pval')] <- round(phen[, c('raw_pval', 'adjust_pval')], 5)


# argument to filter `l_varsy` | `l_varsx`
if(!missing(l_varsy)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(vars_y), all_vars(. %in% l_varsy)) %>%
    droplevels()
}

if(!missing(l_varsx)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(vars_x), all_vars(. %in% l_varsx)) %>%
    droplevels()
}

if(!missing(arg_igg)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(t_igg), all_vars(. %in% arg_igg)) %>%
    droplevels()
}


# group_split
phen_1 <- phen %>%
  select(-one_of(l_vars)) %>%
  filter(complete.cases(.)) %>%
  droplevels() %>%
	group_by(vars_y, vars_x, t2_point_ag) %>%
	group_split()

# list
	l_plots <- list()

	### make all plots
		# loop
		for(i in seq_along(phen_1)) {
		# plot 'point_range'
	 		p_1 <- ggplot(
	 				phen_1[[i]],
		    	aes(x = antigen_ag, y = or)) +
		    geom_linerange(aes(ymin = ci_lower, ymax = ci_upper, colour = or), size = 0.2) +
		    geom_pointrange(aes(ymin = ci_lower, ymax = ci_upper, colour = or), size = 0.2) +
		    scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
		      breaks = seq(
		      	 min(phen_1[[i]] %>% select(or)),
		      	 max(phen_1[[i]] %>% select(or)), by = 1),
		      limits = c(
        	 	 min(phen_1[[i]] %>% select(or)),
		      	 max(phen_1[[i]] %>% select(or)))) +
		   labs(title = unique(paste0(
		   			droplevels(phen_1[[i]][['vars_y']]), ' ~ ',
		   			droplevels(phen_1[[i]][['vars_x']]), ' (',
		   			droplevels(phen_1[[i]][['t2_point_ag']]), ')' )),
		   		x = 'antigen', y = 'odd_ratio', fill = 'odd_ratio') +
		   geom_hline (yintercept = 0, linetype = 3) +
		   coord_flip() +
		   facet_wrap( ~ t_igg, ncol = length(levels(phen_1[[i]][['t_igg']])), scales = 'free_x') +
		    geom_text(
        data = phen_1[[i]] %>% filter(raw_pval <= 0.05),
        	aes(label = paste(or)),
        vjust = 0.7, hjust = -0.7, nudge_x = 0.4,
        size = 3, family = 'Georgia', colour = 'black',
        check_overlap = TRUE)

			# list
				l_plots[[length(l_plots) + 1]] = p_1
				}

# return
	return(l_plots)
}




### geom_point plots (for BETA)-------------------------------------------------------------------------------

# '@export
chmi.stat_plot.pointrange_beta <- function(
	phen,
	l_varsy, l_varsx,
	arg_igg = c('IgG', 'IgG1', 'IgG2', 'IgG3', 'IgG4'),
  stat = c('ggpubr'))
{
# stat
  stat <- match.arg(stat, c('ggpubr'))

# variables to avoid
	l_vars <- c('null_deviance', 'df_null', 'log_lik', 'aic', 'bic', 'deviance', 'df_residual',
		'r_squared', 'adj_r_squared', 'sigma', 'mod_statistic', 'mod_pval', 'df', 'deviance', 'df_residual',
		'std_error', 't_value')

	# fix round variables
	phen[, c('beta', 'ci_lower', 'ci_upper')] <- round(phen[, c('beta', 'ci_lower', 'ci_upper')], 2)
	phen[, c('raw_pval', 'adjust_pval')] <- round(phen[, c('raw_pval', 'adjust_pval')], 5)


# argument to filter `l_varsy` | `l_varsx`
if(!missing(l_varsy)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(vars_y), all_vars(. %in% l_varsy)) %>%
    droplevels()
}

if(!missing(l_varsx)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(vars_x), all_vars(. %in% l_varsx)) %>%
    droplevels()
}

if(!missing(arg_igg)) {
# filter 'phen_1'
	phen <- phen %>%
		filter_at(vars(isotype), all_vars(. %in% arg_igg)) %>%
    droplevels()
}


phen <- phen %>%
	mutate(
		beta_signif_raw = factor(ifelse(raw_pval <= 0.05, beta, '')),
    beta_signif_adj = factor(ifelse(adjust_pval <= 0.05, paste0('*'), '')))

# group_split
phen_1 <- phen %>%
  select(-analyte_block) %>%
  select(-one_of(l_vars)) %>%
  filter(complete.cases(.)) %>%
  droplevels() %>%
	group_by(subset, vars_x) %>%
	group_split()

# list
	l_plots <- list()

	### make all plots
		# loop
		for(i in seq_along(phen_1)) {
		# plot 'point_range'
	 		p_1 <- ggplot(
	 				phen_1[[1]],
		    	aes(x = analyte_label, y = beta)) +
		    geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1) +
		    geom_pointrange(aes(ymin = ci_lower, ymax = ci_upper), size = 0.6) +
		    scale_fill_gradientn(colours = c('#56B4E9', 'white', '#8B0000'),
		      breaks = seq(
		      	 min(phen_1[[1]] %>% select(beta)),
		      	 max(phen_1[[1]] %>% select(beta)), by = 0.5),
		      limits = c(
        	 	 min(phen_1[[1]] %>% select(beta)),
		      	 max(phen_1[[1]] %>% select(beta)))) +
		   labs(title = paste0(phen_1[[1]][['subset']]), x = '', y = 'beta', fill = 'beta', size = 5) +
		   geom_hline (yintercept = 0, linetype = 'dotdash', size = 1.5, color = 'coral') +
           #ylim(-0.8, 1) +
		   coord_flip() +
		   facet_wrap( ~ vars_x_label + isotype, ncol = 5, scales = 'free_x') +
           theme(
                plot.title = element_text(color = "black", size = 22, face = "bold", vjust = 6),
                plot.margin = unit(c(4, 4, 4, 1.8), "cm"),
                axis.title.x = element_text(color = "black", size = 16, face = "bold", vjust = -0.7),
                axis.title.y = element_text(color = "black", size = 16, face = "bold"),
                panel.border = element_rect(color = "black", fill = NA, size = 0.5),
                panel.background = element_rect(fill = "white", colour = "white", size = 0.5, linetype = 'solid'),
                panel.grid.major = element_line(size = 0.5, linetype = 'dashed', colour = "grey"),
                panel.grid.minor = element_line(size = 0.25, linetype = 'dashed', colour = "grey"),
                strip.text = element_text(color = "black", size = 14),
                panel.spacing = unit(0.8, "lines"),
                axis.text = element_text(size = 13)) +
		   geom_text(
        data = phen_1[[1]],
        	aes(label = paste0(beta_signif_raw, beta_signif_adj)),
        vjust = 0.3, hjust = 1.5, nudge_x = -0.4, nudge_y = 0.02,
        size = 4, family = 'Arial', colour = 'black',
        check_overlap = TRUE)

			# list
				l_plots[[length(l_plots) + 1]] = p_1
				}

# return
	return(l_plots)
}
mvazquezs/chmitools documentation built on May 1, 2020, 2:06 a.m.