---
title: "Statistical Bioinformatics as simple as possible"
author: "Guillermo Ayala" 
date: "`r Sys.Date()`"
format: 
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true
    pdf:
       toc: true
       number-sections: false
       colorlinks: true  
---

# Introduction 

## Abstract  
This vignette is concerned with the use of `tami`. 
We propose different ways to do the same things. 
All functionalities are used. This is our purpose. 


# Differential expression analysis for microarrays

## Packages  
```{r}
#| label: ut1
#| echo: false
dirTamiData = "/home/gag/ownCloud/alltami/tami-data/"
doAll = FALSE
```
The following code download and install the package. 
```{r}
#| label: ut2
#| eval: false
url_tami = "http://www.uv.es/ayala/docencia/tami/tami_1.0.tar.gz"
install.packages(url_tami,repos=  NULL,type="source")
```
We load the package. 
```{r}
#| label: ut3
pacman::p_load(tami)
```
We will use data in the package **tamidata**. This package can be installed
with 
```{r}
#| label: ut4
#| eval: false
url_tamidata = "http://www.uv.es/ayala/docencia/tami/tamidata_0.6.tar.gz"
install.packages(url_tamidata,repos=  NULL,type="source")
```
We need some additional packages.
```{r}
#| label: ut5
pacman::p_load("Biobase","hgu133a.db","ggplot2","ggfortify")
```

## How to create an `ExpressionData` starting from scratch

An `ExpressionData`is an S4 class with three slots. 
The first one is 
matrix `exprm` and we use the attribute `rownames` of the 
matrix (`rownames(exprm)`) as annotation.
The second slot is the phenotypic variable, covariable or experimental factor 
with two levels indicating the classification of each column (sample) in one 
experimental group. 

Let us make an `ExpressionData` starting from scratch.
The expression matrix is just a matrix. The expression values (no real interest)
will be random values generated according a standard normal distribution
(null mean and standard deviation one)

```{r}
#| label: ut6
exprm0 = matrix(rnorm(100),ncol=5)
```
Now `exprm0` has no attribute `rownames`.
```{r}
#| label: ut7
rownames(exprm0)
```
We are going to define as `rownames` 
```{r}
#| label: ut8
rownames(exprm0) = letters[1:20]
```
This is the (small) expression matrix.
```{r}
#| label: ut8b
exprm0
```

The second slot is the covariable `groups`, the experimental factor with 
two levels.
The first level will be labelled "A" and the second levels as "B". 
The samples corresponding to group "A" occupies the columns 1, 4 and 5.
The other columns (samples) corresponding to group "B".

```{r}
#| label: ut9
groups0 = factor(c(1,2,2,1,1),levels = 1:2,labels=c("A","B"))
```
Note that the variable or factor `groups0` takes numerical values 1 and 2 but
R knows that it represents categories or groups. The labels of the categories
will be used in plots and summaries. 

Now we can construct the `ExpressionData`. We have two possibilities. 
The simplest and preferred is (we use the `constructor`).
```{r}
#| label: ut10
x  = new("ExpressionData",exprm = exprm0,groups=groups0)
```
A simple summary can be obtained with 
```{r}
#| label: ut11
x
```
The other equivalent procedure to make the `ExpressionData` is
```{r}
#| label: ut12
x = new("ExpressionData",exprm = exprm0, groups = groups0)
```

## From ExpressionSet to ExpressionData  

Usually, we will start our analysis with a data set organized using
a `Biobase::ExpressionSet`. For instance, `tamidata::gse21942`. It is
easy to construct an `ExpressionData`. 
First we load `tamidata::gse21942`.
```{r}
#| label: ut13
data(gse21942,package="tamidata")
x = ExpressionData(exprm = exprs(gse21942), 
                     groups = pData(gse21942)[,"FactorValue..DISEASE.STATE."],
                     type="microarray")
```

## Using `ExpressionData`
We can recover the different slots using the corresponding `accessors`. 
```{r}
#| label: ut14
head(exprm(x)) ## to limit the output
groups(x)
```


## Plots for `ExpressionData`

The first plot shows the mean versus the standard deviation of the gene
expression profile.

```{r}
#| label: ut15
plot(x,1)
```

