Breakout Exercise 2 - Gene Annotations

Estimated time: 15 Minutes

Motivation

The data from from our differential expression test (result_plus_vs_minus) is an important result from our analysis.

head(results_plus_vs_minus)

# A tibble: 6 × 8
  id                 baseMean log2FoldChange lfcSE  stat   pvalue     padj call 
  <chr>                 <dbl>          <dbl> <dbl> <dbl>    <dbl>    <dbl> <fct>
1 ENSMUSG00000032715   1178.          -3.64  0.450 -8.10 5.72e-16 7.08e-12 Down 
2 ENSMUSG00000101645     77.3         10.1   1.33   7.58 3.42e-14 2.12e-10 Up   
3 ENSMUSG00000027313    243.          -2.67  0.410 -6.50 7.96e-11 2.47e- 7 Down 
4 ENSMUSG00000027889   3250.          -0.946 0.145 -6.52 7.00e-11 2.47e- 7 Down 
5 ENSMUSG00000054136     92.6         -3.39  0.526 -6.44 1.19e-10 2.95e- 7 Down 
6 ENSMUSG00000020142    625.          -1.72  0.287 -5.99 2.07e- 9 4.28e- 6 Down 

Looking at the first few rows, we notice that the genes are labelled with their ENSEMBL identifiers, e.g. ENSMUSG00000000056. It is quite rare that anyone will know off the top of their head which gene is associated with any given ENSEMBL identifier; we tend to think of genes by their symbols, e.g. Tp53.

Exercise

Using the biomaRt package, to map ENSEMBL IDs to gene symnbols, add a column with the gene symbols to results_plus_vs_minus.

Instructions

• One group member should share their screen in the breakout room. If nobody volunteers, a helper may randomly select someone.
• The group members should discuss the exercise and work together to find a solution.
• If there is time after a solution is found, allow all members to complete the exercise.

Preliminaries

First, let’s download the mapping of ENSEMBL IDs to gene symbols together. Then we’ll split into small groups to add the gene symbol column we want. We’ll start by loading the biomaRt library and calling the useEnsembl() function to select the database we’ll use to extract the information we need.

library('biomaRt')
ensembl = useEnsembl(dataset = 'mmusculus_gene_ensembl', biomart='ensembl')

Note that this process takes some time and will take up a larger amount of working memory so proceed with caution if you try to run these commands on a laptop with less than 4G of memory

To identify possible filters to restrict our data, we can use the listFilters function. To identify the attributes we want to retrive, we can use the listAttributes function. The best approach is to use list or search functions to help narrow down the available options.

head(listFilters(mart = ensembl), n = 20)
head(listAttributes(ensembl), n = 30)

We can access additional genomic annotations using the bioMart package. To identify we’ll structure our ‘query’ or search of the bioMart resources to use the ENSEMBL id from our alignment to add the gene symbols and gene description for each gene.

id_mapping = getBM(attributes=c('ensembl_gene_id', 'external_gene_name'),
      filters = 'ensembl_gene_id',
      values = row.names(assay(dds)),
      mart = ensembl)
# will take some time to run

# Preview the result
head(id_mapping)

     ensembl_gene_id external_gene_name
1 ENSMUSG00000000001              Gnai3
2 ENSMUSG00000000028              Cdc45
3 ENSMUSG00000000031                H19
4 ENSMUSG00000000037              Scml2
5 ENSMUSG00000000049               Apoh
6 ENSMUSG00000000056               Narf

Now that we have the ENSEMBL information and a gene symbol to match to our results, we can proceed in the smaller groups. As with the previous exercise, we have broken it into small steps with hints as needed.

Note: For additional information regarding bioMart, please consult the ENSEMBL bioMart vignette or the broader Bioconductor Annotation Resources vignette.

Steps

1. Look at the two data frames that are going to be needed for the exercise: id_mapping and results_plus_vs_minus.

Answer

head(id_mapping)

     ensembl_gene_id external_gene_name
1 ENSMUSG00000000001              Gnai3
2 ENSMUSG00000000028              Cdc45
3 ENSMUSG00000000031                H19
4 ENSMUSG00000000037              Scml2
5 ENSMUSG00000000049               Apoh
6 ENSMUSG00000000056               Narf

head(results_plus_vs_minus)

# A tibble: 6 × 8
  id                 baseMean log2FoldChange lfcSE  stat   pvalue     padj call 
  <chr>                 <dbl>          <dbl> <dbl> <dbl>    <dbl>    <dbl> <fct>
