---
title: Gene set analyses of the HUSH knockdown data
author: Aaron Lun
date: 12 December 2016
---

```{r, echo=FALSE, results='hide'}
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
options(width=100)
```

# Introduction

Here we use the `fry` function from _limma_ to perform gene set analyses.
These are self-contained gene set tests, determining if there is _any_ differential expression across genes in each set.
We use the hallmark, C5 gene ontology and C2 curated gene sets from BROAD's MSigDB.

```{r}
fpath <- "~/AaronDocs/Research/DataFiles/miscellaneous"
load(file.path(fpath, "human_H_v5p2.rdata"))
head(names(Hs.H))
load(file.path(fpath, "human_c5_v5p2.rdata"))
head(names(Hs.c5))
load(file.path(fpath, "human_c2_v5p2.rdata"))
head(names(Hs.c2))
```

We also load in the processed `DGEList` object for our genes.

```{r}
library(edgeR)
y <- readRDS("object.rds")$genes
design <- y$design
head(y$samples)
```

We convert the names to Entrez identifiers.

```{r}
library(org.Hs.eg.db)
anno <- select(org.Hs.eg.db, keys=y$genes$ENSEMBL, keytype="ENSEMBL", columns="ENTREZID")
new.names <- anno$ENTREZID[match(y$genes$ENSEMBL, anno$ENSEMBL)]
summary(is.na(new.names))
```

We also create a directory to save the results.

```{r}
dir.out <- "gene_set_results"
dir.create(dir.out, showWarning=FALSE)
```

# Performing the gene set analyses

## For periphilin

Setting up the contrast for periphilin:

```{r}
con <-  makeContrasts((shPHLplustet - shPHLnotet) - 
                      (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r executor}
out.H <- fry(y, design=design, geneid=new.names, index=Hs.H, contrast=con)
head(out.H)
out.c5 <- fry(y, design=design, geneid=new.names, index=Hs.c5, contrast=con)
head(out.c5)
out.c2 <- fry(y, design=design, geneid=new.names, index=Hs.c2, contrast=con)
head(out.c2)
```
          
Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_PHL.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_PHL.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_PHL.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For TASOR:

Setting up the contrast for TASOR:

```{r}
con <- makeContrasts((shTASORplustet - shTASORnotet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_TASOR.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_TASOR.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_TASOR.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For SETDB-1:

Setting up the contrast for SETDB-1:

```{r}
con <- makeContrasts((shSETDBplustet - shSETDBnotet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_SETDB.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_SETDB.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_SETDB.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For MPP8:

Setting up the contrast for MPP8:

```{r}
con <- makeContrasts((shMPP8plustet - shMPP8notet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_MPP8.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_MPP8.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_MPP8.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For MORC2:

Setting up the contrast for MORC2:

```{r}
con <- makeContrasts((shMORCplustet - shMORCnotet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_MORC.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_MORC.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_MORC.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For KAP:

Setting up the contrast for KAP:

```{r}
con <- makeContrasts((shKAPplustet - shKAPnotet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_KAP.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_KAP.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_KAP.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

## For ATF7IP:

Setting up the contrast for ATF7IP:

```{r}
con <- makeContrasts((shATF7IPplustet - shATF7IPnotet) - 
                     (shControlplustet - shControlnotet), levels=design)
```

Testing each of the gene sets.

```{r, ref.label="executor"}
```

Saving them to file.

```{r}
write.table(file=file.path(dir.out, "H_ATF7IP.tsv"), out.H, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c5_ATF7IP.tsv"), out.c5, sep="\t", quote=FALSE, col.names=NA)
write.table(file=file.path(dir.out, "c2_ATF7IP.tsv"), out.c2, sep="\t", quote=FALSE, col.names=NA)
```

# Session information

```{r}
sessionInfo()
```