---
title: "Human Microbiome Cooccurance"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Human Microbiome Cooccurance}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

Load packages
```{#r setup, include=FALSE}
library(CoreMicro)
library(phyloseq)
library(parallelDist)
library(vegan)
library(backbone)
library(tidyverse)
library(igraph)
library(visNetwork)
```

```{#r}
hmbp_md<-read.csv("/Users/gordoncuster/Desktop/Manuscript Submissions/In Review/Core_community/November2020/Data_for_additional_analyses/v35_map_uniquebyPSN.txt", sep = "\t")
rownames(hmbp_md)<-hmbp_md$X.SampleID
hmbp_md$X.SampleID = NULL

hmbp_otu<-read.csv("/Users/gordoncuster/Desktop/Manuscript Submissions/In Review/Core_community/November2020/Data_for_additional_analyses/otu_table_psn_v35.txt", sep = "\t")
#end long so it grabs everything
names(hmbp_otu)<-substr(names(hmbp_otu), start = 2, stop = 20)
rownames(hmbp_otu)<-hmbp_otu$.OTU.ID
hmbp_otu$.OTU_ID = NULL

hmpb_tax<-as.character(hmbp_otu$onsensus.Lineage)
#export and split to meet structure for inclusinon in phyloseq object
```

```{#r}
commonalities<-intersect(names(hmbp_otu), rownames(hmbp_md))
hmbp_otu_sub<-hmbp_otu[,names(hmbp_otu) %in% commonalities]
hmbp_md_sub<-hmbp_md[rownames(hmbp_md) %in% commonalities,]

human_tax<-read.csv("/Users/gordoncuster/Desktop/Manuscript Submissions/In Review/Core_community/November2020/Data_for_additional_analyses/Human_taxonomy.csv", header = T)
rownames(human_tax)<-human_tax$OTU_ID
human_tax$OTU_ID=NULL
human_tax<-as.matrix(human_tax)

human_tab<-otu_table(hmbp_otu_sub, taxa_are_rows = T)
human_sample_data<-sample_data(hmbp_md_sub)
human_taxonomy<-tax_table(human_tax)

human_ps<-phyloseq(human_tab, human_sample_data, human_taxonomy)
#subset to those samples that our core analysis was conducted with
hmbp_working_set<-subset_samples(physeq = human_ps, HMPbodysubsite=="Stool")
human_working_set <- prune_taxa(taxa_sums(hmbp_working_set) >=1, hmbp_working_set)
```

Assign core membership

```{#r}
human_core_df <- core_methods(human_stool) %>% data.frame()
table(human_core_df$name, human_core_df$value)

# Create table all taxa that are never assigned.
wide_human <- pivot_wider(human_core_df)
#table(rowSums(wide_human[,5:8]))
wide_human$num_methods <- (rowSums(wide_human[,5:8]))
#tables number assigned by each method to add in paper. 
#table(wide_human$`Proportion of Sequence Reads`)
#table(wide_human$`Proportion of Sequence Reads and Replicates`)
#table(wide_human$`Hard Cut Off`)
#table(wide_human$`Proportion of Sequence Replicates`)

#pull out those common taxa found by all
core_4 <-(wide_human[wide_human$num_methods==4,])$X
#pull out those taxa included by each method
prop_seq_reads_core <- (wide_human[wide_human$`Proportion of Sequence Reads`==1,])$X
prop_seq_readsnrep_core <-(wide_human[wide_human$`Proportion of Sequence Reads and Replicates`==1,])$X
HC_core <- (wide_human[wide_human$`Hard Cut Off`==1,])$X
prop_rep_core <- (wide_human[wide_human$`Proportion of Sequence Replicates`==1,])$X
```

# Adonis

## Prep Data

```{#r}
md_a<-data.frame(sample_data(human_working_set))
md_a$visitno<-as.factor(md_a$visitno)
otu_a<-data.frame(otu_table(human_working_set))
otu_a_t<-data.frame(t(otu_a))
otu_a_t<-data.frame(otu_a_t)
```

## Create Binary and Bray Functions 

```{#r}
adonis_binary <- function(data) {
  dist <- parallelDist(data.matrix(data), method = "binary")
  adonis(dist ~ visitno + sex + RUNCENTER, data=md_a)
}

adonis_bray <- function(data) {
  dist <- parallelDist(data.matrix(data), method = "bray")
  adonis(dist ~ visitno + sex + RUNCENTER, data=md_a)
}
```

### full human microbiome stool dataset

```{#r}
adonis_binary(otu_a_t)
adonis_bray(otu_a_t)
```

### Core 4

Core 4
```{#r}
human_core_4 <- otu_a_t[,names(otu_a_t) %in% core_4]

adonis_binary(human_core_4)
adonis_bray(human_core_4)
```

