01-individual-network-enrichment.Rmd

---
author: V. Keith Hughitt
date: "`r format(Sys.time(), '%d %B, %Y')`"
params:
  version: "7.0"
  analysis_dirname: "01-individual-network-enrichment"
title: "Functional Enrichment of Individual Co-expression Networks (`r params$version`)"
output:
  html_document:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: hide
  pdf_document:
    toc: true
    latex_engine: xelatex
---

Overview
========

The purpose of this analysis is to compare co-expression networks
and network partitions resulting from alternative parameterizations, and to
gauge the impact of each parameter choice on the resulting network, across
multiple datasets.

The basic approach used here is to construct multiple networks with varying
parameters, and for each network, measure the level of functional enrichment of
GO terms, KEGG and CPDB pathways, and known TF-regulated genes.

In this file, we summarize the findings for such networks. The scripts used
to generate the networks can be found in the `scripts/` directory.

**Parameters Tested**

For each dataset, all possible combinations of the following parameters will be
tested:

- Counts-per-million (CPM) transformation (TRUE|FALSE)
- Log2 transformation (TRUE|FALSE)
- Quantile normalization (TRUE|FALSE)
- Batch adjustment, when available (limma|combat|none)
- Similarity Measure (pearson correlation, spearman correlation, biweight
  mid-correlation, cor-dist)
- Adjacency power (1-14)

This will result in a total of 1344 networks per dataset where batch
information is available, and 448 networks otherwise.

Other parameters previously tested include:

- Network sign (`abs(sim)` vs. `(1 + cor) / 2`)
- Topological overlap transformation (TOM)
- Hybrid tree cut algorithm
- dynamicTreeCut `deepSplit=FALSE` option
- Low-count filtering
- Low-variance filtering
- Module merge correlation
- Voom transformation

In each of the above cases, the transformations either had little impact on the
resulting network, or had a consistently positive or negative effect across
many datasets, and were thus fixed to allow for alternative parameters to
be tested.

Results
=======

## Setup

```{r load_libraries, message=FALSE}
library('DT')
library('plyr')
library('gridExtra')
library('heatmaply')
library('knitr')
library('reshape2')
library('pander')
library('tidyverse')
source('../../2015/00-shared/R/util.R')
```

```{r child='child/00-shared-setup.Rmd'}
```

```{r plot_adjustments, echo=FALSE}
# alternative plot / display dimensions for single plots
fig_height <- 4.5
out_width  <- 480
out_height <- 300
```

## Host networks

<!--
### *H. sapiens* all samples (HsAll)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_hsall.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
## create a copy to create a combined dataset for comparative analysis
#combined <- cbind(df, network='HsAll')
#```
-->

### *H. sapiens* infected with *T. cruzi* (HsTc)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_hstc.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
## create a copy to create a combined dataset for comparative analysis
combined <- cbind(df, network='HsTc')

# add to combined result dataframe
#l <- list(combined, cbind(df, network='HsTc'))
#combined <- do.call(rbind.fill, l)
```

### *H. sapiens* infected with *T. cruzi* - uninfected samples (HsTcUI)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_hstc-uninf.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
# add to combined result dataframe
l <- list(combined, cbind(df, network='HsTcUI'))
combined <- do.call(rbind.fill, l)
```

<!-- 
### *H. sapiens* infected with *L. braziliensis* (HsLb)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_hslb.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
## add to combined result dataframe
#l <- list(combined, cbind(df, network='HsLb'))
#combined <- do.call(rbind.fill, l)
#```
-->

### *H. sapiens* infected with *L. major* (HsLm)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_hslm.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
# add to combined result dataframe
l <- list(combined, cbind(df, network='HsLm'))
combined <- do.call(rbind.fill, l)
```

<!--
### *H. sapiens* infected with *L. major*, including alt parameters (HsLm-full)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/hsapiens_infected_with_lmajor-full-v6.0.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```
-->


### *H. sapiens* infected with *L. major* - uninfected samples (HsLmUI)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_hslm-uninf.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
l <- list(combined, cbind(df, network='HsLmUI'))
combined <- do.call(rbind.fill, l)
```

<!--
### *M. musculus* infected with *L. major* (MmLm)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_mmlm.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```
#```{r}
## add to combined result dataframe
#l <- list(combined, cbind(df, network='MmLm'))
#combined <- do.call(rbind.fill, l)
#```
-->

