---
title: "Gene Set Analysis"
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  
Different procedures for gene set analysis are proposed in this vignette. 

```{r}
#| label: gsa1
#| echo: false
doAll = TRUE
doNot = FALSE
dirTamiData = "/home/gag/ownCloud/alltami/tami-data/"
```


# Packages  

```{r}
#| label: gsa2
#| results: false
pacman::p_load(Biobase,SummarizedExperiment,tami,ggplot2)
```

# Data

We load the data set from `tamidata`. We are going to analyze two different  
data sets. The first data set has been obtained using microarrays.
The second is a RNA-seq data set.

## gse21942
First an experiment with microarrays applied to samples of homo sapiens. 
```{r gsa3}
#| label: gsa3
data(gse21942,package="tamidata")

```
The phenotipic variable classifying the samples is `FactorValue..DISEASE.STATE.`.
```{r gsa4}
#| label: gsa4
pData(gse21942)[,"FactorValue..DISEASE.STATE."]
```

## PRJNA297664 
Second, an experiment with RNA-seq data of Saccharomyces cerevisiae. 
```{r}
#| label: gsa5
data(PRJNA297664,package="tamidata")
```
The phenotipic variable is `treatment`.
```{r}
#| label: gsa6
colData(PRJNA297664)[,"treatment"]
```

# Differential expression analysis for `gse21942`

##  t-tests for `gse21942`  

We perform a differential expression analysis using  t-tests for
`gse21942`.  
```{r}
#| label: gsa7
gse21942_rowt = dema(x=gse21942,y="FactorValue..DISEASE.STATE.",
                test = rowt,correction = "BH",
                fdr = .0001,foutput = "gse21942")
```

A `data.frame` with the results is obtained with
```{r}
#| label: gsa8
gse21942_rowt_df = tidy(gse21942_rowt)
```
The same information is saved in a html file using 
```{r}
#| label: gsa9
#| eval: false
glimpse(gse21942_rowt)
```
The same information can be studied using an html report with
```{r}
#| label: gsa10
#| eval: false
browseURL(glimpse(gse21942_rowt))
```

Note that only the genes such that the adjusted p-value is
lesser than the corresponding `fdr` are contained in the `data.frame`
or the html report provided by `tidy` or `glimpse` applied over
the object returned by `dema`. 
How many genes (features) are considered in our study?
```{r}
#| label: gsa11
nrow(gse21942)
```
How many genes are returned by `dema`?

```{r}
#| label: gsa12
nrow(gse21942_rowt_df)
```
Let us see a summary of the adjusted p-value for these **significant** genes.
```{r}
#| label: gsa13
summary(gse21942_rowt_df[,"adjp"])
```
In particular, the  maximum adjusted p-value is lesser than `fdr = .0001`.

## Moderated t-tests for `gse21942`  
The test used is now `rowtmod`. It is just the method proposed 
in tha package **Limma**.
We perform a differential expression analysis using a moderated t-tests for
`gse21942`.
```{r}
#| label: gsa14
gse21942_rowtmod = dema(x=gse21942,y="FactorValue..DISEASE.STATE.",
                test = rowtmod,correction = "BH",
                fdr = .0001,foutput = "gse21942")
```
A `data.frame` with the results is obtained with
```{r}
#| label: gsa15
gse21942_rowtmod_df = tidy(gse21942_rowtmod)
```
The same information can be studied using an html report with
```{r}
#| label: gsa16
#| eval: false
browseURL(glimpse(gse21942_rowtmod))
```

# Marginal differential expression analysis for `PRJNA297664`

## `edgeR` method for `PRJNA297664`

First the edgeR method with a common dispersion is used. 

```{r}
#| label: gsa17
PRJNA297664_erc=dema(x=PRJNA297664,y="treatment",test=edgercommon,
                       correction = "BH",fdr= 0.01,foutput = "output")
```

Second, the edgeR method with a tagwise dispersion is used. 

```{r}
#| label: gsa18
PRJNA29766_ert = dema(x=PRJNA297664,y="treatment",
                      test = edgertagwise,correction = "BH",
                      fdr= 0.01,foutput = "output")
```