The second plot is the median versus the interquantile range of the 
expression profile.
```{r}
#| label: ut16
plot(x,2)
```

The third plot displays the two first principal components. The package
`ggfortify` has to be loaded. 

```{r}
#| label: ut17
#| eval: false
library(ggfortify)
plot(x,3)
```

The fourth plot displays the **Tukey mean-difference** plot to compare 
the mean
expressión of the gene per condition.
We consider the mean expression per condition for each gene: $(x_i,y_i)$ where
$i$ is the gene and $x_i$ and $y_i$ and the corresponding means per condition.
The `x` axis is the mean $\frac{x_i+y_i}{2}$ and the `y` axis is the difference
$x_i - y_i$. The left plot corresponds to the original scale and the right plot
uses the natural logarithms of the original data. Perhaps, it is more usual in 
the omics data literature to use the logarithm with base 2 but it is more
natural our choice. 

```{r}
#| label: ut18
#| eval: false
plot(x,4)
```


## Differential Expression Input

The S4 class `DifferentialExpressionInput` provides two additional slots
to control the statistical analysis to be performed. 

We are going to apply a t-test per gene with the p-values  adjusted 
using the Benjamini-Hochberg correction. Also the q-values from the original
(raw) p-values are calculated.

The first option is to apply the method `rowtt` to an **ExpressionData**.

We have to define an `DifferentialExpressionInput` object. 
There are two options.

```{r}
#| label: ut19
x_dei = DifferentialExpressionInput(exprm =exprm(x),
                 groups =  groups(x),association ="pvalue",
                 correction = "BH")

```

```{r}
#| label: ut20
x_dei = new("DifferentialExpressionInput",exprm =exprm(x),
                 groups =  groups(x),association ="pvalue",
                 correction = "BH")

```
Again a short review of the object is returned typing its name.

```{r}
#| label: ut21
x_dei
```

```{r}
#| label: ut22
x_t = rowt(x= exprm(x),y=groups(x))
```

```{r}
#| label: ut23
x_t2 = rowtest(x= exprm(x),y=groups(x),test = rowt,correction = "BH")
```

The value is a `data.frame`.
```{r}
#| label: ut24
head(x_t2)
```
Note that the rows are identified using the original Affymetrix identifiers.
They are the row names of the data frame.
```{r}
#| label: ut25
head(rownames(x_t2))
```

They will be used as **primary** gene identifiers.
Instead of the t-tests where the standard error of the difference of means is
estimated using the profile expression of each gene, other possible procedure
is contained in the Limma  R package where an hierarchical model is used. 
Again the same results can be obtained using two procedures implemented in tami.


```{r}
#| label: ut26
x_tm = rowtmod(x= exprm(x),y=groups(x))
head(x_tm)
```
The same analysis and the adjusted p-values and q-values can be obtained using 
tami::rowtest. 

```{r}
#| label: ut27
x_tm2 = rowtest(x= exprm(x),y=groups(x),test = rowtmod,correction = "BH")
```

The first rows of the resulting data frame can be inspected  with
```{r}
#| label: ut28
head(x_tm2)
```

Similarly we can applied the `rowt` option. 
```{r}
#| label: ut29
x_tm2 = rowtest(x= exprm(x),y=groups(x),test = rowt,correction = "BH")
```

The first rows of the resulting data frame can be inspected  with
```{r}
#| label: ut30
head(x_tm2)
```


## Differential Expression Output

We have now a data frame with the annotation data and a data frame with the 
statistical results of the differential expresion analysis.
We create a **DifferentialExpressionOutput**. We can made it with two 
procedures.


```{r}
#| label: ut31
x_t2_deo = DifferentialExpressionOutput(GeneData=fData(gse21942),
                                        GeneStat = x_t2,
                                        foutput ="gse21942_rowtt",fdr = .2)
```
We can access or modify the different slots using 
```{r}
#| label: ut32
#| output: false
GeneData(x_t2_deo)
GeneStat(x_t2_deo)
foutput(x_t2_deo)
fdr(x_t2_deo)
```


```{r}
#| label: ut33
x_t2_deo = new("DifferentialExpressionOutput",
                        GeneData = fData(gse21942), GeneStat = x_t2,
                        foutput = "gse21942_rowtt",fdr = .2)
```

