```{r setup, echo=FALSE}
# Set up the environment
library(knitr) 
opts_chunk$set(cache.path='cache/', fig.path='img/', cache=T, tidy=T, fig.keep='high', echo=T, dpi=100, out.width=700)
options(replace.assign=TRUE, width=220)
options(stringsAsFactors = FALSE)
set.seed(1)
```

```{r}
library(DT)
library(limma)
# source("http://bioconductor.org/biocLite.R") # Import biocLite() function into R environment
# biocLite("limma")
library(limma)
# limmaUsersGuide() # Opens pdf manual for limma
# biocLite("samr")
library(samr) # Load the library
library(affy) # Load affy package 
library(GEOquery)
# Use biomaRt for annotation
# biocLite("biomaRt") 
library(biomaRt)
```

SAM (Significance Analysis of Microarrays)
===

We'll be comparing gene expression changes of lung cells in response to different bacterial strains. Data from [Cobb et.al., 2004](https://www.ncbi.nlm.nih.gov/pubmed/15494517), obtained from [GEO, GDS858](https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS858)

```{r, warning=FALSE}
# Download GDS file, put it in the current directory, and load it:
gds858 <- getGEO('GDS858', destdir="data") 
# gds858<-getGEO(filename="../data/GDS858.soft.gz", destdir="data") # If FTP doesn't work, read in from local file
eset <- GDS2eSet(gds858, do.log2=TRUE) # Convert the data to ESET object
# help(ExpressionSet) # If needed, refresh your memory
```

Let's check which infection statuses do we have:

```{r}
pData(eset)$infection # Let's check the infection status
table(pData(eset)$infection) # How many samples are in each infection status
```

Select indexes of the columns associated with "uninfected" or "FRD440" infection status. Check that the right infection status was selected:

```{r}
selected <- grep("uninfected|FRD440", pData(eset)$infection) # Select indexes
pData(eset)$infection[selected] # Check if the right infection status was selected
```

We selected 4 "uninfected" samples and 4 "FRD440" samples. Let's make a vector of outcome measurements, defining groups of replicates. "Uninfected" samples are labeled as "1" and "FRD440" samples are labeled as "2":

```{r}
y <- c(rep(1,4), rep(2,4)) # Vector of outcome measurements
y # Visualize it
```

From the whole expression matrix we select a subset of columns corresponding to the selected indexes. We then quantile normalize the data:

```{r}
exprs.selected <- exprs(eset)[, selected] # Get a subset of expression values of the selected samples
exprs.selected.q <- normalizeQuantiles(exprs.selected) # Quantile normalize the data
```

In order to perform the analysis, SAM requires a list object containing a data object with expression values in a form of _p_ genes by _n_ samples matrix (missing values allowed),  a vector of length n of outcome measurements, vectors of gene names and gene IDs, both of length _p_, and a boolean indicating whether the data is log2-transformed. This object resembles ExpressionSet with slots having different data, but we assemble it as a list.

Now we have our expression matrix (exprs.selected.q), our vector of outcomes (y). Let's use row names of the expression matrix as both gene names and IDs, and assemble our data:

```{r}
genenames <- rownames(exprs(eset)) # Get row names = gene names or IDs
data <- list(x = exprs.selected.q, y = y, geneid = genenames, genenames = genenames, logged2 = TRUE)
```

SAM can perform a variety types of analyses, specified in "resp.type" parameter. Each analysis type requires specific formatting of outcome vector and expression data. Refer to `help(samr)` for details. Now, we're performing "Two class unpaired" analysis:

```{r, results='hide'}
samr.obj <- samr(data, resp.type = "Two class unpaired", nperms=100)
```

Let's check what we have. Everything obvious, isn't it?

```{r}
names(samr.obj) # Get object names from samr.obj
```

Now, we have to choose the delta value that is able to give us the best compromise in terms of called genes, false genes and False Discovery Rate (FDR). In microarray analysis, is very important to have statistically robust results, but we have to keep in mind that too short gene lists are not able to describe the biological meaning of the experiment. In any case, keeping the FDR < 10% (the number of false positives is < 10%) is pretty safe in most cases.