# Over-representation analysis  for `gse21942`

## Introduction:  exploring a gene set collection

We pretend to perform an over-representation analysis using an unilateral
(or one-tail) Fisher test. The minimun gene set size should be censored to a 
given minimum value, for instance, 10 genes per set and the maximum 
size will be 100 genes per set.

```{r} 
#| label: gsa19
#| results: false
#| echo: false
pacman::p_load(org.Hs.eg.db,GSEABase)
frame = toTable(org.Hs.egGO)
goframeData = data.frame(frame$go_id, frame$Evidence,frame$gene_id)
goFrame = GOFrame(goframeData, organism = "Homo sapiens")
goAllFrame = GOAllFrame(goFrame)
hsa_go = GeneSetCollection(goAllFrame, setType = GOCollection())
hsa_go = geneIds(hsa_go)
save(hsa_go,file=paste0(dirTamiData,"hsa_go.rda"))
```

First we need a `list` of  gene sets.
Son los grupos Gene Ontology para humanos.

```{r}
#| label: gsa20
load(paste0(dirTamiData,"hsa_go.rda"))
```

We have load the gene set collection `hsa_go`.
It is convenient to know what of data we have?

```{r}
#| label: gsa21
class(hsa_go)
```
We can use generic functions for `list`. In particular, it can be know
the number of gene sets with 
```{r}
#| label: gsa22
length(hsa_go)
```
We can access to each element in the `list` using the position
```{r}
#| label: gsa23
hsa_go[1]
```
or the name
```{r}
#| label: gsa24
hsa_go["GO:0000002"]
```
If we want the elements then
```{r}
#| label: gsa25
hsa_go[[1]]
```
or 
```{r}
#| label: gsa26
hsa_go[["GO:0000002"]]
```
The first elements of the `list`.
```{r}
#| label: gsa27
head(hsa_go)
```

## Using `lapply` and `sapply`

If we want to know the number of genes for each gene set. 
```{r}
#| label: gsa28
#| echo: false
ngs = lapply(hsa_go,length)
```
What is `ngs`?
```{r}
#| label: gsa29
#| results: false
class(ngs)
```
Using `unlist` we have a vector with the lengths of the different
gene sets. 
```{r}
#| label: gsa30
#| echo: false
ngs = unlist(ngs)
```
Probably a better choice to calculate the lengths would be
```{r}
#| label: gsa31
#| results: false
ngs = sapply(hsa_go,length)
```
Now `ngs` is 
```{r}
#| label: gsa32
class(ngs)
```
A summary of the lengths is 
```{r}
#| label: gsa33
summary(ngs)
```
We can filter those gene sets with less than 10 genes per set,
```{r}
#| label: gsa34
hsa_go1 = hsa_go[ngs>=10]
```
or those gene sets with a number of genes greater or equal to 10 and
lesser or equal to 10 with
```{r}
#| label: gsa35
hsa_go1 = hsa_go[ngs>=10 & ngs<=10]
```

## `ora` for gse21942  

An over-representation analysis of the data set `gse21942` will be done
using the gene set list `hsa_go` previously loaded. No `userGeneSet`
is provided and this significant gene set is defined as those genes
with an adjusted p-value lesser that `fdr`. The `fdr` value, test
and correction have been given in `dema`.

```{r}
#| label: gsa36
#| results: false
gse21942_ora = overRepresentation(gse21942_rowt,minsize=10,maxsize=100,
                            correction = "BH",
                            GeneSetList = hsa_go,
                            foutput="gse21942_ora")
```

A `data.frame` with the results is obtained with
```{r}
#| label: gsa37
gse21942_ora_df = tidy(gse21942_ora)
```
The report in a html file is obtained as usual as
```{r}
#| label: gsa38
#| eval: false
glimpse(gse21942_ora)
```
It can be opened simply with
```{r}
#| label: gsa39
#| eval: false
browseURL(glimpse(gse21942_ora))
```

## `ora` for PRJNA297664

