DESEQ2_bio.0.Rmd

---
title: "Differential Expression Analysis"
author: "Adam Festa"
date: "5/8/2018"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

## Script objective

The purpose of this script is to conduct differntial expression on the full set of families.  

```{r cars}
suppressPackageStartupMessages(library(DESeq2))
```

## Load data

Load the biological replicate data

```{r pressure, echo=FALSE}
load("/mnt/media/disk6/ARF/DiffExpression/counts/86kSalmon/Step2_Load_Counts.RData")
```


## Remove data with less than 3 biological replicates

```{r filt.dat}
all.fams <- as.character(load.counts$phenos$fam_id)

remove.fams <- names(which(table(all.fams) < 3))
remove.rows <- which(all.fams %in% remove.fams)

load.counts$txi_object$abundance <- load.counts$txi_object$abundance[,-remove.rows]
load.counts$txi_object$counts <- load.counts$txi_object$counts[,-remove.rows]
load.counts$txi_object$length <- load.counts$txi_object$length[,-remove.rows]
load.counts$txi_object$variance <- load.counts$txi_object$variance[,-remove.rows]

load.counts$phenos <- droplevels.data.frame(load.counts$phenos[-remove.rows,])

```

## Create DE object

```{r}
ddsTxi <- DESeqDataSetFromTximport(load.counts$txi_object,
                                   colData = load.counts$phenos,
                                   design = ~ 0 + batch + index_seq + fam_id)
```


## Run DESeq2

```{r}
library(BiocParallel)
register(MulticoreParam(30))
ddsTxi <- DESeq(ddsTxi,parallel = T)
```

Logs:
```
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates: 30 workers
mean-dispersion relationship
final dispersion estimates, fitting model and testing: 30 workers
```

## Save output

```{r }
save(ddsTxi,file = "/mnt/media/disk6/ARF/DiffExpression/ddstxi_bio_0.RData",compress = T)
```


## Generate all pairwise comparison 

```{r}
all.fams <- as.character(unique(ddsTxi$fam_id))
all.de.txpts <- vector("list")
for(each.fam in 1:length(all.fams)) {
  # Specify the family that everything will be compared to 
  first.comp <- all.fams[1]
  # Specify the list of other families that it will be compared against
  other.comp.fams <- all.fams[-1]
  # Identify DE transcripts for each family
  temp.out <- lapply(1:length(other.comp.fams),function(each.other.fam){ 
    other.fam <- other.comp.fams[each.other.fam]
    results(ddsTxi, 
            contrast = c("fam_id",first.comp,other.fam), 
            parallel = T)
  })
  all.de.txpts[[each.fam]] <- temp.out
  names(all.de.txpts)[[each.fam]] <- first.comp
  all.fams <- all.fams[-1]
}
```

### Save Workspace 
 
```{r save.wkspace}
all.de.txpts.0 <- all.de.txpts
save(all.de.txpts.0,file = "all.de.txpts.0.deseq2_bio.Rdata",compress = T)
```