## Parasite networks

<!--
### *L. major* infecting *M. musculus* (LmMm)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_lmmm.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
#l <- list(combined, cbind(df, network='LmMm'))
#combined <- do.call(rbind.fill, l)
#```
-->

### *L. major* all hosts (LmAll)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_lmall.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
l <- list(combined, cbind(df, network='LmAll'))
combined <- do.call(rbind.fill, l)
```

### *L. major* infecting *H. sapiens* (LmHs)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_lmhs.csv') 
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
# add to combined result dataframe
l <- list(combined, cbind(df, network='LmHs'))
combined <- do.call(rbind.fill, l)
```

### *T. cruzi* infecting *H. sapiens* (TcHs)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_tchs.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
# add to combined result dataframe
l <- list(combined, cbind(df, network='TcHs'))
combined <- do.call(rbind.fill, l)
```

### *T. cruzi* all developmental stages (TcAll)

```{r load_data, message=FALSE}
df <- read_csv('input/indiv_network_stats/networks_tcall.csv')
```

```{r child='child/01-indiv-net-load-network-stats.Rmd'}
```

```{r child='child/01-indiv-net-results.Rmd'}
```

```{r}
# add to combined result dataframe
l <- list(combined, cbind(df, network='TcAll'))
combined <- do.call(rbind.fill, l)
```

<!--
## modENCODE networks

### ModENCODE - Fly (Fly)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_modencode_fly.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
## add to combined result dataframe
#l <- list(combined, cbind(df, network='Fly'))
#combined <- do.call(rbind.fill, l)
#```

### ModENCODE - Worm (Worm)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_modencode_worm.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
## add to combined result dataframe
#l <- list(combined, cbind(df, network='Worm'))
#combined <- do.call(rbind.fill, l)
#```

### Illumina BodyMap network

#### *H .sapiens* BodyMap multi-tissue network (Body)

#```{r load_data, message=FALSE}
#df <- read_csv('input/indiv_network_stats/networks_bodymap.csv')
#```

#```{r child='child/01-indiv-net-load-network-stats.Rmd'}
#```

#```{r child='child/01-indiv-net-results.Rmd'}
#```

#```{r}
## add to combined result dataframe
#l <- list(combined, cbind(df, network='Body'))
#combined <- do.call(rbind.fill, l)
#```
-->

## Comparison

### Log2 transformation

```{r summary_log2, fig.height=3.375, out.width=out_width, out.height=out_height}
ggplot(combined, aes(network, enrichment_score)) + 
    geom_boxplot(aes(fill=factor(log2_transform)), size=0.25, outlier.size=0.25) +
    theme_bw_high_res() +
    theme(legend.title = element_blank(), 
          axis.text.x=element_text(angle=90, hjust=1)) +
    ylab('Enrichment score')
```

### CPM transformation

```{r summary_cpm, fig.height=3.375, out.width=out_width, out.height=out_height}
ggplot(combined, aes(network, enrichment_score)) + 
    geom_boxplot(aes(fill=factor(cpm_transform)), size=0.25, outlier.size=0.25) +
    theme_bw_high_res() +
    theme(legend.title = element_blank(), 
          axis.text.x=element_text(angle=90, hjust=1)) +
    ylab('Enrichment score')
```

### Quantile Normalization

```{r summary_qnorm, fig.height=3.375, out.width=out_width, out.height=out_height}
ggplot(combined, aes(network, enrichment_score)) + 
    geom_boxplot(aes(fill=factor(quantile_normalize)), size=0.25, outlier.size=0.25) +
    theme_bw_high_res() +
    theme(legend.title = element_blank(), 
          axis.text.x=element_text(angle=90, hjust=1)) +
    ylab('Enrichment score')
```

### Similarity Measure

```{r summary_sim_meas, fig.height=3.375, out.width=out_width, out.height=out_height}
ggplot(combined, aes(network, enrichment_score)) + 
    geom_boxplot(aes(fill=factor(similarity_measure)), size=0.25, outlier.size=0.25) +
    theme_bw_high_res() +
    theme(legend.title = element_blank(), 
          axis.text.x=element_text(angle=90, hjust=1)) +
    ylab('Enrichment score')
```

### Batch Adjustment