We are going to do an over-representation analysis of the
dataset PRJNA297664. The significant gene set is calculated 
as the set of genes with an adjusted p-value lesser or equal
to `fdr`.

First we need a gene set list.

```{r}
#| label: gsa40
#| echo: false
pacman::p_load(GSEABase,org.Sc.sgd.db)
frame = toTable(org.Sc.sgdGO)
goframeData = data.frame(frame$go_id, frame$Evidence, frame$systematic_name)
goFrame = GOFrame(goframeData, organism = "Saccharomyces cerevisiae")
goAllFrame = GOAllFrame(goFrame)
gscSc = GeneSetCollection(goAllFrame, setType = GOCollection())
sce_go = geneIds(gscSc)
save(sce_go,file=paste0(dirTamiData,"sce_go.rda"))
```

```{r}
#| label: gsa41
#| echo: false
load(paste0(dirTamiData,"sce_go.rda"))
```

```{r}
#| label: gsa42
#| results: false
PRJNA297664_ora = overRepresentation(PRJNA297664_erc,
                                     minsize=10,maxsize=100,
                                     correction = "BH",
                                     GeneSetList = sce_go,
									 foutput="gse21942_ora")
```
The report in a html file is obtained as usual as
```{r}
#| label: gsa43
#| eval: false
glimpse(PRJNA297664_ora)
```
It can be opened simply with
```{r}
#| label: gsa44
#| eval: false
browseURL(glimpse(PRJNA297664_ora))
```
A `data.frame` with the results is given by
```{r}
#| label: gsa45
#| echo: false
PRJNA297664_ora_df = tidy(PRJNA297664_ora)
```

We can provide the significant gene set. Let us evaluate this significant
gene set and will give it `overRepresentation`. We use the slot 
`GeneStat` of the `DifferentialExpressionOutput` object.
```{r}
#| label: gsa46
#| results: false
GeneStat(PRJNA297664_erc)
```
The significant gene set will be composed by those genes with an adjusted
p-value lesser or equal to 0.001. The row indices in the expression
matrix, `assay(PRJNA297664)`, is 
```{r}
#| label: gsa47
(sig0 = which(GeneStat(PRJNA297664_erc)[,"adjp"] < .0001))
```
Note that the gene set list is specified using `ENTREZID` identifiers.
The significant gene set has to be specified using this identifier. 
```{r}
#| label: gsa48
(SigGeneSet0 = GeneData(PRJNA297664_erc)[sig0,"ENTREZID"])
```

Note that some genes have no `ENTREZID` identifier. It does not matter.
The analysis will remove these genes. Now we call `overRepresentation`
indicating the significant gene set.
```{r}
#| label: gsa49
#| results: false
PRJNA297664_ora = overRepresentation(PRJNA297664_erc,minsize=10,maxsize=100,
                            correction = "BH",
                            GeneSetList = sce_go,
                            SigGeneSet = SigGeneSet0,
                            foutput="gse21942_ora")
```

# Gene set analysis

## Introduction  
The over-representation analysis is a gene set analysis procedure. 
In this section we use other procedures.

## `gse21942`
It means that we use `rowtmod` for the marginal differential expression analysis.
The null distribution to be tested is the `randomization` distribution,
`GeneNullDistr = "randomization"`, the gene set null distribution
is the `self-contained` hypothesis, `GeneSetNullDistr = "self-contained"`,
the maximum number of simulations is given by `nmax`.
The name of the identifiers used in the gene set collection `gsc`
is `id`. The function has as an argument the gene set collection
`gsc`. The enrichment method is a function `descriptive`.
Finally, the file `foutput` is used to save
the html report. 

```{r}
#| label: gsa50
#| results: false
gse21942_self_mean = GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "self-contained",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive=mean,
         foutput = "gse21942_self_mean")
```

Now we reproduce the previous analysis using as `descriptive`
the function `maxmean`.
```{r}
#| label: gsa51
#| results: false
gse21942_self_maxmean = 
  GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "self-contained",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive= maxmean,
         foutput = "gse21942_self_maxmean")
```
Or the `median`.
```{r}
#| label: gsa52
#| results: false
gse21942_self_median = 
  GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "self-contained",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive= median,
         foutput = "gse21942_self_maxmean")
```

