RNA-seq Analysis in R

Annotation and Visualisation of RNA-seq results

Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law, Stephane Ballereau, Oscar Rueda, Ashley Sawle Based on the course RNAseq analysis in R delivered on May 11/12th 2016

Resources and data files

This material has been created using the following resources:

• http://monashbioinformaticsplatform.github.io/RNAseq-DE-analysis-with-R/99-RNAseq_DE_analysis_with_R.html
  
• http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
• https://bioconductor.github.io/BiocWorkshops/rna-seq-data-analysis-with-deseq2.html

Before starting this section, we will make sure we have all the relevant objects from the Differential Expression analysis present.

suppressPackageStartupMessages(library(DESeq2))
de_status <- readRDS("Robjects/de_status.rds")
de_gleason <- readRDS("Robjects/de_gleason.rds")
dds <- readRDS("Robjects/dds.rds")

Overview

• Visualising DE results
• Customising plots with ggplot2
• Filtering and exporting tables with dplyr
• Getting annotation
• Retrieving gene models

We can now have a list of genes ordered according to their evidence for being differentially-expressed.

library(dplyr)
library(tibble)
results_status <- as.data.frame(results(de_status)) %>% 
  rownames_to_column("GeneID")
  
results_ordered <- arrange(results_status, padj)

In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.

The log2 fold change for a particular comparison is plotted on the y-axis and the average of the counts normalized by size factor is shown on the x-axis (“M” for minus, because a log ratio is equal to log minus log, and “A” for average). Each gene is represented with a dot. Genes with an adjusted p value below a threshold (here 0.1, the default) are shown in red.

Both limma and DESeq2 have a function called plotMA, and R can sometimes pick the wrong function. To explictly use the DESeq2 function you can use:-

DESeq2::plotMA(de_status)

MA-plots often display a fanning-effect at the left-hand side (genes with low numbers of counts) due to the high variability of the measurements for these genes. For more informative visualization and more accurate ranking of genes by effect size (the log fold change may sometimes be referred to as an effect size), the DESeq2 authors recommend “shrinking” the log fold-changes which is available in DESeq2’s lfcShrink function. This results in more stable fold change values. The p-values are unaffected.

de_shrunk <- lfcShrink(de_status,contrast = c("Status","Tumour","Normal"))
DESeq2::plotMA(de_shrunk)

Another common plot for displaying the results of a differential expression analysis is a volcano plot

results_ordered %>% 
  ggplot(aes(x =log2FoldChange, y = -log10(padj))) + geom_point()

It can also be useful to examine the counts of reads for a single gene across the groups. A simple function for making this plot is plotCounts, which normalizes counts by sequencing depth and adds a pseudocount of 1/2 to allow for log scale plotting. The counts are grouped by the variables in intgroup, where more than one variable can be specified. Here we specify the gene which had the smallest p value from the results table created above. You can select the gene to plot by rowname or by numeric index:-

plotCounts(dds, "MIR4508",intgroup = c("Status"))

If we want greater control over how to visualise the data, we can use the plotCounts function to return the count data, but not actually produce the plot:-

plotCounts(dds, "MIR4508",intgroup = c("Status"),returnData=TRUE)

                                  count Status
TCGA.CH.5753.01A.11R.1580.07  61.435388 Tumour
TCGA.CH.5761.01A.11R.1580.07 204.621241 Tumour
TCGA.EJ.5507.01A.01R.1580.07 183.095180 Tumour
TCGA.FC.A8O0.01A.41R.A37L.07 479.611539 Tumour
TCGA.G9.6371.01A.11R.1789.07  91.793648 Tumour
TCGA.G9.6499.01A.12R.1965.07 328.732132 Tumour
TCGA.H9.A6BX.01A.31R.A311.07  89.189164 Tumour
TCGA.HC.7077.01A.11R.1965.07 202.857879 Tumour
TCGA.HC.7081.01A.11R.1965.07 301.699014 Tumour
TCGA.HC.7209.01A.11R.2118.07 212.938435 Tumour
TCGA.HC.7748.01A.11R.2118.07 126.892571 Tumour
TCGA.HC.A631.01A.11R.A29R.07 228.360080 Tumour
TCGA.HC.A632.01A.11R.A29R.07 123.860760 Tumour
TCGA.J4.A83K.01A.11R.A352.07 694.901430 Tumour
TCGA.KK.A6E8.01A.11R.A31N.07 259.133677 Tumour
TCGA.QU.A6IO.01A.11R.A31N.07 751.216078 Tumour
TCGA.QU.A6IP.01A.11R.A31N.07 153.528481 Tumour
TCGA.VP.A87H.01A.11R.A352.07 559.539652 Tumour
TCGA.YL.A8SJ.01B.11R.A37L.07 154.172447 Tumour
TCGA.ZG.A8QX.01A.11R.A37L.07 264.086903 Tumour
TCGA.CH.5761.11A.01R.1580.07   2.761284 Normal
TCGA.EJ.7314.11A.01R.2118.07   1.599503 Normal
TCGA.EJ.7327.11A.01R.2118.07   2.083821 Normal
TCGA.G9.6356.11A.01R.1789.07   6.504027 Normal
TCGA.G9.6384.11A.01R.1858.07   0.500000 Normal
TCGA.G9.6499.11A.02R.1965.07   6.301939 Normal
TCGA.HC.8260.11A.01R.2263.07   1.446933 Normal

