---
title: 'Lab #7: Unsupervised learning'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
set.seed(24)
```

We will analyze single-cell RNA-seq data from mouse embryonic development from Deng. et al. Science 2014, Vol. 343 no. 6167 pp. 193-196, “Single-Cell RNA-Seq Reveals Dynamic, Random Monoallelic Gene Expression in Mammalian Cells.” Cells were dissociated from oocyte to blastocyst stages of mouse embryos. The cell types, from earliest to latest, are MII oocyte, zygote, early 2-cell, mid 2-cell, late 2-cell, 4-cell, 8-cell, 16-cell, early blastocyst, mid blastocyst, and late blastocyst. Gene expression measurements are provided per cell.


1. Read in the RNA-seq dataset. What are the rows? What are the columns?
```{r}
```

2. We will analyze the log counts. That is, add one to each count and take the base 2 log.
```{r}
```

3. Each cell is associated with a cell stage, which is given in the column name. The cell stage is the string before the underscore. Extract the cell stages.
```{r}
```


## Principal Components Analysis

4. We will perform PCA to understand the major axes of variation in the gene expression values across cells. To do this, create a new dataframe that is the transpose of the original data, where rows are observations and columns are variables.
```{r}
```

5. To run PCA, remove any variables that are constant in their values across all observations.
```{r}
```

6. Run PCA using `prcomp` and name the result as `pr.out`. Set argument `scale=TRUE` to make sure the data will be normalized before performing PCA.
```{r}
```

7. Look at the `center` and `scale` components of `pr.out`. These correspond to the means and standard deviations of the variables that were used for scaling prior to implementing PCA. To visualize these values, plot a histogram.
```{r}
```

8. Look at the `x` component of `pr.out`. This contains the PC scores for each observation. What are the values for PC1 and PC2 for the first observation?
```{r}
```

9. Plot the first two principal components.
```{r}
```

10. Plot all pairwise combinations of the first five principal components as scatterplots (hint: use the `pairs` command).
```{r}
```

11. Now look at the `rotation` matrix component of `pr.out` which gives the principal component loading vector. How many PCs are there? Does that agree with what you would expect?
```{r}
```

12. Look at the `sdev` component of `pr.out`. This gives the standard deviation of each PC. Generate a vector `pr.var` that is the variance explained by each PC.
```{r}
```

13. Generate a new vector `pve` that is the proportion of variance explained by each vector. Generate the `sum` of `pve` as a confirmation that you are on the right lines.
```{r}
```

14. Generate plots of proportion of variance explained and cumulative variance explained against PC number. You will need the function `cumsum` which compute the cumulative sum of elements in a numeric vector.
```{r}
```

15. How many PCs do you need to explain at least 80% of the variance?
```{r}
```

## Hierarchical clustering

We will now create a heatmap of the gene expression values that clusters the cells by their similarity and the top genes by their similarity.

16. Clustering of cells is often based on pairwise correlations when handling gene expression data. To calculate pairwise correlations quickly, we'll need to install the "WGCNA" package. To install this, make sure you have installed the following packages from Bioconductor as well: "impute", "preprocessCore", "GO.db".
```{r}
```

17. Calculate the pairwise correlations and create a distance matrix using (1 - correlation). Use the `cor` function in WGCNA. You'll need to cast the matrix as a distance matrix using the `as.dist` function.
```{r}
```

18. Generate hierarchical clustering with `method="complete"`. Use the command `hc.complete <- hclust(distx, method="complete")`. Plot the tree.
```{r}
```

19. Based on the tree above, cluster the observations into 5 clusters. Use the `cutree` method to get the cluster assignments. How many observations are in each cluster? Was hierarchical clustering able to recover the cell stages?
```{r}
```

20. Calculate the variance of the expression values for each gene. Grab the ten genes whose expression values have the highest variance.
```{r}
```

21. Create a dendrogram between the genes using the same approach we used for the cells.
```{r}
```

22. Create a heatmap that plots the gene expression values for the most variable genes, where the cells and genes are clustered using the constructed dendrograms. Use the `heatmap.2` function in the `gplots` package.
```{r}
```