We can test the `competitive` hypothesis using as enrichment method 
the `maxmean`.

```{r}
#| label: gsa53
#| results: false
gse21942_comp_mean = 
  GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "competitive",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive=maxmean,
         foutput = "gse21942_comp_mean")
```


```{r}
#| label: gsa54
#| results: false
gse21942_comp_maxmean = 
  GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "competitive",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive=maxmean,
         foutput = "gse21942_comp_maxmean")
```


```{r}
#| label: gsa55
#| results: false
gse21942_comp_median = 
  GeneSetTest(x = gse21942,y="FactorValue..DISEASE.STATE.",
         test = rowtmod,association="statistic",correction="BH",
         GeneNullDistr = "randomization",
         GeneSetNullDistr = "competitive",
         alternative="two-sided",nmax = 100,
         id = "ENTREZID",gsc=hsa_go,descriptive=median,
         foutput = "gse21942_comp_maxmean")
```


Now the report can be obtained as a `data.frame`
```{r}
#| label: gsa56
#| results: false
gse21942_self_mean_df = tidy(gse21942_self_mean)
gse21942_comp_mean_df = tidy(gse21942_comp_mean)
gse21942_self_maxmean_df = tidy(gse21942_self_maxmean)
gse21942_comp_maxmean_df = tidy(gse21942_comp_maxmean)
gse21942_self_median_df = tidy(gse21942_self_median)
gse21942_comp_median_df = tidy(gse21942_comp_median)

```
or like a html report
```{r}
#| label: gsa57
#| eval: false
glimpse(gse21942_self_mean)
glimpse(gse21942_comp_mean)
glimpse(gse21942_self_maxmean)
glimpse(gse21942_comp_maxmean)
```
that can be opened with
```{r}
#| label: gsa58
#| eval: false
browseURL(glimpse(gse21942_self_mean))
browseURL(glimpse(gse21942_comp_mean))
browseURL(glimpse(gse21942_self_maxmean))
browseURL(glimpse(gse21942_comp_maxmean))
```

## `PRJNA297664`
It means that we use `rowtmod` for the marginal differential
expression analysis.
The null distribution to be tested is the `randomization`
distribution and the enrichment measure is the mean of the
phenotype-expression association that it has been chosen as
the `statistic` of the hypothesis testing procedure. 


```{r}
#| label: gsa59
load(paste0(dirTamiData,"sce_go.rda"))
```

```{r}
#| label: gsa60
#| results: false
PRJNA297664_self =
    GeneSetTest(x = PRJNA297664,y="treatment",
                test = edgercommon,association="statistic",
                correction="BH",
                GeneNullDistr = "randomization",
                GeneSetNullDistr = "self-contained",
                alternative="two-sided",nmax = 100,
                id = "ORF",gsc=sce_go,descriptive=mean,
                foutput = "PRJNA297664_self")
```	


```{r}
#| label: gsa61
#| results: false
PRJNA297664_comp =
    GeneSetTest(x = PRJNA297664,y="treatment",
                test = edgercommon,association="statistic",
                correction="BH",
                GeneNullDistr = "randomization",
                GeneSetNullDistr = "competitive",
                alternative="two-sided",nmax = 100,
                id = "ORF",gsc=sce_go,descriptive=mean,
                foutput = "PRJNA297664_comp")
```


Now the report can be obtained as a `data.frame`
```{r}
#| label: gsa62
#| results: false
PRJNA297664_self_df = tidy(PRJNA297664_self)
PRJNA297664_comp_df = tidy(PRJNA297664_comp)
```
or like a html report
```{r}
#| label: gsa63
#| eval: false
glimpse(PRJNA297664_self)
glimpse(PRJNA297664_comp)
```
that can be opened with
```{r}
#| label: gsa64
#| eval: false
browseURL(glimpse(PRJNA297664_self))
browseURL(glimpse(PRJNA297664_comp))
```