Prop Seq reads
```{#r}
human_prop_seq_reads_core <- otu_a_t[,names(otu_a_t) %in% prop_seq_reads_core]

adonis_binary(human_prop_seq_reads_core)
adonis_bray(human_prop_seq_reads_core)
```

Prop Read N Rep
```{#r}
human_prop_rnr_core<- otu_a_t[,names(otu_a_t) %in% prop_seq_readsnrep_core]

adonis_binary(human_prop_seq_reads_core)
adonis_bray(human_prop_seq_reads_core)
```

HC
```{#r}
human_HC_core <- otu_a_t[,names(otu_a_t) %in% HC_core]

adonis_binary(human_HC_core)
adonis_bray(human_HC_core)
```

Prop Rep
```{#r}
human_prop_rep_core <- otu_a_t[,names(otu_a_t) %in% prop_rep_core]

adonis_binary(human_prop_rep_core)
adonis_bray(human_prop_rep_core)
```

#Network analysis 
```{#r}
md_a<-data.frame(sample_data(human_working_set))
md_a$visitno<-as.factor(md_a$visitno)
otu_a<-data.frame(otu_table(human_working_set))

#Network of full dataset
#http://pablobarbera.com/big-data-upf/html/02b-networks-descriptive-analysis.html
network_test_dat<-otu_a

network_test_dat2 <- network_test_dat %>% mutate_if(is.numeric, ~1 * (. > 0))

rownames(network_test_dat2)<-rownames(network_test_dat)

network_human<-hyperg(as.matrix(network_test_dat2))

hyperg_human_network_sig <- backbone.extract(network_human, alpha = .0001, class = "igraph", fwer = "bonferroni", narrative = TRUE)

#Gives same output as below.
#sorted_degree_human_nodes<-sort(degree(hyperg_human_network_sig), decreasing = T)

node_extract_full_human<-hyperg_human_network_sig[[]]

nodes_with_sig_full_human<-node_extract_full_human[lapply(node_extract_full_human,length)>0]

full_human_node_w_num_edges<-sort(lengths(nodes_with_sig_full_human), decreasing = T) 
str(full_human_node_w_num_edges)
human_node_list<-names(full_human_node_w_num_edges)
```

Comparison of human nodes with taxa identified by core methods
Intersetion of core assignment methods with global signficance network nodes

```{#r}
uncommon_core_intersect <- function(data) {
  full_human_node_w_num_edges[!names(full_human_node_w_num_edges) %in% intersect(data, human_node_list)]
}

common_core_intersect <- function(data) {
  full_human_node_w_num_edges[names(full_human_node_w_num_edges) %in% intersect(data, human_node_list)]
}
```

```{#r}
#Core_4
uncommon_core_intersect(core_4)
common_core_intersect(core_4)

#Prop_Reads
uncommon_core_intersect(prop_seq_reads_core)
common_core_intersect(prop_seq_reads_core)

#Prop_Reps
uncommon_core_intersect(prop_rep_core)
common_core_intersect(prop_rep_core)

#PropRnR
uncommon_core_intersect(prop_seq_readsnrep_core)
common_core_intersect(prop_seq_readsnrep_core)

#HCs
uncommon_core_intersect(HC_core)
common_core_intersect(HC_core)
```

Barplot of inclusion by group status
```{#r}
taxa_in_core_by_any_method <- unique(
  c(
    as.character(prop_seq_reads_core), 
    as.character(prop_rep_core), 
    as.character(prop_seq_readsnrep_core), 
    as.character(HC_core)
    )
  )

#taxa in network with edge > 1
full_human_node_w_num_edges

#could be accomplished with just the intersect command too. Its currenlty the long way.
shared_core_network_taxa <- common_core_intersect(taxa_in_core_by_any_method)

taxa_only_in_core <- taxa_in_core_by_any_method[!taxa_in_core_by_any_method %in% names(shared_core_network_taxa)]

taxa_only_in_network <- uncommon_core_intersect(taxa_in_core_by_any_method)

# summary of degree of each node. 
# no info is available for taxa only in the core 
# becasue they werent included in the network. 
summary(shared_core_network_taxa)
summary(taxa_only_in_network)

# I think this would be better as a venn diagram actaully. 
counts<-c(
  length(shared_core_network_taxa), 
  length(taxa_only_in_core), 
  length(taxa_only_in_network)
 )

barplot(counts)
```

```{#r, warning=FALSE, message=FALSE}
list(
  core = taxa_in_core_by_any_method,
  network = human_node_list
) %>%
  ggVennDiagram::ggVennDiagram(label = "both", category.names = c("", ""), size = 2, color = "black") +
  annotate("text", x = -1.3, y = 0.8, label = "CORE", fontface =2) +
  annotate("text", x = 5.3, y = 0.8, label = "NETWORK", fontface =2) +
  scale_fill_gradient(high = "grey", low = "white") +
  guides(fill = FALSE) 
```


