---  
title: "Marginal differential expression"
author: "Guillermo Ayala" 
date: "`r Sys.Date()`"
format: 
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true
    pdf:
       toc: true
       number-sections: false
       colorlinks: true  
---  


## Abstract

In this vignette we review the basic functions to perform a marginal
differential expression analysis. 

## Packages  
First, we need to load `Biobase` and `tami` with
```{r}
#| label: mde1
#| eval: false
library(Biobase)
library(tami)
```
or 
```{r}
#| label: mde2
pacman::p_load(Biobase,tami)
```

## Data  
The data set will be loaded from `tamidata` using 
```{r}
#| label: mde3
data(gse21942,package="tamidata")
```

The class (kind of data) of this data set is 
```{r}
#| label: mde4
class(gse21942)
```

Additional information about the data set could be obtained with
```{r}
#| label: mde5
help(gse21942,package="tamidata")
```
or just searching  `gse21942` in the web. 

The phenotipic variables contains the information about the samples. 
Their names are
```{r}
#| label: mde6
varLabels(gse21942)
```
The phenotipic variables can be obtained with
```{r}
#| label: mde7
pData(gse21942)
```

The variable giving the condition is `FactorValue..DISEASE.STATE.`.
```{r}
#| label: mde8
pData(gse21942)[,"FactorValue..DISEASE.STATE."]
```
The absolute frequencies are  
```{r}
#| label: mde9
table(pData(gse21942)[,"FactorValue..DISEASE.STATE."])
```

## t-test for a  given gene (row)

We choose the gene in row 698 in the expression matrix. The t-test comparing
the means assuming a common variance can be applied using stats::t.test.

```{r}
#| label: mde10
x0 = pData(gse21942)[,"FactorValue..DISEASE.STATE."]
y0 = exprs(gse21942)[698,]
t.test(y0 ~ x0,var.equal=TRUE)
```

## t-tests for all genes

It can be used `genefilter::rowttests`.

```{r}
#| label: mde11
tt = genefilter::rowttests(gse21942,pData(gse21942)[,"FactorValue..DISEASE.STATE."])
```
The first rows can be seen with
```{r}
#| label: mde12
head(tt)
```
Note that `tt` is a `data.frame`.

```{r}
#| label: mde13
class(tt)
```
The column headed `statistic` is the statistic of the t-test assuming a common
variance. The column `dm` is the difference of the sample means and
the  column headed `p.value` is the p-value of the t-test. These
are the  **raw p-values**. 

We can evaluated how many genes are significant for a given threshold
$\alpha = 0.05$. 

```{r}
#| label: mde14
#| eval: false
tt[,"p.value"] < 0.001
```

The rows (in the expression matrix) of these genes are 
```{r}
#| label: mde15
#| eval: false
which(tt[,"p.value"] < 0.001)
```

The rows for the minimum and maximum p-value can be found with
```{r}
#| label: mde16
which.min(tt[,"p.value"])
which.max(tt[,"p.value"])
```

A summary of the p-values can be produced using
```{r}
#| label: mde17
summary(tt[,"p.value"])
```
A kernel estimator  of the density of these p-values could be obtained with
```{r}
#| label: mde18
pacman::p_load(ggplot2)
ggplot(tt, aes(x=p.value)) +  geom_density(adjust=1/5)
```

We can adjust the p-values using the Benjamini-Hochberg method 
```{r}
#| label: mde19
p_BH = p.adjust(tt[,"p.value"])
```
or using the Bonferroni correction.
```{r}
#| label: mde20
p_Bonferroni = p.adjust(tt[,"p.value"])
```

Density estimators for Benjamini-Hochberg and Bonferroni adjusted p-values
can be obtained with
```{r}
#| label: mde21
df = data.frame(p_BH)
ggplot(df, aes(x=p_BH)) +  geom_density(adjust=1/5)
```
and
```{r}
#| label: mde22
df = data.frame(p_Bonferroni)
ggplot(df, aes(x=p_Bonferroni)) +  geom_density(adjust=1/5)
```

## q-values

The q-values can be calculated using the package `qvalue`.

