---
title: 'Supplementary code (update) for article: BgeeDB, an R package for retrieval of curated
  expression datasets and for gene list enrichment tests'
author: "Andrea Komljenović, Julien Roux, Julien Wollbrett, Marc Robinson-Rechavi and Frederic Bastian"
date: "June 11, 2018"
output: html_document
---

```{r, message = FALSE, warning = FALSE}
source("https://bioconductor.org/biocLite.R")
biocLite("BgeeDB")
# additional packages for installation
biocLite(c("edgeR", "Mfuzz", "biomaRt"))
```

``` {r load the package, message = FALSE, warning = FALSE}
library(BgeeDB)
```


``` {r list species}
# list the species
listBgeeSpecies()
```

```{r mouse affymetrix}
# specify species and data type
# the examples in this paper are based on Bgee release 14.0
# the following line targets the latest Bgee release. In order to target
# specifically the release 14.0, add the parameter 'release="14.0"'
bgee.affymetrix <- Bgee$new(species="Mus_musculus", dataType="affymetrix")

# retrieve annotation of all mouse affymetrix datasets in Bgee
annotation.bgee.mouse.affymetrix <- getAnnotation(bgee.affymetrix)

```

```{r list experiments}
# retrieve annotations of samples and experiments
sample.annotation <- annotation.bgee.mouse.affymetrix$sample.annotation
experiment.annotation <- annotation.bgee.mouse.affymetrix$experiment.annotation

# list experiments including a zygote sample
selected.experiments <- unique(sample.annotation$Experiment.ID[sample.annotation$Stage.ID == "UBERON:0000106"])
experiment.annotation[experiment.annotation$Experiment.ID %in% selected.experiments,]

```

```{r list sample}
# stages sampled in each of these experiments
unique(sample.annotation[sample.annotation$Experiment.ID %in% selected.experiments, c("Experiment.ID", "Stage.name")])
```


```{r sample annotation in experiment E-MEXP-51}
# List all samples from E-MEXP-51 in Bgee
sample.annotation[sample.annotation$Experiment.ID == "E-MEXP-51",]

```


```{r choose experiment}
data.E.MEXP.51 <- getData(bgee.affymetrix, experimentId="E-MEXP-51")
head(data.E.MEXP.51)

```

```{r format the data }
data.E.MEXP.51.formatted <- formatData(bgee.affymetrix,
                                       data.E.MEXP.51,
                                       callType="all",
                                       stats="intensities")
data.E.MEXP.51.formatted
```

```{r show expression }
# matrix of expression intensities
head(exprs(data.E.MEXP.51.formatted))
```

```{r show annotation }
# annotation of samples
pData(data.E.MEXP.51.formatted)
```


```{r show features }
# annotation of probesets
head(fData(data.E.MEXP.51.formatted))
```


```{r rna seq retrieval}
# specify species and data type
# the examples in this paper are based on Bgee release 14.0
# the following line targets the latest Bgee release. In order to target
# specifically the release 14.0, add the parameter 'release="14.0"'
bgee.rnaseq <- Bgee$new(species="Macaca_mulatta", dataType="rna_seq")


# retrieve annotations of RNA-seq samples and experiments
annotation.bgee.macaque.rna.seq <- getAnnotation(bgee.rnaseq)
sample.annotation <- annotation.bgee.macaque.rna.seq$sample.annotation
experiment.annotation <- annotation.bgee.macaque.rna.seq$experiment.annotation

# list experiments including both brain and liver samples
selected.experiments <- intersect(
    unique(
    sample.annotation$Experiment.ID[sample.annotation$Anatomical.entity.ID == "UBERON:0000955"]),
    unique(
    sample.annotation$Experiment.ID[sample.annotation$Anatomical.entity.ID == "UBERON:0002107"])
    )
experiment.annotation[experiment.annotation$Experiment.ID %in% selected.experiments,]

```


```{r sample annotation for macaca}
# check whether experiments include biological replicates
sample.annotation[sample.annotation$Experiment.ID %in% selected.experiments & (sample.annotation$Anatomical.entity.ID == "UBERON:0000955" | sample.annotation$Anatomical.entity.ID == "UBERON:0002107"), c("Experiment.ID",
                      "Library.ID",
                      "Anatomical.entity.ID",
                      "Anatomical.entity.name",
                      "Stage.ID")]
```


```{r one experiment}
data.GSE41637 <- getData(bgee.rnaseq, experimentId="GSE41637")
data.GSE41637.formatted <- formatData(bgee.rnaseq, data.GSE41637, callType="all", stats="counts")
data.GSE41637.formatted
```