Challenge 1

1. Use the option returnData=TRUE to get a data frame containing the counts of MIR4508 in tumours or normals. Visualise these data using ggplot2 (see plot A below).
2. Repeat the volcano plot from above, but use a different colour to indicate which genes are significant with an adjusted p-value less than 0.05. See plot B below
3. (Optional) The argument intgroup= can be used to retrieve and plot data from multiple variables of interest in the data. Use the value intgroup=c("Status","Batch") and compare the counts between status and batches. See plot C below.

Adding annotation to the DESeq2 results

There are a number of ways to add annotation, but we will demonstrate how to do this using the org.Hs.eg.db package. This package is one of several organism-level packages which are re-built every 6 months. These packages are listed on the annotation section of the Bioconductor, and are installed in the same way as regular Bioconductor packages. An alternative approach is to use biomaRt, an interface to the BioMart resource. BioMart is much more comprehensive, but the organism packages fit better into the Bioconductor workflow.

### Only execute when you need to install the package
library(BiocManager)
install("org.Hs.eg.db")

The packages are larger in size that Bioconductor software pacakges, but essentially they are databases that can be used to make offline queries.

library(org.Hs.eg.db)

First we need to decide what information we want. In order to see what we can extract we can run the columns function on the annotation database.

columns(org.Hs.eg.db)

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"     
 [4] "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
 [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
[10] "GENENAME"     "GO"           "GOALL"       
[13] "IPI"          "MAP"          "OMIM"        
[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"        
[19] "PFAM"         "PMID"         "PROSITE"     
[22] "REFSEQ"       "SYMBOL"       "UCSCKG"      
[25] "UNIGENE"      "UNIPROT"     

We are going to filter the database by a key or set of keys in order to extract the information we want. Valid names for the key can be retrieved with the keytypes function.

keytypes(org.Hs.eg.db)

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"     
 [4] "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
 [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
[10] "GENENAME"     "GO"           "GOALL"       
[13] "IPI"          "MAP"          "OMIM"        
[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"        
[19] "PFAM"         "PMID"         "PROSITE"     
[22] "REFSEQ"       "SYMBOL"       "UCSCKG"      
[25] "UNIGENE"      "UNIPROT"     

We should see SYMBOL, which is the type of key we are going to use in this case. If we are unsure what values are acceptable for the key, we can check what keys are valid with keys

keys(org.Hs.eg.db, keytype="SYMBOL")[1:10]

 [1] "A1BG"     "A2M"      "A2MP1"    "NAT1"    
 [5] "NAT2"     "NATP"     "SERPINA3" "AADAC"   
 [9] "AAMP"     "AANAT"   

For the top gene in our analysis the call to the function would be:-

select(org.Hs.eg.db, keys="MIR4508",
       keytype = "SYMBOL",columns=c("SYMBOL","GENENAME")
)


### In case of errors, try
### AnnotationDBI::select(org.Hs.eg.db, keys="MIR4508",
###       keytype = "SYMBOL",columns=c("SYMBOL","GENENAME")
###)

To annotate our results, we definitely want gene symbols and perhaps the full gene name. Let’s build up our annotation information into a new data frame using the select function.

anno <- AnnotationDbi::select(org.Hs.eg.db,keys=results_ordered$GeneID,
              columns=c("SYMBOL","GENENAME"),
              keytype="SYMBOL")
# Have a look at the annotation
head(anno)

   SYMBOL
1 MIR4508
2  SEMA5B
3   HSPA6
4  NKX2-3
5   HOXC4
6   HOXC5
                                      GENENAME
1                                microRNA 4508
2                                semaphorin 5B
3 heat shock protein family A (Hsp70) member 6
4                               NK2 homeobox 3
5                                  homeobox C4
6                                  homeobox C5

However, we have a problem that the resulting data frame has more rows than our results table. This is due to the one-to-many relationships that often occur when mapping between various identifiers.

dim(anno)

[1] 23372     2

dim(results_ordered)

[1] 23368     7

Such duplicated entries can be identified using the duplicated function.

dup_ids <- anno$SYMBOL[duplicated(anno$SYMBOL)]
filter(anno, SYMBOL %in% dup_ids) %>% 
  arrange(SYMBOL)

  SYMBOL
1    HBD
2    HBD
3  MEMO1
4  MEMO1
5   MMD2
6   MMD2
7    TEC
8    TEC
                                             GENENAME
1                            hemoglobin subunit delta
2                       hypophosphatemic bone disease
3                Methylation modifier for class I HLA
4                         mediator of cell motility 1
5 monocyte to macrophage differentiation associated 2
6                        Miyoshi muscular dystrophy 2
7                         tec protein tyrosine kinase
8           transient erythroblastopenia of childhood

Fortunately, there are not too many so hopefully we won’t lose too much information if we discard the entries that are duplicated. The first occurence of the duplicated ID will still be included in the table.

anno <- AnnotationDbi::select(org.Hs.eg.db,keys=results_ordered$GeneID,
              columns=c("ENSEMBL","SYMBOL","GENENAME","ENTREZID"),
              keytype="SYMBOL") %>% 
  filter(!duplicated(SYMBOL))
dim(anno)

[1] 23368     4

We can bind in the annotation information to the results data frame.

anno <- dplyr::rename(anno, GeneID = SYMBOL)
results_annotated <- left_join(results_ordered, anno,by="GeneID")

We can save the results table using the write.csv function, which writes the results out to a csv file that you can open in excel.

## create a directory for the results. don't give a warning if it already exists
dir.create("results",showWarnings = FALSE)
write.csv(results_annotated,file="tumour_vs_normal_DESeq_annotated.csv",row.names=FALSE)

We have already seen the use of a heatmap as a quality assessment tool to visualise the relationship between samples in an experiment. Another common use-case for such a plot is to visualise the results of a differential expression analysis.

Here we will take the top 10 genes from the differential expression analysis and produce a heatmap. The default colour palette goes from low expression in blue to high expression in red, which is a good alternative to the traditional red/green heatmaps which are not suitable for those with forms of colour-blindness.

library(pheatmap)
top_genes <- results_ordered$GeneID[1:10]
vsd <- vst(dds)
pheatmap(assay(vsd)[top_genes,])

The heatmap is more informative if we add colours underneath the sample dendrogram to indicate which sample group each sample belongs to. This we can do by creating a data frame containing metadata for each of the samples in our dataset. With the DESeq2 workflow we have already created such a data frame. We have to make sure the the rownames of the data frame are the same as the column names of the counts matrix.

sampleInfo <- as.data.frame(colData(dds)[,c("Status","gleason_score")])
pheatmap(assay(vsd)[top_genes,],
         annotation_col = sampleInfo)

Any plot we create in RStudio can be saved as a png or pdf file. We use the png or pdf function to create a file for the plot to be saved into and run the rest of the code as normal. The plot does not get displayed in RStudio, but printed to the specified file.

png("heatmap_top10_genes.png",width=800,height=800)
pheatmap(assay(vsd)[top_genes,],
         annotation_col = sampleInfo)
# dev.off()

Challenge 2

1. Repeat the same heatmap as above, but for the top 75 most differentially-expressed genes in the contrast between Gleason grades 9 and 6
2. Save the plot to a pdf file

Annotating plots with gene names

Now that we have an annotated table of results, we can add the gene names to some of the other plots we have created. This should be straightforward as ggplot2 has a label aesthetic that can be mapped to columns in a data frame. The geom_text plot will then display the labels. However, the following plot is a bit crowded.

## Not a good idea to run this!!
results_annotated %>% 
  ggplot(aes(x = log2FoldChange, y = -log10(padj), label=GeneID)) + geom_point() + geom_text()

The problem here is that ggplot2 is trying to label every point with a name; not quite what we want. The trick is to create a label that is blank for most genes and only labels the points we are interested in. The ifelse function in R is a convenient way to set the entries in a vector based on a logical expression. In this case, make the values in Label the same as the gene symbol if the gene is in our list of “top genes”. Otherwise, points get labeled with a blank string "".

For clarity, we also make the points slightly transparent and use a different colour for the text.

N <- 10
top_genes <- results_annotated$ENSEMBL[1:N]
results_annotated %>% 
  mutate(Label = ifelse(ENSEMBL %in% top_genes, GeneID, "")) %>%  
  ggplot(aes(x = log2FoldChange, y = -log10(padj), label=Label)) + geom_point(alpha=0.4) + geom_text(col="blue")

Finally, a slightly better positioning of text is given by the ggrepel package.

if(!require(ggrepel)) install.packages("ggrepel")
N <- 10
top_genes <- results_annotated$ENSEMBL[1:N]
results_annotated %>% 
  mutate(Label = ifelse(ENSEMBL %in% top_genes, GeneID, "")) %>%  
  ggplot(aes(x = log2FoldChange, y = -log10(padj), label=Label)) + geom_point(alpha=0.4) + geom_text(col="blue")

Annotation with the biomaRt resource

The Bioconductor package have the convenience of being able to make queries offline. However, they are only available for certain organisms. If your organism does not have an org.XX.eg.db package listed on the Bioconductor annotation page (http://bioconductor.org/packages/release/BiocViews.html#___AnnotationData), an alternative is to use biomaRt which provides an interface to the popular biomart annotation resource.

The first step is to find the name of a database that you want to connect to

library(biomaRt)
listMarts()

               biomart               version
1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 99
2   ENSEMBL_MART_MOUSE      Mouse strains 99
3     ENSEMBL_MART_SNP  Ensembl Variation 99
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 99

ensembl=useMart("ENSEMBL_MART_ENSEMBL")
# list the available datasets (species). Replace human with the name of your organism
listDatasets(ensembl) %>% filter(grepl("Human",description))

                dataset              description
1 hsapiens_gene_ensembl Human genes (GRCh38.p13)
     version
1 GRCh38.p13

ensembl = useDataset("hsapiens_gene_ensembl", mart=ensembl)

Queries to biomaRt are constructed in a similar way to the queries we performed with the org.Mm.eg.db package. Instead of keys we have filters, and instead of columns we have attributes. The list of acceptable values is much more comprehensive that for the org.Mm.eg.db package.

listFilters(ensembl) %>% 
    filter(grepl("symbol",name))

               name
1       hgnc_symbol
2 uniprot_gn_symbol
                                 description
1                 HGNC symbol(s) [e.g. A1BG]
2 UniProtKB Gene Name symbol(s) [e.g. 5HT1A]

listAttributes(ensembl) %>% 
    filter(grepl("description",name))

                         name
1                 description
2       phenotype_description
3      goslim_goa_description
4        mim_gene_description
5      mim_morbid_description
6      entrezgene_description
7        wikigene_description
8          family_description
9  interpro_short_description
10       interpro_description
11                description
12                description
13                description
14         source_description
15                description
16 somatic_source_description
17                description
                  description         page
1            Gene description feature_page
2       Phenotype description feature_page
3      GOSlim GOA Description feature_page
4        MIM gene description feature_page
5      MIM morbid description feature_page
6       NCBI gene description feature_page
7        WikiGene description feature_page
8  Ensembl Family Description feature_page
9  Interpro Short Description feature_page
10       Interpro Description feature_page
11           Gene description    structure
12           Gene description     homologs
13           Gene description          snp
14 Variant source description          snp
15           Gene description  snp_somatic
16 Variant source description  snp_somatic
17           Gene description    sequences

An advantage over the org.. packages is that positional information can be retrieved

attributeNames <- c('ensembl_gene_id', 'entrezgene_id','description')
getBM(attributes = attributeNames,
      filters = "hgnc_symbol",
      values=top_genes,
      mart=ensembl)

[1] ensembl_gene_id entrezgene_id   description    
<0 rows> (or 0-length row.names)

Challenge 3

1. Use biomaRt to create an data frame containing the entrezgene, Ensembl ID and genomic coordinates (chromosome, start, end) for the Ensembl IDs in the DESeq2 results
2. Remove duplicates entries from the new data frame
3. Join the biomaRt annotation to the DESeq2 results to produce a data frame with differential expression results and annotation
4. Write the joined data frame to a csv file

Interactive graphs and tables

It is often useful to be able to explore our results in an interactive manner; searching for our favourite genes of interest and plotting on-the-fly whether they are statistically significant in our dataset or not.

Such a visualisation is possible with the Glimma Bioconductor package.

It takes our DESeq2 results object, annotation table and normalized counts, and produces a HTML page including a sortable results table, MA-plot and scatter plot. Particular genes can be searched among the table and their expression patterns can be displayed. Alternatively we can click on particular point in the plot and display their stats.

## Repeat the annotation, as the previous annotation table was created using an ordered results table
anno <- AnnotationDbi::select(org.Hs.eg.db,keys=rownames(dds),
              columns=c("SYMBOL","GENENAME","ENSEMBL"),
              keytype="SYMBOL") %>% 
  dplyr::rename(GeneID = SYMBOL) %>% 
  filter(!duplicated(GeneID))

## Make sure we have normalised counts before proceeding
dds <- estimateSizeFactors(dds)

## Load the Glimma package and create the report
library(Glimma)
glMDPlot(de_status,
         anno,
         groups = colData(dds)$Status,
         counts = counts(dds,normalized=TRUE),
         transform = TRUE,
         side.main = "GeneID")