```{r summary_batch_adjust, fig.height=3.375, out.width=out_width, out.height=out_height}
# exclude networks without batch information
missing_batch <- names(table(combined$network)[table(combined$network) < 1344])

batch_dat <- combined %>%
    filter(!network %in% missing_batch)

ggplot(batch_dat, aes(network, enrichment_score)) + 
    geom_boxplot(aes(fill=factor(batch_adjustment)), size=0.25, outlier.size=0.25) +
    theme_bw_high_res() +
    theme(legend.title = element_blank(), 
          axis.text.x=element_text(angle=90, hjust=1)) +
    ylab('Enrichment score')
```

### Summary

```{r dataset_similarity, message=FALSE, dpi=192}
combined_wide <- combined %>% 
    arrange(network_id, network) %>%
    select(network_id, network, enrichment_score) %>%
    spread(key=network, value=enrichment_score) %>%
    select(-network_id)

# compare datasets with and without batch separately
batch_mask <- sapply(combined_wide, function(x) { sum(is.na(x)) }) == 0

# datasets with batch
cor_mat_batch <- cor(combined_wide[, batch_mask, drop = FALSE], method='spearman')
diag(cor_mat_batch) <- 0
heatmaply(cor_mat_batch)

# without batch
#cor_mat_no_batch <- cor(combined_wide[complete.cases(combined_wide), !batch_mask, drop = FALSE], method='spearman')
#diag(cor_mat_no_batch) <- 0
#heatmaply(cor_mat_no_batch)
```

```{r summary_table}
combined %>% 
    group_by(network) %>% 
    summarise(
        max_go_uniq=max(unique_go_terms),
        max_kegg_uniq=max(unique_kegg_terms),
        max_cpdb_uniq=max(unique_cpdb_terms),
        max_marbach_uniq=max(unique_marbach_terms),
        max_go_pval=max(total_go_pval),
        max_kegg_pval=max(total_kegg_pval),
        max_cpdb_pval=max(total_cpdb_pval),
        max_marbach_pval=max(total_marbach_pval)
    ) %>% arrange(desc(max_go_uniq)) %>%
    kable
```

```{r summary_overview}
summary_dat <- cbind(
    combined %>% 
        group_by(network, log2_transform) %>% 
        summarize(mean_score=mean(enrichment_score)) %>%
        spread(log2_transform, mean_score) %>%
        summarise(log2=`TRUE` - `FALSE`),
    cpm=combined %>% 
        group_by(network, cpm_transform) %>% 
        summarize(mean_score=mean(enrichment_score)) %>%
        spread(cpm_transform, mean_score) %>%
        summarise(cpm=`TRUE` - `FALSE`) %>%
        pull(cpm),
    qnorm=combined %>% 
        group_by(network, quantile_normalize) %>% 
        summarize(mean_score=mean(enrichment_score)) %>%
        spread(quantile_normalize, mean_score) %>%
        summarise(qnorm=`TRUE` - `FALSE`) %>%
        pull(qnorm)
)

rownames(summary_dat) <- summary_dat$network
summary_dat <- summary_dat %>%
    select(-network) %>%
    t %>%
    melt %>%
    setNames(c('Parameter', 'Network', 'Value'))

ggplot(summary_dat, aes(Network, Parameter)) +
    geom_tile(aes(fill=Value)) + 
    geom_text(aes(label=round(Value, 2)), size=3) +
    scale_fill_gradient(low="green", high="red") + 
    theme_bw_high_res() +
    ylab("Parameter")
```

System Information
==================

```{r sysinfo}
sessionInfo()
```

```{r git_commit}
# manuscript-specific code (https://github.com/elsayed-lab/manuscript-coex-networks)
system('git rev-parse --short HEAD', intern = TRUE)

# shared code (https://github.com/elsayed-lab/manuscript-shared-rnaseq)
system('git --git-dir=$UMD/2015/00-shared/.git rev-parse --short HEAD', intern = TRUE)
```

```{r save_output, include=FALSE}
# Give rmarkdown some time to finish convern and then copy original rmarkdown
# along with output images and HTML to archive location
html_filepath <- file.path(output_prefix, 'individual-network-enrichment.html')
pdf_filepath <- sub('html', 'pdf', html_filepath)

if (opts_knit$get("rmarkdown.pandoc.to") == 'latex') {
    system(sprintf('(sleep 30 && mv 01-individual-network-enrichment.pdf %s) &', pdf_filepath))
} else {
    system(sprintf('(sleep 30 && mv 01-individual-network-enrichment.html %s) &', html_filepath))
}
```