```{r formats}
# for simplicity, keep only one sample per condition
data.E.MEXP.51.formatted <- data.E.MEXP.51.formatted[, !duplicated(pData(data.E.MEXP.51.formatted)[
    c("Anatomical.entity.ID",
      "Anatomical.entity.name",
      "Stage.ID", "Stage.name")])]

# order developmental stages
stages <- c("GVoocyte1","MIIoocyte1","Zygote1",
            "Early2-cell1","4Cell1","16cell1",
            "EarlyBlastocyst1","MidBlastocyst1","LateBlastocyst1")
data.E.MEXP.51.formatted <- data.E.MEXP.51.formatted[, stages]
# filter out rows with no variance
data.E.MEXP.51.formatted <-
data.E.MEXP.51.formatted [apply(exprs(data.E.MEXP.51.formatted), 1, sd) != 0, ]

```


```{r mfuzz analysis}
# Mfuzz clustering
library(Mfuzz)
# standardize matric of expression data
z.mat <- standardise(data.E.MEXP.51.formatted)
# cluster data into 16 clusters
clusters <- mfuzz(z.mat, centers=16, m=1.25)

```

```{r fig.width=7, fig.height=6}
# visualizing clusters
mfuzz.plot2(z.mat,
            cl=clusters,
            mfrow=c(4,4),
            colo="fancy",
            time.labels=row.names(pData(z.mat)),
            las=2, xlab="", ylab="Standardized expression level",
            x11 = FALSE)
```


```{r edger analysis}
# differential expression analysis with edgeR
library(edgeR)

# subset the dataset to brain and liver
brain.liver <- data.GSE41637.formatted[, pData(data.GSE41637.formatted)$Anatomical.entity.name %in% c("brain", "liver")]

# filter out very lowly expressed genes
brain.liver.filtered <- brain.liver[rowSums(cpm(brain.liver) > 1) > 3, ]

# create edgeR DGElist object
dge <- DGEList(counts= brain.liver.filtered, group=pData(brain.liver.filtered)$Anatomical.entity.name)
dge <- calcNormFactors(dge)
dge <- estimateCommonDisp(dge)
dge <- estimateTagwiseDisp(dge)
de <- exactTest(dge, pair=c("brain","liver"))
de.genes <- topTags(de, n=nrow(de))$table
```

```{r fig.width=7, fig.height=6}
# MA plot with DE genes highlighted
plotSmear(dge, de.tags=rownames(de.genes)[de.genes$FDR < 0.01], cex=0.3)
```


```{r get zebrafish anatomical structures}
# the examples in this paper are based on Bgee release 14.0
# the following line targets the latest Bgee release. In order to target
# specifically the release 14.0, add the parameter 'release="14.0"'
bgee.topanat <- Bgee$new(species="Danio_rerio")
top.anat.data <- loadTopAnatData(bgee.topanat)
```

```{r biomart annotation}
library(biomaRt)

# zebrafish data in Ensembl 84 (stable link)
ensembl <- useMart("ENSEMBL_MART_ENSEMBL",
                   dataset="drerio_gene_ensembl",
                   host="mar2016.archive.ensembl.org")

# get the mapping of Ensembl genes to phenotypes
genes.to.phenotypes <- getBM(filters=c("phenotype_source"),
                           value=c("ZFIN"),
                           attributes=c("ensembl_gene_id","phenotype_description"),
                           mart=ensembl)

# select phenotypes related to pectoral fin
phenotypes <- grep("pectoral fin", unique(genes.to.phenotypes$phenotype_description), value=T)

# select the genes annotated to select phenotypes
genes <- unique(genes.to.phenotypes$ensembl_gene_id[genes.to.phenotypes$phenotype_description
                                                    %in% phenotypes])

# prepare the gene list vector
gene.list <- factor(as.integer(unique(genes.to.phenotypes$ensembl_gene_id) %in% genes))
names(gene.list) <- unique(genes.to.phenotypes$ensembl_gene_id)
summary(gene.list)

```


```{r prepare topAnat}
# prepare the topAnat object based on topGO
top.anat.object <- topAnat(top.anat.data, gene.list)
top.anat.object
```


```{r run topAnat}
results <- runTest(top.anat.object, algorithm='weight', statistic='fisher')
results
```


```{r make table}
# retrieve anatomical structures enriched at a 1% FDR threshold, sorted by decreasing fold enrichment
table.over <- makeTable(top.anat.data, top.anat.object, results, cutoff=0.01)
head(table.over)
```


```{r retrieve significant terms}
# In order to retrieve significant genes mapped to the term " paired limb/fin bud"
term <- "UBERON:0004357"
termStat(top.anat.object, term)

# 198 genes mapped to this term for Bgee 14.0 and Ensembl 84
genesInTerm(top.anat.object, term)
# 48 significant genes mapped to this term for Bgee 14.0 and Ensembl 84
annotated <- genesInTerm(top.anat.object, term)[["UBERON:0004357"]]
annotated[annotated %in% sigGenes(top.anat.object)]
```


```{r}
sessionInfo()
```