In general, defining the cut-off is a subjective choice and there is no absolute best way to do it.

```{r, results='hide'}
delta.table <- samr.compute.delta.table(samr.obj, min.foldchange = 1.5) # Compute thresholds for different deltas
datatable(delta.table) # Look at the whole range
```

Let's select delta with median FDR <10% - subset the whole delta table and take the first row.

```{r}
delta.table[delta.table[,"median FDR"] < 0.1,][1,] # Check delta corresponding to median FDR ~0.1
delta <- 1.5 # Select the delta
```

Let's select the delta, and have a look at SAM plot:

```{r}
samr.plot(samr.obj, delta) # Check SAM plot
```

Do we have larger number of upregulated genes? Or downregulated? Let's have a look at them:

```{r}
siggenes.table <- samr.compute.siggenes.table(samr.obj, delta, data, delta.table, min.foldchange = 1.5) # Summarize significant genes
names(siggenes.table) # What data we have in the summary list
nrow(siggenes.table$genes.up) # How many upregulated genes
nrow(siggenes.table$genes.lo) # How many downregulated
# Or
siggenes.table$ngenes.up
siggenes.table$ngenes.lo
```

Let's have a look at actual differentially expressed genes:

```{r}
datatable(siggenes.table$genes.up) # Check how table with the results look like
```

We can write them in a clipboard (Windows only), then paste into Excel. Or export them to files:

```{r}
# Write the results in clipboard
# write.table(siggenes.table$genes.up,"clipboard-128",sep='\t')
# Write the results in file
write.table(siggenes.table$genes.up,"results/genes.up.txt",sep='\t') 
write.table(siggenes.table$genes.lo,"results/genes.dn.txt",sep='\t')
```

These results show pretty much everything we need to know. Except the Affy probe IDs do not tell any biological story up front - it would be nice to have at least gene names. Let's extract those IDs and annotate them:

```{r}
# Extract up- and downregulated IDs
up.ids <- siggenes.table$genes.up[, "Gene ID"] # Upregulated
dn.ids <- siggenes.table$genes.lo[, "Gene ID"] # Downregulated
```

We can simply get the platform data from GEO, and extract subsets of IDs:

```{r}
# Annotate Affy probe IDs
Meta(gds858)$platform # Check which platform do we have
# gpl96 <- getGEO('GPL96', destdir = "data") # Get the data for this platform from GEO
gpl96 <- getGEO(filename = "../data/GPL96.soft", destdir = "data")
datatable(Table(gpl96)[Table(gpl96)[,"ID"] %in% up.ids,c("ID","Gene Symbol","Gene Title")]) # Extract annotation for up.ids
```

Or, better, we can use biomaRt. It is a more flexible approach:

```{r}
datatable(listDatasets(useMart("ensembl"))) # List available datasets
mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl") # Load BIOMART dataset for homo sapiens
# Information - lists of filters and attributes
datatable(listFilters(mart)) # Filters, these are our IDs we'll be subsetting the BIOMART annotations on
datatable(listAttributes(mart)) # Attributes, these are annotations we would like to get for our IDs
Meta(gpl96)$title # Check which microarray we have, to select the right attributes
attr <- listAttributes(mart) # Get all attributes as a table
attr[grep("affy", attr[,1]),] # Parse them for anything that looks like from affymetrix
```

Now we know that our IDs correspond to "affy_hg_u133a" attributes in Biomart. Let's extract annotations for them:

```{r}
# Get annotations from Biomart
genes.up <- getBM(attributes = c('affy_hg_u133a', 'external_gene_name', 'description'), filters = 'affy_hg_u133a', values = up.ids, mart = mart)#, uniqueRows = T)
datatable(genes.up) # Check your work
```

Do that for the downregulated genes, and write the results to file.

Clean your workspace

```{r}
# Clean workspace
# unlink(c("GDS858.soft.gz","genes.up.txt","genes.dn.txt","GPL96.soft"))
```

