# Purpose of this analysis

The purpose of this notebook is to provide an example of mapping gene IDs for microarray data obtained from refine.bio using `AnnotationDbi` packages [@Pages2020-package].

⬇️ [**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/gene-id-annotation_microarray_01_ensembl.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/GSE13490/cancer-stem-cells-are-enriched-in-the-side-population-cells-in-a-mouse-model-of-glioma).

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 [mouse glioma stem cells dataset](https://www.refine.bio/experiments/GSE13490/cancer-stem-cells-are-enriched-in-the-side-population-cells-in-a-mouse-model-of-glioma).

This dataset has 15 microarrays measuring gene expression in a transgenic mouse model of glioma.
The authors compared cells from side populations and non-side populations in both tumor samples and normal neural stem cells.


## 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 `GSE13490` 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 `GSE13490` 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", "GSE13490")

# 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, "GSE13490.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_GSE13490.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. 

refine.bio data comes with gene level data identified by Ensembl IDs.
Although this example notebook uses Ensembl IDs from Mouse, (<i>Mus musculus</i>), to obtain gene symbols this notebook can be easily converted for use with different species or annotation types <i>e.g.</i> protein IDs, gene ontology, accession numbers.

<b>For different species</b>, wherever the abbreviation `org.Mm.eg.db` or `Mm` is written, it must be replaced with the respective species abbreviation <i>e.g.</i> for <i>Homo sapiens</i> `org.Hs.eg.db` or `Hs` would be used.
In the case of our [RNA-seq gene identifier annotation example notebook](https://alexslemonade.github.io/refinebio-examples/03-rnaseq/gene-id-annotation_rnaseq_01_ensembl.html), a Zebrafish (<i>Danio rerio</i>) dataset is used, meaning `org.Dr.eg.db` or `Dr` would also need to be used there.
A full list of the annotation R packages from Bioconductor is at this [link](https://bioconductor.org/packages/release/BiocViews.html#___AnnotationData).

***

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


# Obtaining Annotation for Ensembl IDs - Microarray

refine.bio uses Ensembl IDs as the primary gene identifier in its data sets.
While this is a consistent and useful identifier, a string of apparently random letters and numbers is not the most user-friendly or informative for interpretation.
Luckily, we can use the Ensembl IDs that we have to obtain various different annotations at the gene/transcript level.
Let's get ready to use the Ensembl IDs from our mouse dataset to obtain the associated gene symbols.

## 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.

In this analysis, we will be using the `org.Mm.eg.db` R package [@Carlson2019-mouse], which is part of the Bioconductor `AnnotationDbi` framework [@Pages2020-package].
Bioconductor compiles annotations from various sources, and these packages provide convenient methods to access and translate among those annotations.
[Other species can be used](#using-a-different-refinebio-dataset-with-this-analysis).

```{r}
# Install the mouse annotation package
if (!("org.Mm.eg.db" %in% installed.packages())) {
  # Install this package if it isn't installed yet
  BiocManager::install("org.Mm.eg.db", update = FALSE)
}
```

Attach the packages we need for this analysis.
Note that attaching `org.Mm.eg.db` will automatically also attach `AnnotationDbi`.

```{r message=FALSE}
# Attach the library
library(org.Mm.eg.db)

# 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
expression_df <- readr::read_tsv(data_file) %>%
  # Tuck away the Gene ID column as row names
  tibble::column_to_rownames("Gene")
```

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

```{r}
# Make the data in the order of the metadata
expression_df <- expression_df %>%
  dplyr::select(metadata$geo_accession)

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

# Bring back the "Gene" column in preparation for mapping
expression_df <- expression_df %>%
  tibble::rownames_to_column("Gene")
```


## Map Ensembl IDs to gene symbols

The main work of translating among annotations will be done with the the `AnnotationDbi` function `mapIds()`.
The `mapIds()` function has a `multiVals` argument which denotes what to do when there are multiple mapped values for a single gene identifier.
The default behavior is to return just the first mapped value.
It is good to keep in mind that various downstream analyses may benefit from varied strategies at this step. 
Use `?mapIds` to see more options or strategies.

In the next chunk, we will run the `mapIds()` function and supply the `multiVals` argument with the `"list"` option in order to get a large list with all the mapped values found for each gene identifier.

```{r}
# Map ensembl IDs to their associated gene symbols
mapped_list <- mapIds(
  org.Mm.eg.db, # Replace with annotation package for your organism
  keys = expression_df$Gene,
  keytype = "ENSEMBL", # Replace with the gene identifiers used in your data
  column = "SYMBOL", # The type of gene identifiers you would like to map to
  multiVals = "list"
)
```

## Explore gene ID conversion

Now, let's take a look at our mapped object to see how the mapping went.

```{r}
# Let's use the `head()` function for a preview of our mapped list
head(mapped_list)
```

It looks like we have gene symbols that were successfully mapped to the Ensembl IDs we provided.
However, the data is now in a `list` object, making it a little more difficult to explore.
We are going to turn our list object into a data frame object in the next chunk.

```{r}
# Let's make our list a bit more manageable by turning it into a data frame
mapped_df <- mapped_list %>%
  tibble::enframe(name = "Ensembl", value = "Symbol") %>%
  # enframe() makes a `list` column; we will simplify it with unnest()
  # This will result in a data frame with one row per list item
  tidyr::unnest(cols = Symbol)
```

Now let's take a peek at our data frame.

```{r}
head(mapped_df)
```

We can see that our data frame has a new column `Symbol`.
Let's get a summary of the gene symbols in the `Symbol` column of our mapped data frame.

```{r}
# Use the `summary()` function to show the distribution of Symbol values
# We need to use `as.factor()` here to get the counts of unique values
# `maxsum = 10` limits the summary to 10 distinct values
summary(as.factor(mapped_df$Symbol), maxsum = 10)
```

There are 942 `NA`s in the `Symbol` column, which means that 942 out of the 17918 Ensembl IDs did not map to gene symbols.
942 out of 17918 is not too bad a rate, in our opinion, but note that different gene identifier types will have different mapping rates and that is to be expected.
Regardless, it is always good to be aware of how many genes you are potentially "losing" if you rely on this new gene identifier you've mapped to for downstream analyses.

However, if you have almost all `NA`s it is possible that the function was executed incorrectly or you may want to consider using a different gene identifier, if possible.

Now let's check to see if we have any genes that were mapped to multiple symbols.

```{r}
multi_mapped <- mapped_df %>%
  # Let's count the number of times each Ensembl ID appears in `Ensembl` column
  dplyr::count(Ensembl, name = "gene_symbol_count") %>%
  # Arrange by the genes with the highest number of symbols mapped
  dplyr::arrange(desc(gene_symbol_count)) %>%
  # Filter to include only the rows with multi mappings
  dplyr::filter(gene_symbol_count > 1)

# Let's look at the first 6 rows of our `multi_mapped` object
head(multi_mapped)
```

Looks like we have some cases where 3 gene symbols mapped to a single Ensembl ID.
We have a total of 130 out of 17984 Ensembl IDs with multiple mappings to gene symbols.
If we are not too worried about the 130 IDs with multiple mappings, we can filter them out for the purpose of having 1:1 mappings for our downstream analysis.

## Map Ensembl IDs to gene symbols -- filtering out multi mappings

The next code chunk we will rerun the `mapIds()` function, this time supplying the `"filter"` option to the `multiVals` argument.
This will remove all instances of multiple mappings and return a list of only the gene identifiers and symbols that had 1:1 mapping.
Use `?mapIds` to see more options or strategies.

```{r}
# Map Ensembl IDs to their associated gene symbols
filtered_mapped_df <- data.frame(
  "Symbol" = mapIds(
    org.Mm.eg.db, # Replace with annotation package for your organism
    keys = expression_df$Gene,
    keytype = "ENSEMBL", # Replace with the gene identifiers used in your data
    column = "SYMBOL", # The type of gene identifiers you would like to map to
    multiVals = "filter" # This will drop any genes that have multiple matches
  )
) %>%
  # Make an `Ensembl` column to store the rownames
  tibble::rownames_to_column("Ensembl") %>%
  # Join the remaining data from `expression_df` using the Ensembl IDs
  dplyr::inner_join(expression_df, by = c("Ensembl" = "Gene"))
```

Now, let's write our filtered and mapped results to file!

## Write mapped results to file

```{r Write to file}
# Write mapped and annotated data frame to output file
readr::write_tsv(filtered_mapped_df, file.path(
  results_dir,
  "GSE13490_Gene_Symbols.tsv" # Replace with a relevant output file name
))
```


# Resources for further learning

- Marc Carlson has prepared a nice [Introduction to Bioconductor Annotation Packages](https://bioconductor.org/packages/release/bioc/vignettes/AnnotationDbi/inst/doc/IntroToAnnotationPackages.pdf) [@Carlson2020-vignette]
- See our [RNA-seq gene ID conversion notebook](https://alexslemonade.github.io/refinebio-examples/03-rnaseq/gene-id-annotation_rnaseq_01_ensembl.html) as another applicable example, since the steps for this workflow do not change with technology.

# 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
