---
title: "gse28619"
author: " Maria Teresa Rubio Martinez-Abarca and Guillermo Ayala" 
date: "`r Sys.Date()`"
format: 
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true
    pdf:
       toc: true
       number-sections: false
       colorlinks: true  
---  

## Introducción

The dataset has been downloaded from
[GEO](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28619)

The summary given there is
    > Alcoholic hepatitis (AH) is the most severe form of alcoholic liver 
    disease and occurs in patients with excessive alcohol intake It is
    characterized by marked hepatocellular damage, steatosis and pericellular
    fibrosis. Patients with severe AH have a poor short-term prognosis.
    Unfortunately, current therapies (i.e. corticosteroids and pentoxyphylline) 
    are not effective in many patients and novel targeted therapies are 
    urgently needed. The development of such therapies is hampered by a 
    poor knowledge of the underlying molecular mechanisms. Based on studies 
    from animal models, TNF alfa was proposed to play a pivotal role in the
    mechanisms of AH. Consequently, drugs interfering TNF alfa were tested 
    in these patients. The results were disappointing due to an increased
    incidence of severe infections. Unluckily, there are not experimental 
    models that mimic the main findings of AH in humans. To overcome this
    limitation, translational studies with human samples are required. We
    previously analyzed samples from patients with biopsy-proven AH. In these
    previous studies, we identified CXC chemokines as a potential therapeutic
    target for these patients. We expanded these previous observations by
    performing a high-throughout transcriptome analysis.

## Downloading and preprocessing the dataset 
The vignette `tamidata::gse28619` details how to download and preprocessing
the dataset. The final result obtained using the code given there
is an `ExpressionSet` used from now on. 

## Marginal differential expression 

- Loading basic packages needed later. 

```{r}
#| label: gse2
library(Biobase)
library(tami)
```

- Loading the dataset  tamidata::gse28619.  

```{r}
#| label: gse3
data(gse28619,package="tamidata")
```

- The number of features (genes) and samples are 
```{r}
#| label: gse4
dim(gse28619)
```


- Phenotipic variable?
```{r}
#| label: gse5
data(gse28619,package="tamidata")
```

- The number of features (genes) and samples are 
```{r}
#| label: gse6
dim(gse28619)
```


- Phenotipic variable?
```{r}
#| label: gse7
pData(gse28619)[,"type"]
```
- Gene identifiers (not shown)
```{r}
#| label: gse8
#| eval: false
fData(gse28619)
```

## Marginal differential expression  

The test used is `rowt` corresponding to the t-test assuming a common variance.

```{r}
#| label: gse9
x1 = dema(x=gse28619, y = "type", test = rowt, correction = "BH",
fdr = 0.01,foutput = "gse28619")
```
Note that only those genes with an adjusted p-value lesser than
`fdr` are returned. If we want the analysis for all genes then we have
to set `fdr = 1`.

A `data.frame` with the results can be obtained. 
```{r}
#| label: gse10
x1_df = tidy(x1)
```
These results can be examined in a html file. 
```{r}
#| label: gse11
browseURL(glimpse(x1))
```
How many genes have an adjusted p-value lesser than `fdr = 0.01`?
```{r}
#| label: gse12
nrow(x1_df)
```

Which is the highest adjusted p-value for the returned genes?
```{r}
#| label: gse13
max(x1_df[,"adjp"])
```
We want to know all the information about the gene with `ENTREZID` identifier
`6898`. First we need to know the row within the `data.frame` `x1_df`.
```{r}
#| label: gse14
grep("^6898",x1_df[,"ENTREZID"])
```
The expression "^6898" looks for `ENTREZID` beginning with this code.
Let us see these rows
```{r}
#| label: gse15
x1_df[grep("^6898",x1_df[,"ENTREZID"]),]
```
Why we have two rows? Note that the `PROBEID`'s are different for each row. 
They are different probe sets but they corresponds to the same gene. 

Note that the `ENSEMBL` and `GO` identifiers are `ENSG00000198650`
and `GO:0004838` respectively. 



## Gene set collection  

```{r,include=FALSE}
#| label: gse16
dirTamiData = "~/ownCloud/alltami/tami-data/"
load(paste0(dirTamiData,"hsa_go.rda"))
```

```{r,eval=FALSE}
#| label: gse17
load("hsa_go.rda")
hsa_go
```

- How many gene sets?
```{r}
#| label: gse18
length(hsa_go)
```
The first gene set is 
```{r}
#| label: gse19
hsa_go[1]
```
We have the name and the (ENTREZID) identifiers of the genes belonging to
this gene set.
The elements of the set can be accessed with
```{r}
#| label: gse20
hsa_go[[1]]
```
We can know the number of elements in this groups with
```{r}
#| label: gse21
length(hsa_go[[1]])
```
It is not necessary to evaluate the `length` (cardinal) of each group. 
It can be known the cardinality of each set simultaneously using 
`lapply`.
```{r}
#| label: gse22
ngs = lapply(hsa_go,length)
```
What kind of object has been returned?
```{r}
#| label: gse23
class(ngs)
```
The (previously calculated) length of the first gene set is 
```{r}
#| label: gse24
ngs[1]
```
or
```{r}
#| label: gse25
ngs[[1]]
```
The lengths of the first gene sets can be obtained with
```{r}
#| label: gse26
ngs[1:10]
```
The largest gene set is 
```{r}
#| label: gse27
which.max(ngs)
```
How many gene sets have exactly 10 genes?
```{r}
#| label: gse28
table(ngs == 10)
```
Other possibility is 
```{r}
#| label: gse29
sum(ngs == 10)
```
Note that the lengths calculated `ngs` is a `list`. Sometimes it 
is useful to have a vector. 
```{r}
#| label: gse30
ngs0 = unlist(ngs)
class(ngs0)
```
We can evaluate a table of absolute frequencies with 
```{r}
#| label: gse31
table(ngs0)
```
How many gene sets have more than 5 genes?
```{r}
#| label: gse32
sum(ngs0 > 5)
```
Or greater or equal to 5? 
```{r}
#| label: gse33
sum(ngs0 >= 5)
```
How many groups have a number of genes greater than 7 and less than 90?
```{r}
#| label: gse34
sum(ngs0 > 7 & ngs0 < 90)
```
How many groups have less than 34 or more than 67 genes?
```{r}
#| label: gse35
sum(ngs0 < 34 | ngs0 > 67)
```

## Over representation analysis 

We have loaded the gene set collection `hsa_go` previously downloaded.
```{r}
#| label: gse36
x1_ora = overRepresentation(x1,minsize=5,maxsize = 100,
                            correction = "BH",
                            GeneSetList = hsa_go,
                            foutput="x1_ora")
```
A `data.frame` with the results. 
```{r}
#| label: gse37
x1_ora_df = tidy(x1_ora)
```
A html report
```{r}
#| label: gse38
glimpse(x1_ora)
```
that can be opened with
```{r}
#| label: gse39
browseURL(glimpse(x1_ora))
```




## Gene set analysis 


```{r}
#| label: gse40
x1_self_mean = GeneSetTest(x = gse28619 ,y="type",
         test = rowt,association="pvalue",correction="BH",
         GeneNullDistr = "randomization",minsize = 5,
         GeneSetNullDistr = "self-contained",
         alternative="less",nmax = 1000,
         id = "ENTREZID",gsc=hsa_go,descriptive=mean,
         foutput = "x1_self_mean")
```

```{r}
#| label: gse41
x1_self_mean_df = tidy(x1_self_mean)
glimpse(x1_self_mean)
browseURL(glimpse(x1_self_mean))
```