```{r}
#| label: mde23
pacman::p_load(qvalue)
qval = qvalue(tt[,"p.value"])
```
Some plots for these q-values just calculated are obtained with

```{r}
#| label: mde24
plot(qval)
```
In particular the bottom-left plot shows the q-values against the raw p-values. 


## t-test using tami::dema   

The function `tami::dema` performs a differential expression analysis
in one step. 

The first option is the usual t-test with a Benjamini-Hochberg 
correction. 

```{r}
#| label: mde25
x1 = tami::dema(x=gse21942,y="FactorValue..DISEASE.STATE.",test = rowt,
          correction = "BH",fdr=.0001,foutput = "gse21942_rowt")
```
A tidy report can be obtained using 
```{r}
#| label: mde26
x1_df = tidy(x1)
```
We have obtained a `data.frame`. The first rows can be obtained with
```{r}
#| label: mde27
head(x1_df)
```
In particular, the results for the gene in row 698 previously considered
are obtained with
```{r}
#| label: mde28
x1_df[698,]
```

The same report in a html file is obtained with

```{r}
#| label: mde29
glimpse(x1)
```

The html file just generated can be open with
```{r}
#| label: mde30
browseURL(glimpse(x1))
```


## Using moderated t-tests  

We are use the moderated t-tests using the package `limma`.
The test used in tami::dema have to be replaced by `rowtmod`.
```{r}
#| label: mde31
x2 = dema(x=gse21942,y="FactorValue..DISEASE.STATE.",
                test = rowtmod,correction = "BH",
                fdr = .0001,foutput = "gse21942_rowtmod")
```
A `data.frame` with the results is obtained with 
```{r}
#| label: mde32
x2_df= tidy(x2)
```
The same report in a html file is obtained with

```{r}
#| label: mde33
glimpse(x2)
```
We can open the file with
```{r}
#| label: mde34
browseURL(glimpse(x2))
```


## Comparing procedures

Let us assume that we have tested the differential expression using two
procedures `rowt` and `rowtmod`. We want to compare the results obtained.

Note that the `data.frame`s  generated using `tidy` provides information
about those genes with a adjusted p-value lesser than the used `fdr`.
If we want all genes in the study we have to set `fdr` to 1. 

How many genes are called significant for a `fdr` equal to .0001 using 
`rowt`?

```{r}
#| label: mde35
nrow(x1_df)
```
And using `rowtmod`?
```{r}
#| label: mde36
nrow(x2_df)
```
We don't want to set previously given `fdr`'s. The `fdr` is set to 1 and
the function `dema` is executed again for both procedures. Note that the values
of each functions is saved in the same objects `x1` and `x2`.

```{r}
#| label: mde37
x1 = dema(x=gse21942,y="FactorValue..DISEASE.STATE.",
                test = rowt,correction = "BH",
                fdr = 1,foutput = "gse21942_rowt")

x2 = dema(x=gse21942,y="FactorValue..DISEASE.STATE.",
                test = rowtmod,correction = "BH",
                fdr = 1,foutput = "gse21942_rowtmod")
```


The new `data.frame`s are generated. 
```{r}
#| label: mde38
x1_df= tidy(x1)
x2_df= tidy(x2)
```
What information contains these `data.frame`s just generated?
```{r}
#| label: mde39
names(x1_df)
```
They are different gene identifiers and the usual `statistic` of the test used
(`rowt` or `rowtmod`), the raw p-value, the corresponding adjusted p-value
(using the `correction` method for multiple testing) and the Storey's  q-value.

Now we set a value for `fdr` 
```{r}
#| label: mde40
fdr0 = .000001
```
and the number of significant for each method would be
```{r}
#| label: mde41
table(x1_df[,"adjp"] < fdr0)
```
and 
```{r}
#| label: mde42
table(x2_df[,"adjp"] < fdr0)
```
We can compare the results for both methods using 
```{r}
#| label: mde43
res1 = x1_df[,"adjp"] < fdr0
res2 = x2_df[,"adjp"] < fdr0
table(res1,res2)
```