We generate the corresponding tidy report where only genes with an adjusted
p-value lower than `FDR` are included in the report.
```{r,results='hide'}
#| label: ut34
#| output: false
tami::tidy(x_t2_deo)
```
We can see the report at [this file](`r paste(getwd(),"/reports/",x_t2_deo@foutput,"_tidy.html",sep="")`).

However, the report with all genes can be obtained with
```{r}
#| label: ut35
#| eval: false
tami::augment(x_t2_deo)
```
The corresponding report can be found [here](`r paste(getwd(),"/reports/",x_t2_deo@foutput,"_augment.html",sep="")`).


### Using rowttmod 
We can repeat the previous analysis by replacing tami::rowt with the moderated
t-tests (from limma package). Note that we have to modify the statistical results
contained in the data frame `gse21942_rowt` by the data frame
`gse21942_rowtmod`.

```{r}
#| label: ut36
x_tm2_deo = new("DifferentialExpressionOutput",
                        GeneData = fData(gse21942), GeneStat = x_tm2,
                        foutput = "gse21942_rowttmod",fdr = .2)
```
Generate the report.
```{r}
#| label: ut37
#| output: false
tami::tidy(x_tm2_deo)
```

The results are included in
[this new file](`r paste0(getwd(),"/reports/",foutput(x_tm2_deo),"_tidy.html")`)

```{r}
#| label: ut38
#| eval: false
dema(x=gse21942,y="FactorValue..DISEASE.STATE.",correction = "BH",test = rowt,
     foutput = "gse21942",fdr = .01)
```


## Gene set analysis 

We have to choose a gene set collection. The package `EnrichmentBrowser`
have some useful functions to construct these gene set collections.
Note that we are going to analyze human and yeast data sets.


### KEGG  
First, we can use the gene sets from KEGG. We download the collections and
save them. The species names can be found at [](http://www.genome.jp/kegg/catalog/org_list.html).

```{r }
#| label: ut39
pacman::p_load(EnrichmentBrowser)
```


```{r}
#| label: ut40
#| eval: false
hsa_kegg = get.kegg.genesets("hsa")
```

```{r}
#| label: ut41
#| eval: false
sce_kegg = get.kegg.genesets("sce")
```
### Gene Ontology  

```{r}
#| label: ut42
#| eval: false
hsa_go = get.go.genesets(org="hsa", onto="BP", mode="GO.db")
```

```{r}
#| label: ut43
#| eval: false
sce_go = get.go.genesets(org="sce", onto="BP", mode="GO.db")
```
As we are going to see the genes are coded using the `ENTREZID`.
## Exploring a gene set collection  
```{r,include=FALSE}
#| label: ut44
#| echo: false
load(paste0(dirTamiData,"hsa_go.rda"))
```


What kind of data we have?
```{r}
#| label: ut45
class(hsa_go)
```
We can use generic functions for `list`.
```{r}
#| label: ut46
length(hsa_go)
```
We can access to each element in the `list` using the position
```{r}
#| label: ut47
hsa_go[1]
```
or the name
```{r}
#| label: ut48
hsa_go["GO:0000002_mitochondrial_genome_maintenance"]
```
If we want the elements then
```{r}
#| label: ut49
hsa_go[[1]]
```
or 
```{r}
#| label: ut50
hsa_go[["GO:0000002_mitochondrial_genome_maintenance"]]
```
The first elements of the `list`.
```{r}
#| label: ut51
head(hsa_go)
```

## Using `lapply` and `sapply`

If we want to know the number of genes for each gene set.
```{r}
#| label: ut52
ngs = lapply(hsa_go,length)
```

What is `ngs`?
```{r}
#| label: ut53
class(ngs)
```
Probably a better choice to calculate the lengths would be
```{r}
#| label: ut54
ngs = sapply(hsa_go,length)
```

Now `ngs` is 
```{r}
#| label: ut55
class(ngs)
```

A summary of the lengths is 
```{r}
#| label: ut56
summary(ngs)
```
Perhaps a kernel density estimator could show the gene set length distribution.

```{r}
#| label: ut57
df = data.frame(ngs)
ggplot(data=df,aes(x=ngs)) + geom_density()
```
