Networks.Rmd

---
title: "Taxonomic network visualization"
author: "Leo Lahti, Sudarshan Shetty et al."
bibliography: 
- bibliography.bib
output:
  BiocStyle::html_document:
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float: true
    self_contained: true
    thumbnails: true
    lightbox: true
    gallery: true
    use_bookdown: false
    highlight: haddock
---
<!--
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{microbiome tutorial - networks}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->


Load [example data](Data.html):

```{r networks1, message=FALSE, warning=FALSE}
library(microbiome)
data(dietswap)
pseq <- dietswap

# Keep only the prevalent taxa to speed up examples
pseq <- core(pseq, detection = 5^2, prevalence = 80/100)
pseq <- subset_samples(pseq, nationality == "AFR" & group == "DI" & bmi_group == "lean")
```


### Taxonomic network reconstruction 

See the [phyloseq tutorial](http://joey711.github.io/phyloseq/plot_network-examples) for
additional network visualization tools.

The widely reported compositionality bias in similarity measures can
be fixed with SpiecEasi or SparCC; the implementations are available
via the [SpiecEasi package](https://github.com/zdk123/SpiecEasi). Note
that the execution is slow.

```{r networks4, warning=FALSE, message=FALSE, fig.width=10, fig.height=10}
# Pick the OTU table
library(phyloseq)
otu <- abundances(pseq)
```


```{r spieceasi, warning=FALSE, message=FALSE, fig.width=10, fig.height=10, eval=TRUE}
# SPIEC-EASI network reconstruction
# In practice, use more repetitions
library(devtools)
#install_github("zdk123/SpiecEasi")
library(SpiecEasi) #install_github("zdk123/SpiecEasi")
net <- spiec.easi(t(otu), method='mb', lambda.min.ratio=1e-2, nlambda=5, icov.select.params=list(rep.num=1))

## Create graph object
n <- net$refit
#colnames(n) <- rownames(n) <- rownames(otu)

# Network format
library(network)
#netw <- network(as.matrix(n), directed = FALSE)

# igraph format
library(igraph)
# ig <- graph.adjacency(n, mode='undirected', add.rownames = TRUE)

# Network layout
# coord <- layout.fruchterman.reingold(ig)

## set size of vertex to log2 mean abundance 
# vsize <- log2(rowMeans(otu))

# Visualize the network
# print(plot(ig, layout = coord, vertex.size = vsize, vertex.label = names(vsize)))
```


Investigate degree distribution with the following:

```{r degree, warning=FALSE, message=FALSE, fig.width=10, fig.height=7, eval=FALSE}
#dd <- degree.distribution(ig)
#plot(0:(length(dd)-1), dd, ylim = c(0,.35), type = 'b',  ylab = "Frequency", xlab = "Degree", main = "Degree Distributions")
```


Visualize the network with [ggnet2](https://briatte.github.io/ggnet):

```{r networks5, warning=FALSE, message=FALSE, fig.width=12, fig.height=7, eval = TRUE}
library(GGally)
#library(ggnet)
library(network)
library(sna)
library(ggplot2)
library(intergraph) # ggnet2 works also with igraph with this

phyla <- map_levels(rownames(otu),
           from = "Genus", to = "Phylum",
           tax_table(pseq))

#netw %v% "Phylum" <- phyla
#p <- ggnet2(netw, color = "Phylum", label = TRUE, label.size = 2)
```


```{r networks5_plot, warning=FALSE, message=FALSE, fig.width=12, fig.height=7, eval = TRUE}
print(p)
```