1 ENSMUSG00000032715   1178.          -3.64  0.450 -8.10 5.72e-16 7.08e-12 Down 
2 ENSMUSG00000101645     77.3         10.1   1.33   7.58 3.42e-14 2.12e-10 Up   
3 ENSMUSG00000027313    243.          -2.67  0.410 -6.50 7.96e-11 2.47e- 7 Down 
4 ENSMUSG00000027889   3250.          -0.946 0.145 -6.52 7.00e-11 2.47e- 7 Down 
5 ENSMUSG00000054136     92.6         -3.39  0.526 -6.44 1.19e-10 2.95e- 7 Down 
6 ENSMUSG00000020142    625.          -1.72  0.287 -5.99 2.07e- 9 4.28e- 6 Down 

2. Think about what column of results_plus_vs_minus we want to match to id_mapping, and then what column we want to extract values from in id_mapping.

Answer

We want to match the id column of results_plus_vs_minus to the ensembl_gene_id column of id_mapping, and once that match is found, we want to extract the external_gene_name column of id_mapping to get the gene symbol.

3. Look at the documentation for dplyr::left_join() and try to merge the id_mapping table into the results_plus_vs_minus table on the columns you identified above.

Answer

results_plus_vs_minus_annotated = results_plus_vs_minus %>%
    left_join(id_mapping, by = c('id' = 'ensembl_gene_id'))
head(results_plus_vs_minus_annotated)

# A tibble: 6 × 9
  id                 baseMean log2FoldChange lfcSE  stat   pvalue     padj call  external_gene_name
  <chr>                 <dbl>          <dbl> <dbl> <dbl>    <dbl>    <dbl> <fct> <chr>             
1 ENSMUSG00000032715   1178.          -3.64  0.450 -8.10 5.72e-16 7.08e-12 Down  Trib3             
2 ENSMUSG00000101645     77.3         10.1   1.33   7.58 3.42e-14 2.12e-10 Up    Gm28635           
3 ENSMUSG00000027313    243.          -2.67  0.410 -6.50 7.96e-11 2.47e- 7 Down  Chac1             
4 ENSMUSG00000027889   3250.          -0.946 0.145 -6.52 7.00e-11 2.47e- 7 Down  Ampd2             
5 ENSMUSG00000054136     92.6         -3.39  0.526 -6.44 1.19e-10 2.95e- 7 Down  Adm2              
6 ENSMUSG00000020142    625.          -1.72  0.287 -5.99 2.07e- 9 4.28e- 6 Down  Slc1a4            

4. Optionally, how could we use some of the tidyverse functions we’ve encountered previously to rename the external_gene_name column to symbol and to move it into the second column position? Hint: Because of the order of the packages we may have loaded, consider using dplyr::rename() and dplyr::select() instead of just the select() function. We can discuss why this is the case together.

Answer

results_plus_vs_minus_annotated = results_plus_vs_minus_annotated %>%
    dplyr::rename('symbol' = 'external_gene_name') %>%
    dplyr::select(id, symbol, everything())
results_plus_vs_minus_annotated

# A tibble: 16,249 × 9
   id                 symbol   baseMean log2FoldChange lfcSE  stat   pvalue     padj call 
   <chr>              <chr>       <dbl>          <dbl> <dbl> <dbl>    <dbl>    <dbl> <fct>
 1 ENSMUSG00000032715 Trib3      1178.          -3.64  0.450 -8.10 5.72e-16 7.08e-12 Down 
 2 ENSMUSG00000101645 Gm28635      77.3         10.1   1.33   7.58 3.42e-14 2.12e-10 Up   
 3 ENSMUSG00000027313 Chac1       243.          -2.67  0.410 -6.50 7.96e-11 2.47e- 7 Down 
 4 ENSMUSG00000027889 Ampd2      3250.          -0.946 0.145 -6.52 7.00e-11 2.47e- 7 Down 
 5 ENSMUSG00000054136 Adm2         92.6         -3.39  0.526 -6.44 1.19e-10 2.95e- 7 Down 
 6 ENSMUSG00000020142 Slc1a4      625.          -1.72  0.287 -5.99 2.07e- 9 4.28e- 6 Down 
 7 ENSMUSG00000040280 Ndufa4l2     85.9         -3.62  0.628 -5.75 8.68e- 9 1.54e- 5 Down 
 8 ENSMUSG00000023951 Vegfa       350.          -2.44  0.432 -5.64 1.67e- 8 2.59e- 5 Down 
 9 ENSMUSG00000027737 Slc7a11     119.          -2.32  0.429 -5.41 6.47e- 8 8.02e- 5 Down 
10 ENSMUSG00000115420 Rmrp         16.0         21.1   3.91   5.41 6.42e- 8 8.02e- 5 Up   
# … with 16,239 more rows
# ℹ Use `print(n = ...)` to see more rows