04-consensus-network-analysis.Rmd

---
author: V. Keith Hughitt
date: "`r format(Sys.time(), '%d %B, %Y')`"
params:
  settings: ""
  title: ""
  version: ""
  analysis_dirname: "04-consensus-network-analysis"
title: "`r params$title` (`r params$version`)"
output:
  html_document:
    toc: true
    number_sections: true
  pdf_document:
    toc: true
    number_sections: true
    latex_engine: xelatex
---

Introduction
============

**Overview**

This goal of this analysis is to combine information from multiple alternative
network parameterizations for the same dataset into a single consensus network.
While it's not possible to know which, if any, single network parameterization
is optimal, by aggregating across a large number of networks, we may be able to
arrive at a single more robust network topology or cluster partitioning.

**Consensus co-expression network construction**

To construct a consensus co-expression network, we will load the stored
adjacency matrices for each alternative network constructed, and stack them to
create a single combined adjacency matrix. Edges with support from a large
number of networks will thus have larger edge weights while those with less
support will have lower combined edge weights. Hierarchical clustering and
tree-cutting can then be used to derrive the usual network partitioning from
this combined data.

One approach to generating a robust filtered network is to prune edges from the
network that fall below a certain limit. By applying a sufficiently stringent
cutoff, many of the non-robust edges can be removed. Afterwards, genes which
become completely disconnected can also be removed to reduce the size of the
dataset.


Methods
=======

## Setup

```{r load_libraries, message=FALSE, results='hide'}
library('biomaRt')
library('clusteval')
library('doParallel')
library('DT')
library('dynamicTreeCut')
library('flashClust')
library('foreach')
library('goseq')
library('gplots')
library('gridExtra')
library('hpgltools')
library('infotheo')
library('igraph')
library('knitr')
library('knitcitations')
library('matrixStats')
library('RColorBrewer')
library('reshape2')
library('tools')
library('venneuler')
library('viridis')
library('WGCNA')
library('xlsx')
library('tidyverse')

source('../../2015/00-shared/R/annotations.R')
source('../../2015/00-shared/R/count_tables.R')
source('../../2015/00-shared/R/enrichment_analysis.R')
source('../../2015/00-shared/R/filtering.R')
source('../../2015/00-shared/R/plots.R')
source('../../2015/00-shared/R/util.R')
source('../../2015/00-shared/R/wgcna.R')

options(stringsAsFactors=FALSE)
```

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

```{r child='../../2015/00-shared/Rmd/init/load_counts.Rmd'}
```

```{r child='../../2015/00-shared/Rmd/init/load_pathogen_annotations.Rmd', eval=CONFIG$target == 'pathogen'}
```

```{r child='../../2015/00-shared/Rmd/init/load_host_annotations.Rmd', eval=CONFIG$target == 'host'}
```

```{r load_adj_mat}
infile <- sub('04-consensus-network-analysis', '03-consensus-network-construction',
              file.path(output_datadir, paste0('adjmat_', output_suffix, '.rda')))
load(infile)
```

## Filtered consensus co-expression network

```{r child='child/04-create-filtered-consensus-net.Rmd'}
```

```{r child='child/04-detect-consensus-net-modules.Rmd'}
```

```{r child='child/04-module-expression-plots.Rmd'}
```

```{r create_net_annotation_df}
# create dataframe to use for network node annotations
result <- cbind(gene_info, color=module_colors)
annot <- cbind(gene_info, color=module_colors, cluster=module_labels)
```

## Functional enrichment

```{r child='../../2015/00-shared/Rmd/results/go_enrichment_network.Rmd'}
```

```{r child='../../2015/00-shared/Rmd/results/kegg_enrichment_network.Rmd'}
```

```{r child='child/04-kegg-enrichment-results.Rmd'}
```

```{r child='../../2015/00-shared/Rmd/results/hsapiens_marbach2016_tf_enrichment.Rmd', eval=CONFIG$target == 'host' && CONFIG$host == 'H. sapiens'}
```

```{r child='child/04-tf-regulon-enrichment-results.Rmd', eval=CONFIG$target == 'host' && CONFIG$host == 'H. sapiens'}
```

```{r child='../../2015/00-shared/Rmd/results/cpdb_enrichment_network.Rmd', eval=CONFIG$target == 'host'}
```

```{r child='../../2015/00-shared/Rmd/results/secreted_proteins.Rmd', eval=CONFIG$target == 'pathogen'}
```

```{r child='child/00-shared-summary-of-enrichment-results.Rmd'}
```

```{r child='child/04-top-enriched-modules.Rmd'}
```

### Network Visualization

```{r child='child/04-consensus-network-visualization.Rmd'}
```

### Save results

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

```{r save_rdata_results, message=FALSE}
#save(adjmat, file=file.path(output_datadir, 'adjmat_unfiltered.rda'))
message("Saving RData objects")
#save(filtered_adjmat, file=file.path(output_datadir, paste0('adjmat_filtered_', output_suffix, '.rda')))
save(gene_info, file=file.path(output_datadir, paste0('gene_info_', output_suffix, '.rda')))
save(module_go_enrichment, file=file.path(output_datadir, paste0('module_go_enrichment_', output_suffix, '.rda')))
save(module_kegg_enrichment, file=file.path(output_datadir, paste0('module_kegg_enrichment_', output_suffix, '.rda')))

if (CONFIG$target == 'host') {
    save(module_cpdb_enrichment, file=file.path(output_datadir, paste0('module_cpdb_enrichment_', output_suffix, '.rda')))
    save(module_coreg_enrichment, file=file.path(output_datadir, paste0('module_coreg_enrichment_', output_suffix, '.rda')))
}
```

```{r save_network, message=FALSE}
message("Saving consensus network GraphML")

# add secretion / GPI-anchor status cols (pathogen only)
if (CONFIG$target == 'pathogen') {
    annot$gpi_anchored <- gpi_anchored_gene_status
    annot$secreted <- secreted_gene_status
}

# rescale edge weights to 0 - 1 and save network
filtered_adjmat <- log1p(filtered_adjmat)
filtered_adjmat <- filtered_adjmat / max(filtered_adjmat)

filtered_net_filepath <- file.path(output_datadir,
                                   paste0('edge_weight_filtered_network_', 
                                          output_suffix, '.graphml'))

g <- export_network_to_graphml(filtered_adjmat,
                               filtered_net_filepath,
                               threshold=0, max_edge_ratio=5,
                               nodeAttrDataFrame=annot)
```

## Downloads

**Tables**

```{r results='asis'}
# Tables
files <- list.files(output_tabledir)
for (x in files[grepl(output_suffix, files)]) {
    uri <- file.path('table', sub(output_prefix, '', x))
    cat(sprintf('- [%s](%s)\n', basename(x), uri))
}
```

**RData objects**

```{r results='asis'}
# RData objects
files <- list.files(output_datadir)
for (x in files[grepl(output_suffix, files)]) {
    uri <- file.path('data', sub(output_prefix, '', x))
    cat(sprintf('- [%s](%s)\n', basename(x), uri))
}

# free up memory for pandoc
rm(filtered_adjmat)
suppressMessages(gc())
```

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}
# 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, paste0('consesus_network_',
                                                 MANUSCRIPT_CONFIG$output_prefix,
                                                 '.html'))

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