```{r global options, include = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

## Loading gene expression datasets from GEO

> The Gene Expression Omnibus (GEO) project was initiated at NCBI in 1999 in response to the growing demand for a public repository for data generated from high-throughput microarray experiments. GEO has a flexible and open design that allows the submission, storage, and retrieval of many types of data sets, such as those from high-throughput gene expression, genomic hybridization, and antibody array experiments. 

Currently, GEO contains all sorts of high-throughput experiments, although microarray and RNA-seq are the most popular. Explore and search on https://www.ncbi.nlm.nih.gov/geo/

While you can manually download and load the text file into R, it is easier to use the API/R package called GEOquery to do this automatically. To download the dataset on GEO directly into R, install these packages, and their dependencies:

```{r load_hidden, echo=FALSE, results="hide", warning=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("Biobase")
BiocManager::install("GEOquery")
BiocManager::install("limma")
BiocManager::install("edge")
BiocManager::install("genefilter")
BiocManager::install("qvalue")

suppressPackageStartupMessages({
  library(devtools)
  library(Biobase)
  library(limma)
  library(edge)
  library(genefilter)
  library(qvalue)
  library(GEOquery)
})
```

We are interested in the following dataset called GDS39 coming from the classic experiments on cell cycle regulation in yeast.
> Spellman PT, Sherlock G, Zhang MQ, Iyer VR et al. Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol Biol Cell 1998 Dec;9(12):3273-97. PMID: 9843569
Explore GEO and download directly from R:

The following code is used to download and look at meta data and gene expression data:

```{r geo}
#Loading a GDS file with GEOquery
dat <- getGEO('GDS39', destdir=".")

# look at meta data
Meta(dat)

# look at a small portion of gene expression 
Table(dat)[1:10,1:5]
colnames(Table(dat))
```

We look at one time point (one column) of this dataset. Particularly, we want to see the summary statistics. Then, a histogram, a density plot, and a boxplot are made. Please try on other columns, and multiple columns:

```{r vis1}
## look at one column (first tp) GSM940
GSM940 = as.numeric(Table(dat)[,3])
summary(GSM940)

hist(GSM940)
hist(GSM940, breaks=100, xlim=c(-5,5), main="Histogram of 1st Time Point in the Yeast Study")

## density estimation and visualization
density.est = density(GSM940, kernel="gaussian", na.rm=T)
plot(density.est, xlim=c(-5,5))

## look at one column (first tp) GSM940
par(mfrow=c(1,2))
boxplot(GSM940, ylim=c(-5,5))
```

## Tidy data and ggplot2

We will re-structure this gene expression data into a tidy data. For further and well-written summary of tidy data, please see http://vita.had.co.nz/papers/tidy-data.html
Using this tidy data, ggplot2 can be used to create some beautiful visualizations fast. You need to install necessary R packages, which are packaged into tidyverse: https://tidyverse.tidyverse.org/:

```{r ggplot2}
install.packages("tidyverse")
library(tidyverse)

geneexp <- Table(dat)
geneexp.tidy <- gather(geneexp, key="Samples", value="GeneExp", -c(1,2))

# let's make a box plot for all samples combined
ggplot(geneexp.tidy, aes(x=Samples,y=GeneExp)) + geom_boxplot() + ylim(-5,5)

# how about a violin plot for each samples (e.g., time points)
# a violin plot:  https://en.wikipedia.org/wiki/Violin_plot
ggplot(geneexp.tidy, aes(x=Samples,y=GeneExp)) + geom_violin() + ylim(-5,5) 

# how about a series of histogram?
ggplot(geneexp.tidy, aes(x=GeneExp)) + geom_histogram(bins=30) + xlim(-5,5) + facet_grid(~Samples)

# save ggplot, especially if your monitor is too small to see many facets
# use ggsave()
```

The tidyverse is an extensive network of packages that use common data structures and lingos to make exploration and visualization simple. It is a learning process that will benefit you in a long run. There are a lot of intuitive and helpful documentations online. 
See the cheatsheet on data import and structure manipulation: https://github.com/rstudio/cheatsheets/blob/master/data-import.pdf
Or the datacamp has an interactive R course on tidyverse: https://www.datacamp.com/courses/introduction-to-the-tidyverse

## Making a heatmap using gplots

A heatmap represents the individual values in a matrix as colors. For us, it usually allow the large table of numeric values to be visualized. Often, clustering or other organization of rows and columns help reveal patterns. Look and clean up any missing values. For this, we will simply remove any row that has missing values.

```{r missing}
dat.geneexp <- Table(dat)[,-c(1,2)]
dat.geneexp <- as.data.frame(lapply(dat.geneexp, as.numeric))
dat.geneexp <- apply(dat.geneexp,2,as.numeric)

# see if it has any missing values
sum(is.na(dat.geneexp))

# find rows with missing values
rows_missing <- apply(dat.geneexp, 1, function(x) any(is.na(x)))
dat.geneexp.complete <- dat.geneexp[!rows_missing,]
```

Let's create a heatmap with ``gplots``.

```{r heatmap}
install.packages(c("RColorBrewer","gplots"))
library(RColorBrewer)
library(gplots)
my_palette <- colorRampPalette(c("blue", "white", "yellow"))(n = 299)

pdf(file="geneexp_heatmap.pdf", width=7, height=7)
heatmap.2(as.matrix(dat.geneexp.complete),
          main = "Heatmap", # heat map title
          notecol="black",      # change font color of cell labels to black
          density.info="none",  # turns off density plot inside color legend
          trace="none",         # turns off trace lines inside the heat map
          margins =c(12,9),     # widens margins around plot
          col=my_palette,       # use on color palette defined earlier 
          dendrogram="none",     # only draw a row dendrogram
          scale = "row",
          Colv=FALSE)
dev.off()


pdf(file="geneexp_heatmap_nocluster.pdf", width=7, height=7)
heatmap.2(as.matrix(dat.geneexp.complete),
          main = "Heatmap", # heat map title
          notecol="black",      # change font color of cell labels to black
          density.info="none",  # turns off density plot inside color legend
          trace="none",         # turns off trace lines inside the heat map
          margins =c(12,9),     # widens margins around plot
          col=my_palette,       # use on color palette defined earlier 
          dendrogram="none",     # only draw a row dendrogram
          scale = "row",
          Colv=FALSE,
          Rowv=FALSE)
dev.off()
```

Please try creating a heatmap using ggplot2.

## Using Biobase and the ExpressionSet class

> Biobase contains standardized data structures to represent genomic data.  TheExpressionSetclass is designed to combine several different sources of information into a single convenient structure. An ExpressionSet can be manipulated (e.g., subsetted, copied) conveniently, and is the inputor output from many Bioconductor functions

The ExpressionSet class is complex, but can be very handy when several different data types must be combined, stored, and linked. These includes:

expression data in ``assayData``
metadata describing samples in ``phenoData``
metadata about the features in ``featureData``
annotations in ``annotation``
protocols in ``protocolData``
description of experiments in ``experimentData``

For a starter, let's create an ExpressionSet from GDS39 that we downloaded earlier from GEO:

```{r ExpressionSet}
#turning into ExpressionSet
eset <- GDS2eSet(dat, do.log2=FALSE) #log is already done
eset
#geneNames(eset)[1:10]
sampleNames(eset)
#pdat(eset)

#Loading a GPL (Annotation) file with GEOquery
Meta(dat)$platform
platform <- getGEO(Meta(dat)$platform, destdir=".")

Meta(platform)$title
colnames(Table(platform))
Table(platform)[1:10,1:4]
```
