# Purpose of this analysis

This notebook illustrates one way that you can use microarray data from refine.bio to perform Principal Component Analysis (PCA) and plot the scores using `ggplot2`. 

⬇️ [**Jump to the analysis code**](#analysis) ⬇️

# How to run this example

For general information about our tutorials and the basic software packages you will need, please see our ['Getting Started' section](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#how-this-tutorial-is-structured).
We recommend taking a look at our [Resources for Learning R](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#resources-for-learning-r) if you have not written code in R before. 

## Obtain the `.Rmd` file

To run this example yourself, [download the `.Rmd` for this analysis by clicking this link](https://alexslemonade.github.io/refinebio-examples/02-microarray/dimension-reduction_microarray_01_pca.Rmd).

Clicking this link will most likely send this to your downloads folder on your computer. 
Move this `.Rmd` file to where you would like this example and its files to be stored.

You can open this `.Rmd` file in RStudio and follow the rest of these steps from there. (See our [section about getting started with R notebooks](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#how-to-get-and-use-rmds) if you are unfamiliar with `.Rmd` files.)

## Set up your analysis folders 

Good file organization is helpful for keeping your data analysis project on track!
We have set up some code that will automatically set up a folder structure for you. 
Run this next chunk to set up your folders! 

If you have trouble running this chunk, see our [introduction to using `.Rmd`s](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#how-to-get-and-use-rmds) for more resources and explanations. 

```{r}
# Create the data folder if it doesn't exist
if (!dir.exists("data")) {
  dir.create("data")
}

# Define the file path to the plots directory
plots_dir <- "plots"

# Create the plots folder if it doesn't exist
if (!dir.exists(plots_dir)) {
  dir.create(plots_dir)
}

# Define the file path to the results directory
results_dir <- "results"

# Create the results folder if it doesn't exist
if (!dir.exists(results_dir)) {
  dir.create(results_dir)
}
```

In the same place you put this `.Rmd` file, you should now have three new empty folders called `data`, `plots`, and `results`!

## Obtain the dataset from refine.bio

For general information about downloading data for these examples, see our ['Getting Started' section](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#how-to-get-the-data). 

Go to this [dataset's page on refine.bio](https://www.refine.bio/experiments/GSE37382).

Click the "Download Now" button on the right side of this screen.  

<img src="https://github.com/AlexsLemonade/refinebio-examples/raw/40e47f4d3f39283effbd9843a457168061be9680/template/screenshots/download-now.png" width=200>  

Fill out the pop up window with your email and our Terms and Conditions:  

<img src="https://github.com/AlexsLemonade/refinebio-examples/raw/40e47f4d3f39283effbd9843a457168061be9680/template/screenshots/download-email.png" width=500>    

It may take a few minutes for the dataset to process.
You will get an email when it is ready.

## About the dataset we are using for this example

For this example analysis, we will use this [medulloblastoma dataset](https://www.refine.bio/experiments/GSE37382).

The data that we downloaded from refine.bio for this analysis has 285 microarray samples obtained from patients with medulloblastoma.
The purpose of the experiment is to identify the predominant regions with somatic copy number aberrations in each medulloblastoma subgroup.
We will use the `subgroup` and `histology` variables from our metadata to annotate our plot later.
These variables hold the information on the one of the defined subgroups of medulloblastoma to which it sample belongs to and the molecular histology of each sample, respectively.


## Place the dataset in your new `data/` folder

refine.bio will send you a download button in the email when it is ready. 
Follow the prompt to download a zip file that has a name with a series of letters and numbers and ends in `.zip`. 
Double clicking should unzip this for you and create a folder of the same name.  

<img src="https://github.com/AlexsLemonade/refinebio-examples/raw/40e47f4d3f39283effbd9843a457168061be9680/template/screenshots/download-folder-structure.png" width=400> 

For more details on the contents of this folder see [these docs on refine.bio](http://docs.refine.bio/en/latest/main_text.html#downloadable-files).  

The `<experiment_accession_id>` folder has the data and metadata TSV files you will need for this example analysis.
Experiment accession ids usually look something like `GSE1235` or `SRP12345`. 

Copy and paste the `GSE37382` folder into your newly created `data/` folder.

## Check out our file structure!

Your new analysis folder should contain: 

- The example analysis `.Rmd` you downloaded  
- A folder called "data" which contains:
  - The `GSE37382` folder which contains:
    - The gene expression  
    - The metadata TSV  
- A folder for `plots` (currently empty)  
- A folder for `results` (currently empty)
    
Your example analysis folder should now look something like this (except with respective experiment accession ID and analysis notebook name you are using): 

<img src="https://github.com/AlexsLemonade/refinebio-examples/raw/40e47f4d3f39283effbd9843a457168061be9680/template/screenshots/analysis-folder-structure.png" width=400>

In order for our example here to run without a hitch, we need these files to be in these locations so we've constructed a test to check before we get started with the analysis. 
These chunks will declare your file paths and double check that your files are in the right place.

First we will declare our file paths to our data and metadata files, which should be in our data directory.
This is handy to do because if we want to switch the dataset (see next section for more on this) we are using for this analysis, we will only have to change the file path here to get started. 

```{r}
# Define the file path to the data directory
# Replace with the path of the folder the files will be in
data_dir <- file.path("data", "GSE37382")

# Declare the file path to the gene expression matrix file
# inside directory saved as `data_dir`
# Replace with the path to your dataset file
data_file <- file.path(data_dir, "GSE37382.tsv")

# Declare the file path to the metadata file
# inside the directory saved as `data_dir`
# Replace with the path to your metadata file
metadata_file <- file.path(data_dir, "metadata_GSE37382.tsv")
```

Now that our file paths are declared, we can use the `file.exists()` function to check that the files are where we specified above. 

```{r}
# Check if the gene expression matrix file is at the path stored in `data_file`
file.exists(data_file)

# Check if the metadata file is at the file path stored in `metadata_file`
file.exists(metadata_file)
```

If the chunk above printed out `FALSE` to either of those tests, you won't be able to run this analysis _as is_ until those files are in the appropriate place.

If the concept of a "file path" is unfamiliar to you; we recommend taking a look at our [section about file paths](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#an-important-note-about-file-paths-and-Rmds). 

# Using a different refine.bio dataset with this analysis?

If you'd like to adapt an example analysis to use a different dataset from [refine.bio](https://www.refine.bio/), we recommend placing the files in the `data/` directory you created and changing the filenames and paths in the notebook to match these files (we've put comments to signify where you would need to change the code).
We suggest saving plots and results to `plots/` and `results/` directories, respectively, as these are automatically created by the notebook.
From here you can customize this analysis example to fit your own scientific questions and preferences. 

***

<!-- Do not delete this line --> <a name="analysis" style="padding-top:56px;margin-top:-56px;">&nbsp;</a>


# PCA Visualization - Microarray

## Install libraries

See our Getting Started page with [instructions for package installation](https://alexslemonade.github.io/refinebio-examples/01-getting-started/getting-started.html#what-you-need-to-install) for a list of the other software you will need, as well as more tips and resources.

Attach the packages we need for this analysis:

```{r message=FALSE}
# Attach the library
library(ggplot2)

# We will need this so we can use the pipe: %>%
library(magrittr)
```

## Import and set up data

Data downloaded from refine.bio include a metadata tab separated values (TSV) file and a data TSV file. 
This chunk of code will read the both TSV files and add them as data frames to your environment.

We stored our file paths as objects named `metadata_file` and `data_file` in [this previous step](#check-out-our-file-structure).

```{r}
# Read in metadata TSV file
metadata <- readr::read_tsv(metadata_file)

# Read in data TSV file
df <- readr::read_tsv(data_file) %>%
  # Tuck away the gene ID column as row names, leaving only numeric values
  tibble::column_to_rownames("Gene")
```

Let's ensure that the metadata and data are in the same sample order. 

```{r}
# Make the sure the columns (samples) are in the same order as the metadata
df <- df %>%
  dplyr::select(metadata$geo_accession)

# Check if this is in the same order
all.equal(colnames(df), metadata$geo_accession)
```


Now we are going to use a combination of functions from base R and the `ggplot2` package to perform and visualize the results of the Principal Component Analysis (PCA) dimension reduction technique on our medulloblastoma samples.

## Perform Principal Components Analysis

In this code chunk, we are going to perform Principal Component Analysis (PCA) on our data and create a data frame using the PCA scores and the variables from our metadata that we are going to use to annotate our plot later.
We are using the base R `prcomp()` function to perform Principal Component Analysis (PCA) here.
The `prcomp()` function calculates principal component scores for each row of a matrix, but our data is arranged with each sample in a column, so we will need to transpose the data frame first.
In most cases, we will want to use the `scale = TRUE` argument so that all of the expression measurements have the same variance.
This prevents the PCA results from being dominated by a few highly variable genes.

```{r}
# Perform Principal Component Analysis (PCA) using the `prcomp()` function
pca <- prcomp(
  t(df), # transpose our data frame to obtain PC scores for samples, not genes
  scale = TRUE # we want the data scaled to have unit variance for each gene
)
```

Let's take a preview of the PCA results.
Each row will be a sample, as in the transposed data matrix we used as input, and each column is one of the new principal component (PC) values.
We are using indexing to only print out the first 5 PC columns: `[, 1:5]`.

```{r}
# We can access the results from our `pca` object using `$x`
head(pca$x[, 1:5])
```

In total, we have 285 principal component values, because we provided 285 samples' data (we will always have as many PCs as the smaller dimension of the input data matrix).

## Explore Variance in PCA Results

Before visualizing and interpreting the results, it can be useful to understand the proportion of variance explained by each principal component.
The principal components are automatically ordered by the variance they explain, meaning PC1 would always be the principal component that explains the most variance in your dataset.
If the largest variance component, PC1, explained 96% of the variance in your data and very clearly showed a difference between sample batches you would be very concerned about your dataset!
On the other hand, if a separation of batches was apparent in a different principal component that explained a low proportion of variance and the first few PCs explained most of the variance and appeared to correspond to something like tissue type and treatment, you would be less concerned [@CCDL2020-cluster-validation].

The `summary()` function reports the proportion of variance explained by each principal component.

```{r}
# Save the summary of the PCA results using the `summary()` function
pca_summary <- summary(pca)
```

The `importance` element of the summary object contains the proportion of variance explained by each principal component along with other statistics, with `pca_summary$importance`, we can use indexing to only look at the first `n` PCs.

```{r}
# Now access the importance information for the first 5 PCs
pca_summary$importance[, 1:5]
```

Now that we've seen the proportion of variance for the first set of PCs, let's prepare and plot the PC scores for the first two principal components, the components that explain the largest proportion of the expression variance in our dataset.
(Note though, that in this case, they explain less than 15% of the total variance!)

## Prepare a final data frame with PCA results for plotting

In the next chunk, we are going to extract the first two principal components from our `pca` object to prepare a data frame for plotting.

```{r}
# Make the first two PCs into a data frame for plotting with `ggplot2`
pca_df <- data.frame(pca$x[, 1:2]) %>%
  # Turn samples IDs stored as row names into a column
  tibble::rownames_to_column("refinebio_accession_code") %>%
  # Bring only the variables that we want from the metadata into this data frame
  # here we are going to join by `refinebio_accession_code` values
  dplyr::inner_join(
    dplyr::select(metadata, refinebio_accession_code, histology, subgroup),
    by = "refinebio_accession_code"
  )
```

## Plot PCA Results

Now let's plot the PC scores for the first two principal components.

Let's also label the data points based on their genotype subgroup since medulloblastoma has been found to comprise of subgroups that each have molecularly distinct profiles [@Northcott2012].

```{r}
# Make a scatterplot using `ggplot2` functionality
pca_plot <- ggplot(
  pca_df,
  aes(
    x = PC1,
    y = PC2,
    color = subgroup # label points with different colors for each `subgroup`
  )
) +
  geom_point() + # Plot individual points to make a scatterplot
  theme_classic() # Format as a classic-looking plot with no gridlines

# Print out the plot here
pca_plot
```

Looks like Group 4 and SHH groups cluster with each other somewhat, but Group 3 seems to be less distinct, as there are some samples clustering with Group 4 as well.
Most of the differences that we see between the groups are along the first axis of variation, PC1.

We can add another label to our plot to get more information about our dataset.
Let's also label the data points based on the histological subtype that each sample belongs to.

```{r}
# Make a scatterplot with ggplot2
pca_plot <- ggplot(
  pca_df,
  aes(
    x = PC1,
    y = PC2,
    color = subgroup, # Draw points with different colors for each `subgroup`
    shape = histology # Use a different shape for each `histology` group
  )
) +
  geom_point() +
  theme_classic()

# Print out the plot here
pca_plot
```

Adding the histological subtype label to our plot made our plot more informative, but the diffuse Group 3 data doesn't appear to be related to a histology subtype.
We could test out other variables as annotation labels to get a further understanding of the cluster behavior of each subgroup, or plot other PC values to see if they might also reveal some structure in the data.

## Save annotated PCA plot as a PNG

Now that we have an annotated PCA plot, let's save it!

You can easily switch this to save to a JPEG or TIFF by changing the file name within the `ggsave()` function to the respective file suffix.

```{r}
# Save plot using `ggsave()` function
ggsave(
  file.path(
    plots_dir,
    "GSE37382_pca_scatterplot.png" # Replace with a good file name for your plot
  ),
  plot = pca_plot # The plot object that we want saved to file
)
```

# Resources for further learning

- [Overall PCA Explanation by Matt Brems](https://towardsdatascience.com/a-one-stop-shop-for-principal-component-analysis-5582fb7e0a9c) [@Brems2017]
- [A visual explanation of PCA](http://setosa.io/ev/principal-component-analysis/) [@pca-visually-explained]
- [Guidelines on choosing dimension reduction methods](https://journals.plos.org/ploscompbiol/article/file?id=10.1371/journal.pcbi.1006907&type=printable) [@Nguyen2019]
- [More on `ggplot2`](http://r-statistics.co/Complete-Ggplot2-Tutorial-Part1-With-R-Code.html) [@Prabhakaran2016]

# Session info

At the end of every analysis, before saving your notebook, we recommend printing out your session info. 
This helps make your code more reproducible by recording what versions of software and packages you used to run this.

```{r}
# Print session info
sessioninfo::session_info()
```

# References
