DE_varroa.Rmd

---
title: "Varroa microbiome"
author: "Sasha Mikheyev"
date: "1/23/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = T)
```

```{r library}
library(tidyverse)
library(sleuth)
```


```{r sleuth}
bacteria <- read_tsv("mondet/gene_calls_nr.bacteria.tsv", col_names = c("contig", "protein_id", "taxname"), col_types = "-c---c-c") %>% filter(grepl("Bacteria", taxname)) %>% mutate(wolbachia = ifelse(grepl("Wolbachia", taxname), "yes", "no")) %>% separate(taxname, sep = ";\\s+", into = c("Bacteria",LETTERS[1:6]), extra = 'drop') %>% mutate(F = sub(";", "", F))
samples <- tibble(id = dir(file.path("mondet/kallisto"))) %>% mutate(path = paste0("mondet/kallisto/", id)) %>%
  left_join(read_tsv("mondet/treatments.tsv", col_types =  cols())) %>% mutate(sample = id, condition = shortName)
```


### Comparing males vs young

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("young", "males")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)


```
```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionyoung')
results_table = sleuth_results(so, 'conditionyoung')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
male_young_ids = as.vector(DE_genes$target_id)
```

###Young vs Phoretic

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("young", "phoretic")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionyoung')
results_table = sleuth_results(so, 'conditionyoung')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```
```{r}
young_phoretic_ids = as.vector(DE_genes$target_id)
```

### Phoretics vs caged

There were nio

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("phoretic", "phoreticCaged")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionphoreticCaged')
results_table = sleuth_results(so, 'conditionphoreticCaged')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
phor_caged_ids = as.vector(DE_genes$target_id)
```
### Phoretic vs arresting

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("phoretic", "arresting")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```
```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionphoretic')
results_table = sleuth_results(so, 'conditionphoretic')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```

Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
phor_arrest_ids = as.vector(DE_genes$target_id)
```
### Аrresting vs. pre-laying

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("arresting", "prelaying")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionprelaying')
results_table = sleuth_results(so, 'conditionprelaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
arrest_prelaying_ids = as.vector(DE_genes$target_id)
```
### Pre-laying vs. Laying

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("prelaying", "laying")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionprelaying')
results_table = sleuth_results(so, 'conditionprelaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
prelaying_laying_ids = as.vector(DE_genes$target_id)
```
### Laying vs. post-laying

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("laying", "postlaying")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionpostlaying')
results_table = sleuth_results(so, 'conditionpostlaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
laying_postlaying_ids = as.vector(DE_genes$target_id)
```

### Post-laying vs. Emerging

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("postlaying", "emerging")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)

sleuth_significant %>% left_join(bacteria, by = c("target_id" = "contig")) %>% dplyr::filter(!is.na(A)) %>% ggplot(aes(A))+geom_histogram(stat = "count") + scale_y_log10() + coord_flip()
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionpostlaying')
results_table = sleuth_results(so, 'conditionpostlaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
postlaying_emerging_ids = as.vector(DE_genes$target_id)
```
### Post-laying vs. non-laying

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("postlaying", "nonlaying")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionpostlaying')
results_table = sleuth_results(so, 'conditionpostlaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
postlaying_nonlaying_ids = as.vector(DE_genes$target_id)
```
### Laying vs. non-laying

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("laying", "nonlaying")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionnonlaying')
results_table = sleuth_results(so, 'conditionnonlaying')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```
```{r}
laying_nonlaying_ids = as.vector(DE_genes$target_id)
```
### Emerging vs. Phoretic

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("emerging", "phoretic")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionphoretic')
results_table = sleuth_results(so, 'conditionphoretic')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```
```{r}
emerging_phoretic_ids = as.vector(DE_genes$target_id)
```
### Emerging vs. Young

```{r}
so <- sleuth_prep(dplyr::filter(samples, shortName %in% c("emerging", "young")), extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
dim(sleuth_significant)
head(sleuth_significant)
```

```{r}
models(so)
```
Number of Differentially expressed genes (TRUE):
```{r}
so = sleuth_wt(so, 'conditionyoung')
results_table = sleuth_results(so, 'conditionyoung')
Number_DE_genes = table(results_table$qval <= 0.05)
Number_DE_genes
```
Number of up-regulated (TRUE) and down-regulated (FALSE) genes:
```{r}
results_table
DE_genes = subset(results_table, qval <= 0.05)
DE_genes
up_down_regulated = table(DE_genes$b > 0)
up_down_regulated

```

```{r}
emerging_young_ids = as.vector(DE_genes$target_id)
```

```{r}
#fwrite(as.list(DEGs_unique), "Microbe_DEG_IDs.txt")
write.csv(DEGs_unique, "Microbe_DEG_IDs.csv",  row.names = F)

```