RNA-seq Analysis in R

Annotation and Visualisation of RNA-seq results

library(biomaRt)
library(DESeq2)
library(tidyverse)

package ‘dplyr’ was built under R version 3.5.1

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

load("Robjects/DE.RData")

Overview

• Getting annotation
• Visualising DE results
• Retrieving gene models
• Exporting browser tracks

Adding annotation to the DESeq2 results

We have a list of significantly differentially expressed genes, but the only annotation we can see is the Ensembl Gene ID, which is not very informative.

There are a number of ways to add annotation. One method is to do this using the org.Mm.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. This is the method we will use today.

Select BioMart database and dataset

The first step is to select the Biomart database we are going to access and which data set we are going to use.

# view the available databases
listMarts()

## set up connection to ensembl database
ensembl=useMart("ENSEMBL_MART_ENSEMBL")
# list the available datasets (species)
listDatasets(ensembl) %>% 
    filter(str_detect(description, "Mouse"))

# specify a data set to use
ensembl = useDataset("mmusculus_gene_ensembl", mart=ensembl)

Query the database

Now we need to set up a query. For this we need to specify three things:

1. What type of information we are going to search the dataset on - called filters. In our case this is Ensembl Gene IDs
2. A vector of the values for our filter - the Ensembl Gene IDs from our DE results table
3. What columns (attributes) of the dataset we want returned.

Returning data from Biomart can take time, so it’s always a good idea to test your query on a small list of values first to make sure it is doing what you want. We’ll just use the first 1000 genes for now.

# check the available "filters" - things you can filter for
listFilters(ensembl) %>% 
    filter(str_detect(name, "ensembl"))

# Set the filter type and values
filterType <- "ensembl_gene_id"
filterValues <- rownames(resLvV)[1:1000]
# check the available "attributes" - things you can retreive
listAttributes(ensembl) %>% 
    head(20)

# Set the list of attributes
attributeNames <- c('ensembl_gene_id', 'entrezgene', 'external_gene_name')
# run the query
annot <- getBM(attributes=attributeNames, 
               filters = filterType, 
               values = filterValues, 
               mart = ensembl)

Batch submitting query [============================] 100% eta:  0s
                                                                   

One-to-many relationships

Let’s inspect the annotation.

head(annot)

dim(annot) # why are there more than 1000 rows?

[1] 1001    3

length(unique(annot$ensembl_gene_id)) # why are there less than 1000 Gene ids?

[1] 999

isDup <- duplicated(annot$ensembl_gene_id)
dup <- annot$ensembl_gene_id[isDup]
annot[annot$ensembl_gene_id%in%dup,]

There are a couple of genes that have multiple entries in the retrieved annotation. This is becaues there are multiple Entrez IDs for a single Ensembl gene. These one-to-many relationships come up frequently in genomic databases, it is important to be aware of them and check when necessary.

We will need to do a little work before adding the annotation to out results table. We could decide to discard one or both of the Entrez ID mappings, or we could concatenate the Entrez IDs so that we don’t lose information.

Retrieve full annotation

Challenge 1

That was just 1000 genes. We need annotations for the entire results table. Also, there may be some other interesting columns in BioMart that we wish to retrieve.

1. Search the attributes and add the following to our list of attributes:
   1. The gene description
      
   2. The genomic position - chromosome, start, end, and strand (4 columns)
   3. The gene biotype
      
2. Query BioMart using all of the genes in our results table (resLvV)
   
3. How many Ensembl genes have multipe Entrez IDs associated with them?
   
4. How many Ensembl genes in resLvV don’t have any annotation? Why is this?

Add annotation to the results table

We can now add the annotation to the results table and then save the results using the write_tsv function, which writes the results out to a tab separated file. To save time we have created an annotation table in which we have modified the cumbersome Biomart column names, added median transcript length (we’ll need this in a later session), and dealt with the one-to-many issues for Entrez IDs.

load("Robjects/Ensembl_annotations.RData")
colnames(ensemblAnnot)

 [1] "GeneID"         "Entrez"         "Symbol"         "Description"   
 [5] "Biotype"        "Chr"            "Start"          "End"           
 [9] "Strand"         "medianTxLength"

annotLvV <- as.data.frame(resLvV) %>% 
    rownames_to_column("GeneID") %>% 
    left_join(ensemblAnnot, "GeneID") %>% 
    rename(logFC=log2FoldChange, FDR=padj)

Finally we can output the annotation DE results using write_csv.

write_tsv(annotLvV, "data/VirginVsLactating_Results_Annotated.txt")

Challenge 2

Have a look at gene symbols for most significant genes by adjusted p-value. Do they make biological sense in the context of comparing gene expression in mammary gland tissue between lactating and virgin mice? You may want to do a quick web search of your favourite gene/protein database

Visualisation

DESeq2 provides a functon called lfcShrink that shrinks log-Fold Change (LFC) estimates towards zero using and empirical Bayes procedure. The reason for doing this is that there is high variance in the LFC estimates when counts are low and this results in lowly expressed genes appearing to be show greater differences between groups that highly expressed genes. The lfcShrink method compensates for this and allows better visualisation and ranking of genes. We will use it for our visualisations of the data.

ddsShrink <- lfcShrink(ddsObj, coef="Status_lactate_vs_virgin")
shrinkLvV <- as.data.frame(ddsShrink) %>%
    rownames_to_column("GeneID") %>% 
    left_join(ensemblAnnot, "GeneID") %>% 
    rename(logFC=log2FoldChange, FDR=padj)

P-value histogram

A quick and easy “sanity check” for our DE results is to generate a p-value histogram. What we should see is a high bar in the 0 - 0.05 and then a roughly uniform tail to the right of this. There is a nice explanation of other possible patterns in the histogram and what to when you see them in this post.

hist(shrinkLvV$pvalue)

MA plots

MA plots are a common way to visualize the results of a differential analysis. We met them briefly towards the end of Session 2. This plot shows the log-Fold Change for each gene against its average expression across all samples in the two conditions being contrasted.

DESeq2 has a handy function for plotting this…

…this is fine for a quick look, but it is not easy to make changes to the way it looks or add things such as gene labels. Perhaps we would like to add labels for the top 20 most significantly differentially expressed genes. Let’s use the package ggplot2 instead.

A Brief Introduction to ggplot2

The ggplot2 package has emerged as an attractive alternative to the traditional plots provided by base R. A full overview of all capabilities of the package is available from the cheatsheet.

In brief:-

• shrinkLvV is our data frame containing the variables we wish to plot
• aes creates a mapping between the variables in our data frame to the aesthetic proprties of the plot:
  • the x-axis will be mapped to log2(baseMean)
  • the y-axis will be mapped to the logFC
• geom_point specifies the particular type of plot we want (in this case a bar plot)
• geom_text allows us to add labels to some or all of the points
  • see the cheatsheet for other plot types

The real advantage of ggplot2 is the ability to change the appearance of our plot by mapping other variables to aspects of the plot. For example, we could colour the points based on a the sample group. To this we can add metadata from the sampleinfo table to the data. The colours are automatically chosen by ggplot2, but we can specifiy particular values.

# add a column with the names of only the top 10 genes
cutoff <- sort(shrinkLvV$pvalue)[10]
shrinkLvV <- shrinkLvV %>% 
    mutate(TopGeneLabel=ifelse(pvalue<=cutoff, Symbol, ""))
ggplot(shrinkLvV, aes(x = log2(baseMean), y=logFC)) + 
    geom_point(aes(colour=FDR < 0.05), pch=20, size=0.5) +
    geom_text(aes(label=TopGeneLabel)) +
    labs(x="mean of normalised counts", y="log fold change")

Volcano plot

Another common visualisation is the volcano plot which displays a measure of significance on the y-axis and fold-change on the x-axis. In this case we use the log2 fold change (logFC) on the x-axis, and on the y-axis we’ll use -log10(FDR). This -log10 transformation is commonly used for p-values as it means that more significant genes have a higher scale. We should first remove the genes that we excluded by the independent filtering process of DESeq2

# first remove the filtered genes (FDR=NA) and create a -log10(FDR) column
filtTab <- shrinkLvV %>% 
    filter(!is.na(FDR)) %>% 
    mutate(`-log10(FDR)` = -log10(FDR))
ggplot(filtTab, aes(x = logFC, y=`-log10(FDR)`)) + 
    geom_point(aes(colour=FDR < 0.05), size=2)

We could limit the values at the top of the plot so that we can see the lower portion more clearly.

filtTab <- filtTab %>% 
    mutate(`-log10(FDR)`=pmin(`-log10(FDR)`, 51))
ggplot(filtTab, aes(x = logFC, y=`-log10(FDR)`)) + 
    geom_point(aes(colour=FDR < 0.05, shape = `-log10(FDR)` > 50), size=2)

Strip chart for gene expression

Before following up on the DE genes with further lab work, a recommended sanity check is to have a look at the expression levels of the individual samples for the genes of interest. We can quickly look at grouped expression by using plotCounts function of DESeq2 to retrieve the normalised expression values from the ddsObj object and then plotting with ggplot2.

# Let's look at the most significantly differentially expressed gene
topgene <- filter(shrinkLvV, Symbol=="Wap")
geneID <- topgene$GeneID
plotCounts(ddsObj, gene = geneID, intgroup = c("CellType", "Status"),
           returnData = T) %>% 
    ggplot(aes(x=Status, y=log2(count))) +
      geom_point(aes(fill=Status), pch=21, size=2) +
      facet_wrap(~CellType) +
      expand_limits(y=0)

Interactive StripChart with Glimma

An interactive version of the volcano plot above that includes the raw per sample values in a separate panel is possible via the glXYPlot function in the Glimma package.

save(annotLvV, shrinkLvV, file="results/Annotated_Results_LvV.RData")

library(Glimma)

group <- str_remove_all(sampleinfo$Group, "[aeiou]")

de <- as.integer(filtTab$FDR <= 0.05)

normCounts <- log2(counts(ddsObj))
filtCounts <- normCounts[filtTab$GeneID,]

glXYPlot(
  x = filtTab$logFC,
  y = -log10(filtTab$FDR),
  xlab = "logFC",
  ylab = "FDR",
  main = "Lactating v Virgin",
  counts = filtCounts,
  groups = group,
  status = de,
  anno = filtTab[, c("GeneID", "Symbol", "Description")],
  folder = "volcano"
)

This function creates an html page (./volcano/XY-Plot.html) with a volcano plot on the left and a plot showing the log-CPM per sample for a selected gene on the right. A search bar is available to search for genes of interest.

Heatmap

We’re going to use the package ComplexHeatmap (Z. Gu, Eils, and Schlesner 2016). We’ll also use circlize to generate a colour scale (Z. Gu et al. 2014).

library(ComplexHeatmap)

Loading required package: grid
========================================
ComplexHeatmap version 1.18.1
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://bioconductor.org/packages/ComplexHeatmap/

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.
========================================

library(circlize)

========================================
circlize version 0.4.4
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: http://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization 
  in R. Bioinformatics 2014.
========================================

We can’t plot the entire data set, let’s just select the top 200 by FDR. We’ll also z-transform the counts.

# get the top genes
sigGenes <- as.data.frame(resLvV) %>% 
    rownames_to_column("GeneID") %>% 
    top_n(150, wt=-padj) %>% 
    pull("GeneID")
# filter the data for the top 200 by padj in the LRT test
plotDat <- vst(ddsObj)[sigGenes,] %>% 
    assay()
z.mat <- t(scale(t(plotDat), center=TRUE, scale=TRUE))

# colour palette
myPalette <- c("red3", "ivory", "blue3")
myRamp = colorRamp2(c(-2, 0, 2), myPalette)

Heatmap(z.mat, name = "z-score",
        col = myRamp,            
        show_row_name = FALSE,
        cluster_columns = FALSE)

we can also split the heat map into clusters and add some annotation.

# cluster the data and split the tree
hcDat <- hclust(dist(z.mat))
cutGroups <- cutree(hcDat, h=4)
ha1 = HeatmapAnnotation(df = colData(ddsObj)[,c("CellType", "Status")])
Heatmap(z.mat, name = "z-score",
        col = myRamp,            
        show_row_name = FALSE,
        cluster_columns = FALSE,
        split=cutGroups,
        rect_gp = gpar(col = "darkgrey", lwd=0.5),
        top_annotation = ha1)

Additional Material

Working with Genomic Locations - Transcript database packages

There is a whole suite of annotation packages that can be used to access and for perform advanced queries on information about the genomic location of genes, trancripts and exons. These are listed on the Bioconductor annotation page and have the prefix TxDb. (where “tx” is “transcript”). In addition there are a large number of packages that make use of these annotations for downstream analyses and visualisations.

Unfortunately, these packages do not cover all species and tend only to be available for UCSC genomes. Thankfully, there is a way to build your own database from either a GTF file or from various online resources such as Biomart using the package GenomicFeatures.

Creating a TxDb database from a GTF file

The created database is only loaded into the current R session. You will need to run this command each time - it can be a little slow.

library(GenomicFeatures)
txMm <- makeTxDbFromGFF(file = "counts/Mus_musculus.GRCm38.80.gtf", format = "gtf")

Creating a TxDb database from Biomart

The created database is only loaded into the current R session. You will need to run this command each time.

library(GenomicFeatures)
txMm <- makeTxDbFromBiomart(dataset="mmusculus_gene_ensembl")

Creating a TxDb package from Biomart

This creates an R package that can be installed just like a package that you might download from Bioconductor or CRAN. This can then loaded as normal whenever it is needed, saving you having to build the database each time.

A little extra work is needed at the command line to build the package from the files produced by this method. Feel free to skip this section if you want - the previous two methods are adequate if you can tolderate the short wait each time you create the database.

library(GenomicFeatures)
makeTxDbPackageFromBiomart(version="0.01",
                           destDir = "~",
                           maintainer="Some One <so@someplace.org>",
                           author="Some One <so@someplace.com>",
                           dataset="mmusculus_gene_ensembl")

This creates a new folder in your home directory (~ is a shortcut for home).

We will not go into detail about how to contruct an R package, or the contents of the package directory. This method generates all the files you need. More information on contructing R packages can be found in Hadely Wickham’s “R Packages” book.

The directory created will be something like TxDb.Mmusculus.BioMart.ENSEMBLMARTENSEMBL.GRCm38.p6. This going to be the packaged name and is referenced in various files in the package directory. We recommend changing it to something more manageable such as TxDb.Mmusculus.Ens.GRCm38. We need to change the directory name, the database file name and each reference in the DESCRIPTION and man/package.Rd files.

To this run the following at the command line (in the terminal not R).

OldName=TxDb.Mmusculus.BioMart.ENSEMBLMARTENSEMBL.GRCm38.p6
NewName=TxDb.Mmusculus.Ens.GRCm38

cd ~
# rename the package directory
mv ${OldName} ${NewName}

cd ${NewName}
# rename the database file
mv inst/extdata/${OldName}.sqlite inst/extdata/${NewName}.sqlite

# replace the references in the old directory
sed -i s/${OldName}/${NewName}/ DESCRIPTION 
sed -i s/${OldName}/${NewName}/ man/package.Rd 

cd ~
# Build the package from the directory of files created by the above command
R CMD build TxDb.Mmusculus.Ens.GRCm38
# Install the package from the tarball created
R CMD INSTALL TxDb.Mmusculus.Ens.GRCm38.p6_0.01.tar.gz

library(TxDb.Mmusculus.Ens.GRCm38)
txMm <- TxDb.Mmusculus.Ens.GRCm38

Retrieving information from TxDb packages

Accessing the information in these TxDb databases is similar to the way in which we accessed information using biomaRt except that filters (the information we are filtering on) are now called keys and attributes (things we want to retrieve) are columns.

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.

library(GenomicFeatures)
columns(txMm)

 [1] "CDSCHROM"   "CDSEND"     "CDSID"      "CDSNAME"    "CDSPHASE"  
 [6] "CDSSTART"   "CDSSTRAND"  "EXONCHROM"  "EXONEND"    "EXONID"    
[11] "EXONNAME"   "EXONRANK"   "EXONSTART"  "EXONSTRAND" "GENEID"    
[16] "TXCHROM"    "TXEND"      "TXID"       "TXNAME"     "TXSTART"   
[21] "TXSTRAND"   "TXTYPE"    

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(txMm)

[1] "CDSID"    "CDSNAME"  "EXONID"   "EXONNAME" "GENEID"   "TXID"    
[7] "TXNAME"  

To extract information we use the select function. Let’s get transcript information for our most highly differentially expressed gene.

keyList <- ensemblAnnot$GeneID[ensemblAnnot$Symbol=="Wap"]
select(txMm, 
       keys=keyList,
       keytype = "GENEID",
       columns=c("TXNAME", "TXCHROM", "TXSTART", "TXEND", "TXSTRAND", "TXTYPE")
      )

'select()' returned 1:many mapping between keys and columns

Challenge 3

Use the txMm to retrieve the exon coordinates for the genes: + ENSMUSG00000021604 + ENSMUSG00000022146 + ENSMUSG00000040118

Overview of GenomicRanges

One of the real strengths of the txdb.. databases is the ability to interface with GenomicRanges, which is the object type used throughout Bioconductor to manipulate Genomic Intervals.

These object types permit us to perform common operations on intervals such as overlapping and counting. We can define the chromosome, start and end position of each region (also strand too, but not shown here).

library(GenomicRanges)
simple_range <- GRanges(seqnames = "1", ranges = IRanges(start=1000, end=2000))
simple_range

GRanges object with 1 range and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        1 1000-2000      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

We don’t have to have all our ranges located on the same chromosome

chrs <- c("13", "15", "5")
start <- c(73000000, 6800000, 15000000)
end <- c(74000000, 6900000, 16000000)
my_ranges <- GRanges(seqnames = rep(chrs, 3),
                     ranges = IRanges(start = rep(start, each = 3),
                                      end = rep(end, each = 3))
                     )
my_ranges

GRanges object with 9 ranges and 0 metadata columns:
      seqnames            ranges strand
         <Rle>         <IRanges>  <Rle>
  [1]       13 73000000-74000000      *
  [2]       15 73000000-74000000      *
  [3]        5 73000000-74000000      *
  [4]       13   6800000-6900000      *
  [5]       15   6800000-6900000      *
  [6]        5   6800000-6900000      *
  [7]       13 15000000-16000000      *
  [8]       15 15000000-16000000      *
  [9]        5 15000000-16000000      *
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

There are a number of useful functions for calculating properties of the data (such as coverage or sorting). Not so much for RNA-seq analysis, but GenomicRanges are used throughout Bioconductor for the analysis of NGS data.

For instance, we can quickly identify overlapping regions between two GenomicRanges.

keys <- c("ENSMUSG00000021604", "ENSMUSG00000022146", "ENSMUSG00000040118")
genePos <- select(txMm,
                  keys = keys,
                  keytype = "GENEID",
                  columns = c("EXONCHROM", "EXONSTART", "EXONEND")
                  )

'select()' returned 1:many mapping between keys and columns

geneRanges <- GRanges(genePos$EXONCHROM, 
                      ranges = IRanges(genePos$EXONSTART, genePos$EXONEND), 
                      GENEID = genePos$GENEID)
geneRanges

GRanges object with 96 ranges and 1 metadata column:
       seqnames            ranges strand |             GENEID
          <Rle>         <IRanges>  <Rle> |        <character>
   [1]       13 73260479-73260653      * | ENSMUSG00000021604
   [2]       13 73264848-73264979      * | ENSMUSG00000021604
   [3]       13 73265458-73265709      * | ENSMUSG00000021604
   [4]       13 73266596-73266708      * | ENSMUSG00000021604
   [5]       13 73267504-73267832      * | ENSMUSG00000021604
   ...      ...               ...    ... .                ...
  [92]        5 16327973-16329883      * | ENSMUSG00000040118
  [93]        5 16326151-16326383      * | ENSMUSG00000040118
  [94]        5 16340707-16341059      * | ENSMUSG00000040118
  [95]        5 16361395-16361875      * | ENSMUSG00000040118
  [96]        5 16362265-16362326      * | ENSMUSG00000040118
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

findOverlaps(my_ranges, geneRanges)

Hits object with 40 hits and 0 metadata columns:
       queryHits subjectHits
       <integer>   <integer>
   [1]         1           1
   [2]         1           2
   [3]         1           3
   [4]         1           4
   [5]         1           5
   ...       ...         ...
  [36]         9          36
  [37]         9          75
  [38]         9          84
  [39]         9          85
  [40]         9          87
  -------
  queryLength: 9 / subjectLength: 96

However, we have to pay attention to the naming convention used for each object. seqlevelsStyle can help.

seqlevelsStyle(simple_range)

[1] "NCBI"    "Ensembl" "MSU6"    "AGPvF"  

seqlevelsStyle(my_ranges)

[1] "NCBI"    "Ensembl" "JGI2.F" 

seqlevelsStyle(geneRanges)

[1] "NCBI"    "Ensembl" "JGI2.F" 

Exporting tracks

It is also possible to save the results of a Bioconductor analysis in a browser to enable interactive analysis and integration with other data types, or sharing with collaborators. For instance, we might want a browser track to indicate where our differentially-expressed genes are located. We shall use the bed format to display these locations. We will annotate the ranges with information from our analysis such as the fold-change and significance.

First we create a data frame for just the DE genes.

sigGenes <- filter(shrinkLvV, FDR <= 0.01)
message("Number of significantly DE genes: ", nrow(sigGenes))

Number of significantly DE genes: 4279

head(sigGenes)

Create a genomic ranges object

Several convenience functions exist to retrieve the structure of every gene from a given TxDb object in one list. The output of exonsBy is a list, where each item in the list is the exon co-ordinates of a particular gene, however, we do not need this level of granularity for the bed output, so we will collapse to a single region for each gene.

First we use the range function to obtain a single range for every gene and tranform to a more convenient object with unlist.

exoRanges <- exonsBy(txMm, "gene") %>% 
    range() %>% 
    unlist()
sigRegions <- exoRanges[na.omit(match(sigGenes$GeneID, names(exoRanges)))]
sigRegions

GRanges object with 4271 ranges and 0 metadata columns:
                     seqnames          ranges strand
                        <Rle>       <IRanges>  <Rle>
  ENSMUSG00000025903        1 4807788-4848410      +
  ENSMUSG00000103280        1 4905751-4906861      -
  ENSMUSG00000033793        1 5070018-5162529      +
  ENSMUSG00000051285        1 7088920-7173628      +
  ENSMUSG00000103509        1 7148110-7152137      +
                 ...      ...             ...    ...
  ENSMUSG00000064354       MT       7013-7696      +
  ENSMUSG00000064357       MT       7927-8607      +
  ENSMUSG00000064363       MT     10167-11544      +
  ENSMUSG00000064367       MT     11742-13565      +
  ENSMUSG00000064368       MT     13552-14070      -
  -------
  seqinfo: 139 sequences (1 circular) from an unspecified genome

For visualisation purposes, we are going to restrict the data to genes that are located on chromosomes 1 to 19 and the sex chromosomes. This can be done with the keepSeqLevels function.

seqlevels(sigRegions)

  [1] "CHR_CAST_EI_MMCHR11_CTG4"  "CHR_CAST_EI_MMCHR11_CTG5" 
  [3] "CHR_MG104_PATCH"           "CHR_MG117_PATCH"          
  [5] "CHR_MG132_PATCH"           "CHR_MG153_PATCH"          
  [7] "CHR_MG171_PATCH"           "CHR_MG184_PATCH"          
  [9] "CHR_MG190_MG3751_PATCH"    "CHR_MG191_PATCH"          
 [11] "CHR_MG209_PATCH"           "CHR_MG3172_PATCH"         
 [13] "CHR_MG3231_PATCH"          "CHR_MG3251_PATCH"         
 [15] "CHR_MG3490_PATCH"          "CHR_MG3496_PATCH"         
 [17] "CHR_MG3530_PATCH"          "CHR_MG3561_PATCH"         
 [19] "CHR_MG3562_PATCH"          "CHR_MG3609_PATCH"         
 [21] "CHR_MG3618_PATCH"          "CHR_MG3627_PATCH"         
 [23] "CHR_MG3648_PATCH"          "CHR_MG3656_PATCH"         
 [25] "CHR_MG3683_PATCH"          "CHR_MG3686_PATCH"         
 [27] "CHR_MG3699_PATCH"          "CHR_MG3700_PATCH"         
 [29] "CHR_MG3712_PATCH"          "CHR_MG3714_PATCH"         
 [31] "CHR_MG3829_PATCH"          "CHR_MG3833_MG4220_PATCH"  
 [33] "CHR_MG3835_PATCH"          "CHR_MG3836_PATCH"         
 [35] "CHR_MG3999_PATCH"          "CHR_MG4136_PATCH"         
 [37] "CHR_MG4138_PATCH"          "CHR_MG4151_PATCH"         
 [39] "CHR_MG4162_PATCH"          "CHR_MG4180_PATCH"         
 [41] "CHR_MG4198_PATCH"          "CHR_MG4200_PATCH"         
 [43] "CHR_MG4209_PATCH"          "CHR_MG4211_PATCH"         
 [45] "CHR_MG4212_PATCH"          "CHR_MG4213_PATCH"         
 [47] "CHR_MG4214_PATCH"          "CHR_MG4222_MG3908_PATCH"  
 [49] "CHR_MG4243_PATCH"          "CHR_MG4248_PATCH"         
 [51] "CHR_MG4249_PATCH"          "CHR_MG4254_PATCH"         
 [53] "CHR_MG4255_PATCH"          "CHR_MG4259_PATCH"         
 [55] "CHR_MG4261_PATCH"          "CHR_MG4264_PATCH"         
 [57] "CHR_MG4265_PATCH"          "CHR_MG4266_PATCH"         
 [59] "CHR_MG4281_PATCH"          "CHR_MG4288_PATCH"         
 [61] "CHR_MG4308_PATCH"          "CHR_MG4310_MG4311_PATCH"  
 [63] "CHR_MG51_PATCH"            "CHR_MG65_PATCH"           
 [65] "CHR_MG74_PATCH"            "CHR_MG89_PATCH"           
 [67] "CHR_MMCHR1_CHORI29_IDD5_1" "CHR_PWK_PHJ_MMCHR11_CTG1" 
 [69] "CHR_PWK_PHJ_MMCHR11_CTG2"  "CHR_PWK_PHJ_MMCHR11_CTG3" 
 [71] "CHR_WSB_EIJ_MMCHR11_CTG1"  "CHR_WSB_EIJ_MMCHR11_CTG2" 
 [73] "CHR_WSB_EIJ_MMCHR11_CTG3"  "1"                        
 [75] "2"                         "3"                        
 [77] "4"                         "5"                        
 [79] "6"                         "7"                        
 [81] "8"                         "9"                        
 [83] "10"                        "11"                       
 [85] "12"                        "13"                       
 [87] "14"                        "15"                       
 [89] "16"                        "17"                       
 [91] "18"                        "19"                       
 [93] "X"                         "Y"                        
 [95] "MT"                        "GL456210.1"               
 [97] "GL456211.1"                "GL456212.1"               
 [99] "GL456213.1"                "GL456216.1"               
[101] "GL456219.1"                "GL456221.1"               
[103] "GL456233.1"                "GL456239.1"               
[105] "GL456350.1"                "GL456354.1"               
[107] "GL456359.1"                "GL456360.1"               
[109] "GL456366.1"                "GL456367.1"               
[111] "GL456368.1"                "GL456370.1"               
[113] "GL456372.1"                "GL456378.1"               
[115] "GL456379.1"                "GL456381.1"               
[117] "GL456382.1"                "GL456383.1"               
[119] "GL456385.1"                "GL456387.1"               
[121] "GL456389.1"                "GL456390.1"               
[123] "GL456392.1"                "GL456393.1"               
[125] "GL456394.1"                "GL456396.1"               
[127] "JH584292.1"                "JH584293.1"               
[129] "JH584294.1"                "JH584295.1"               
[131] "JH584296.1"                "JH584297.1"               
[133] "JH584298.1"                "JH584299.1"               
[135] "JH584300.1"                "JH584301.1"               
[137] "JH584302.1"                "JH584303.1"               
[139] "JH584304.1"               

sigRegions <- keepSeqlevels(sigRegions, 
                            value = c(1:19,"X","Y"),
                            pruning.mode="tidy")
seqlevels(sigRegions)

 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
[15] "15" "16" "17" "18" "19" "X"  "Y" 

Add metadata to GRanges object

A useful propery of GenomicRanges is that we can attach metadata to each range using the mcols function. The metadata can be supplied in the form of a data frame.

mcols(sigRegions) <- sigGenes[match(names(sigRegions), sigGenes$GeneID), ]
sigRegions

GRanges object with 4263 ranges and 17 metadata columns:
                     seqnames            ranges strand |             GeneID
                        <Rle>         <IRanges>  <Rle> |        <character>
  ENSMUSG00000025903        1   4807788-4848410      + | ENSMUSG00000025903
  ENSMUSG00000103280        1   4905751-4906861      - | ENSMUSG00000103280
  ENSMUSG00000033793        1   5070018-5162529      + | ENSMUSG00000033793
  ENSMUSG00000051285        1   7088920-7173628      + | ENSMUSG00000051285
  ENSMUSG00000103509        1   7148110-7152137      + | ENSMUSG00000103509
                 ...      ...               ...    ... .                ...
  ENSMUSG00000033478       19 57361009-57389594      + | ENSMUSG00000033478
  ENSMUSG00000040022       19 59902884-59943654      - | ENSMUSG00000040022
  ENSMUSG00000024993       19 60811585-60836227      + | ENSMUSG00000024993
  ENSMUSG00000024997       19 60864051-60874556      - | ENSMUSG00000024997
  ENSMUSG00000074733       19 61053840-61140840      - | ENSMUSG00000074733
                             baseMean             logFC             lfcSE
                            <numeric>         <numeric>         <numeric>
  ENSMUSG00000025903 724.446609497753 0.647871377826558 0.144229304891497
  ENSMUSG00000103280 11.0727087099247 -1.58750612880154  0.43450320307594
  ENSMUSG00000033793 1263.66334600512 0.877213503228418 0.106454855140448
  ENSMUSG00000051285  1483.9749407736  1.29960059994024 0.176033709290525
  ENSMUSG00000103509 25.8677212181917  1.18134725817845 0.299145205180626
                 ...              ...               ...               ...
  ENSMUSG00000033478 604.940764140757 0.485457965666852 0.143244108885764
  ENSMUSG00000040022 420.277348654596   1.0490335717025 0.177594139808139
  ENSMUSG00000024993 273.706275359975  0.57020888233623 0.163411073404219
  ENSMUSG00000024997 1155.47226515427 0.896987748934961 0.184738597369506
  ENSMUSG00000074733 151.521609493818 0.831352686057738 0.159578210641852
                                  stat               pvalue
                             <numeric>            <numeric>
  ENSMUSG00000025903  4.50342169141735 6.68680141269572e-06
  ENSMUSG00000103280 -3.55360486193329 0.000379989673727633
  ENSMUSG00000033793  8.25203716657208 1.55716974899267e-16
  ENSMUSG00000051285  7.35737302702125 1.87564671729837e-13
  ENSMUSG00000103509  3.84551924990212 0.000120297419390124
                 ...               ...                  ...
  ENSMUSG00000033478  3.38170660330153 0.000720370394458763
  ENSMUSG00000040022  5.87606192865726 4.20141204161541e-09
  ENSMUSG00000024993  3.45286302118214 0.000554670584413306
  ENSMUSG00000024997  4.88409037206884 1.03907412917225e-06
  ENSMUSG00000074733  5.14752752399225 2.63942285560448e-07
                                      FDR      Entrez      Symbol
                                <numeric> <character> <character>
  ENSMUSG00000025903 6.82491399552232e-05       18777      Lypla1
  ENSMUSG00000103280  0.00227338281227062        <NA>     Gm37277
  ENSMUSG00000033793 1.34050472735681e-14      108664     Atp6v1h
  ENSMUSG00000051285 9.64437264773011e-12      319263      Pcmtd1
  ENSMUSG00000103509 0.000849854587410736        <NA>     Gm38372
                 ...                  ...         ...         ...
  ENSMUSG00000033478  0.00387098901405385      226252    Fam160b1
  ENSMUSG00000040022 9.31606807595958e-08       74998   Rab11fip2
  ENSMUSG00000024993  0.00311011136706842       67894      Fam45a
  ENSMUSG00000024997 1.31438732097752e-05       11757       Prdx3
  ENSMUSG00000074733 3.88962198511814e-06      414758      Zfp950
                                                                                                                          Description
                                                                                                                          <character>
  ENSMUSG00000025903                                                          lysophospholipase 1 [Source:MGI Symbol;Acc:MGI:1344588]
  ENSMUSG00000103280                                                        predicted gene, 37277 [Source:MGI Symbol;Acc:MGI:5610505]
  ENSMUSG00000033793                              ATPase, H+ transporting, lysosomal V1 subunit H [Source:MGI Symbol;Acc:MGI:1914864]
  ENSMUSG00000051285 protein-L-isoaspartate (D-aspartate) O-methyltransferase domain containing 1 [Source:MGI Symbol;Acc:MGI:2441773]
  ENSMUSG00000103509                                                        predicted gene, 38372 [Source:MGI Symbol;Acc:MGI:5611600]
                 ...                                                                                                              ...
  ENSMUSG00000033478                               family with sequence similarity 160, member B1 [Source:MGI Symbol;Acc:MGI:2147545]
  ENSMUSG00000040022                                 RAB11 family interacting protein 2 (class I) [Source:MGI Symbol;Acc:MGI:1922248]
  ENSMUSG00000024993                                 family with sequence similarity 45, member A [Source:MGI Symbol;Acc:MGI:1915144]
  ENSMUSG00000024997                                                                peroxiredoxin 3 [Source:MGI Symbol;Acc:MGI:88034]
  ENSMUSG00000074733                                                      zinc finger protein 950 [Source:MGI Symbol;Acc:MGI:2652824]
                            Biotype         Chr     Start       End
                        <character> <character> <integer> <integer>
  ENSMUSG00000025903 protein_coding           1   4807788   4848410
  ENSMUSG00000103280            TEC           1   4905751   4906861
  ENSMUSG00000033793 protein_coding           1   5070018   5162529
  ENSMUSG00000051285 protein_coding           1   7088920   7173628
  ENSMUSG00000103509            TEC           1   7148110   7152137
                 ...            ...         ...       ...       ...
  ENSMUSG00000033478 protein_coding          19  57361009  57389594
  ENSMUSG00000040022 protein_coding          19  59902884  59943654
  ENSMUSG00000024993 protein_coding          19  60811585  60836227
  ENSMUSG00000024997 protein_coding          19  60864051  60874556
  ENSMUSG00000074733 protein_coding          19  61053840  61140840
                        Strand medianTxLength TopGeneLabel
                     <integer>      <numeric>  <character>
  ENSMUSG00000025903         1          903.5             
  ENSMUSG00000103280        -1           1111             
  ENSMUSG00000033793         1           1930             
  ENSMUSG00000051285         1           1906             
  ENSMUSG00000103509         1           4028             
                 ...       ...            ...          ...
  ENSMUSG00000033478         1           5630             
  ENSMUSG00000040022        -1           4482             
  ENSMUSG00000024993         1            826             
  ENSMUSG00000024997        -1           1478             
  ENSMUSG00000074733        -1          758.5             
  -------
  seqinfo: 21 sequences from an unspecified genome

Scores and colour on exported tracks

The .bed file format is commonly used to store genomic locations for display in genome browsers (e.g. the UCSC browser or IGV) as tracks. Rather than just representing the genomic locations, the .bed format is also able to colour each range according to some property of the analysis (e.g. direction and magnitude of change) to help highlight particular regions of interest. A score can also be displayed when a particular region is clicked-on.

For the score we can use the \(-log_{10}\) of the adjusted p-value and colour scheme for the regions based on the fold-change

colorRampPalette is a useful function in base R for constructing a palette between two extremes. When choosing colour palettes, make sure they are colour blind friendly. The red / green colour scheme traditionally-applied to microarrays is a bad choice.

We will also truncate the fold-changes to between -5 and 5 to and divide this range into 10 equal bins

rbPal <- colorRampPalette(c("red", "blue"))
logFC <- pmax(sigRegions$logFC, -5)
logFC <- pmin(logFC , 5)
Cols <- rbPal(10)[as.numeric(cut(logFC, breaks = 10))]

The colours and score have to be saved in the GRanges object as score and itemRgb columns respectively, and will be used to construct the browser track. The rtracklayer package can be used to import and export browsers tracks.

Now we can export the signifcant results from the DE analysis as a .bed track using rtracklayer. You can load the resulting file in IGV, if you wish.

mcols(sigRegions)$score <- -log10(sigRegions$FDR)
mcols(sigRegions)$itemRgb <- Cols
sigRegions

GRanges object with 4263 ranges and 19 metadata columns:
                     seqnames            ranges strand |             GeneID
                        <Rle>         <IRanges>  <Rle> |        <character>
  ENSMUSG00000025903        1   4807788-4848410      + | ENSMUSG00000025903
  ENSMUSG00000103280        1   4905751-4906861      - | ENSMUSG00000103280
  ENSMUSG00000033793        1   5070018-5162529      + | ENSMUSG00000033793
  ENSMUSG00000051285        1   7088920-7173628      + | ENSMUSG00000051285
  ENSMUSG00000103509        1   7148110-7152137      + | ENSMUSG00000103509
                 ...      ...               ...    ... .                ...
  ENSMUSG00000033478       19 57361009-57389594      + | ENSMUSG00000033478
  ENSMUSG00000040022       19 59902884-59943654      - | ENSMUSG00000040022
  ENSMUSG00000024993       19 60811585-60836227      + | ENSMUSG00000024993
  ENSMUSG00000024997       19 60864051-60874556      - | ENSMUSG00000024997
  ENSMUSG00000074733       19 61053840-61140840      - | ENSMUSG00000074733
                             baseMean             logFC             lfcSE
                            <numeric>         <numeric>         <numeric>
  ENSMUSG00000025903 724.446609497753 0.647871377826558 0.144229304891497
  ENSMUSG00000103280 11.0727087099247 -1.58750612880154  0.43450320307594
  ENSMUSG00000033793 1263.66334600512 0.877213503228418 0.106454855140448
  ENSMUSG00000051285  1483.9749407736  1.29960059994024 0.176033709290525
  ENSMUSG00000103509 25.8677212181917  1.18134725817845 0.299145205180626
                 ...              ...               ...               ...
  ENSMUSG00000033478 604.940764140757 0.485457965666852 0.143244108885764
  ENSMUSG00000040022 420.277348654596   1.0490335717025 0.177594139808139
  ENSMUSG00000024993 273.706275359975  0.57020888233623 0.163411073404219
  ENSMUSG00000024997 1155.47226515427 0.896987748934961 0.184738597369506
  ENSMUSG00000074733 151.521609493818 0.831352686057738 0.159578210641852
                                  stat               pvalue
                             <numeric>            <numeric>
  ENSMUSG00000025903  4.50342169141735 6.68680141269572e-06
  ENSMUSG00000103280 -3.55360486193329 0.000379989673727633
  ENSMUSG00000033793  8.25203716657208 1.55716974899267e-16
  ENSMUSG00000051285  7.35737302702125 1.87564671729837e-13
  ENSMUSG00000103509  3.84551924990212 0.000120297419390124
                 ...               ...                  ...
  ENSMUSG00000033478  3.38170660330153 0.000720370394458763
  ENSMUSG00000040022  5.87606192865726 4.20141204161541e-09
  ENSMUSG00000024993  3.45286302118214 0.000554670584413306
  ENSMUSG00000024997  4.88409037206884 1.03907412917225e-06
  ENSMUSG00000074733  5.14752752399225 2.63942285560448e-07
                                      FDR      Entrez      Symbol
                                <numeric> <character> <character>
  ENSMUSG00000025903 6.82491399552232e-05       18777      Lypla1
  ENSMUSG00000103280  0.00227338281227062        <NA>     Gm37277
  ENSMUSG00000033793 1.34050472735681e-14      108664     Atp6v1h
  ENSMUSG00000051285 9.64437264773011e-12      319263      Pcmtd1
  ENSMUSG00000103509 0.000849854587410736        <NA>     Gm38372
                 ...                  ...         ...         ...
  ENSMUSG00000033478  0.00387098901405385      226252    Fam160b1
  ENSMUSG00000040022 9.31606807595958e-08       74998   Rab11fip2
  ENSMUSG00000024993  0.00311011136706842       67894      Fam45a
  ENSMUSG00000024997 1.31438732097752e-05       11757       Prdx3
  ENSMUSG00000074733 3.88962198511814e-06      414758      Zfp950
                                                                                                                          Description
                                                                                                                          <character>
  ENSMUSG00000025903                                                          lysophospholipase 1 [Source:MGI Symbol;Acc:MGI:1344588]
  ENSMUSG00000103280                                                        predicted gene, 37277 [Source:MGI Symbol;Acc:MGI:5610505]
  ENSMUSG00000033793                              ATPase, H+ transporting, lysosomal V1 subunit H [Source:MGI Symbol;Acc:MGI:1914864]
  ENSMUSG00000051285 protein-L-isoaspartate (D-aspartate) O-methyltransferase domain containing 1 [Source:MGI Symbol;Acc:MGI:2441773]
  ENSMUSG00000103509                                                        predicted gene, 38372 [Source:MGI Symbol;Acc:MGI:5611600]
                 ...                                                                                                              ...
  ENSMUSG00000033478                               family with sequence similarity 160, member B1 [Source:MGI Symbol;Acc:MGI:2147545]
  ENSMUSG00000040022                                 RAB11 family interacting protein 2 (class I) [Source:MGI Symbol;Acc:MGI:1922248]
  ENSMUSG00000024993                                 family with sequence similarity 45, member A [Source:MGI Symbol;Acc:MGI:1915144]
  ENSMUSG00000024997                                                                peroxiredoxin 3 [Source:MGI Symbol;Acc:MGI:88034]
  ENSMUSG00000074733                                                      zinc finger protein 950 [Source:MGI Symbol;Acc:MGI:2652824]
                            Biotype         Chr     Start       End
                        <character> <character> <integer> <integer>
  ENSMUSG00000025903 protein_coding           1   4807788   4848410
  ENSMUSG00000103280            TEC           1   4905751   4906861
  ENSMUSG00000033793 protein_coding           1   5070018   5162529
  ENSMUSG00000051285 protein_coding           1   7088920   7173628
  ENSMUSG00000103509            TEC           1   7148110   7152137
                 ...            ...         ...       ...       ...
  ENSMUSG00000033478 protein_coding          19  57361009  57389594
  ENSMUSG00000040022 protein_coding          19  59902884  59943654
  ENSMUSG00000024993 protein_coding          19  60811585  60836227
  ENSMUSG00000024997 protein_coding          19  60864051  60874556
  ENSMUSG00000074733 protein_coding          19  61053840  61140840
                        Strand medianTxLength TopGeneLabel            score
                     <integer>      <numeric>  <character>        <numeric>
  ENSMUSG00000025903         1          903.5              4.16590281703525
  ENSMUSG00000103280        -1           1111              2.64332742777773
  ENSMUSG00000033793         1           1930              13.8727316501172
  ENSMUSG00000051285         1           1906              11.0157260173192
  ENSMUSG00000103509         1           4028              3.07065537697694
                 ...       ...            ...          ...              ...
  ENSMUSG00000033478         1           5630              2.41217806121594
  ENSMUSG00000040022        -1           4482              7.03076734657514
  ENSMUSG00000024993         1            826              2.50722405944991
  ENSMUSG00000024997        -1           1478              4.88127663890353
  ENSMUSG00000074733        -1          758.5              5.41009260375256
                         itemRgb
                     <character>
  ENSMUSG00000025903     #71008D
  ENSMUSG00000103280     #AA0055
  ENSMUSG00000033793     #71008D
  ENSMUSG00000051285     #5500AA
  ENSMUSG00000103509     #5500AA
                 ...         ...
  ENSMUSG00000033478     #71008D
  ENSMUSG00000040022     #5500AA
  ENSMUSG00000024993     #71008D
  ENSMUSG00000024997     #71008D
  ENSMUSG00000074733     #71008D
  -------
  seqinfo: 21 sequences from an unspecified genome

library(rtracklayer)
export(sigRegions , con = "results/topHits.bed")

Extracting Reads

As we have been using counts as our starting point, we haven’t investigated the aligned reads from our experiment, and how they are represented. As you may be aware, aligned reads are usually stored in a bam file that can be manipulated with open-source command-line tools such as samtools and picard. Bioconductor provide a low-level interface to data/bam/sam files in the form of the Rsamtools package. The GenomicAlignments package can also be used to retrieve the reads mapping to a particular genomic region in an efficient manner.

library(GenomicAlignments)

In the directory small_bams there should be .bam files for some of the samples in the example study. The workflow to produce these files is described in a supplmentary page for the course. In brief, the raw reads (fastq) were downloaded from the Short Read Archive (SRA) and aligned with hisat2. Each bam file was named according to the file name in SRA, but we have renamed the files according to their name in the study. An index file (.bai) has been generated for each bam file. In order to reduce the size, the bam files used here only contain a subset of the reads that were aligned in the region chr15:101707000-101713000.

list.files("counts/small_bams/")

 [1] "MCL1.DG.15.sm.bam"                "MCL1.DG.15.sm.bam.bai"           
 [3] "MCL1.DH.15.sm.bam"                "MCL1.DH.15.sm.bam.bai"           
 [5] "MCL1.DI.15.sm.bam"                "MCL1.DI.15.sm.bam.bai"           
 [7] "MCL1.DJ.15.sm.bam"                "MCL1.DJ.15.sm.bam.bai"           
 [9] "MCL1.DK.15.sm.bam"                "MCL1.DK.15.sm.bam.bai"           
[11] "MCL1.DL.15.sm.bam"                "MCL1.DL.15.sm.bam.bai"           
[13] "MCL1.LA.15.sm.bam"                "MCL1.LA.15.sm.bam.bai"           
[15] "MCL1.LB.15.sm.bam"                "MCL1.LB.15.sm.bam.bai"           
[17] "MCL1.LC.15.sm.bam"                "MCL1.LC.15.sm.bam.bai"           
[19] "MCL1.LD.15.sm.bam"                "MCL1.LD.15.sm.bam.bai"           
[21] "MCL1.LE.15.sm.bam"                "MCL1.LE.15.sm.bam.bai"           
[23] "MCL1.LF.15.sm.bam"                "MCL1.LF.15.sm.bam.bai"           
[25] "Mus_musculus.GRCm38.80.chr15.gtf"

The readGAlignments function provides a simple interface to interrogate the aligned reads for a particular sample. It can also utilise the index file in order to retrieve only the reads that correspond to a specific region in an efficient manner. The output includes the genomic location of each aligned read and the CIGAR (Compact Idiosyncratic Gapped Alignment Report); where M denotes an match to the genome and I, D correspond to insertions and deletions.

exo <- exonsBy(txMm, "gene") 
generegion <- exo[["ENSMUSG00000022146"]] %>% 
    keepSeqlevels(value = 15, pruning.mode="tidy")
my.reads <- readGAlignments(file="counts/small_bams/MCL1.DG.15.sm.bam",
                       param=ScanBamParam(which=generegion))
my.reads

GAlignments object with 25419 alignments and 0 metadata columns:
          seqnames strand          cigar    qwidth     start       end
             <Rle>  <Rle>    <character> <integer> <integer> <integer>
      [1]       15      + 81M53311N11M8S       100   6799340   6852742
      [2]       15      +           100M       100   6813575   6813674
      [3]       15      +          3S97M       100   6813579   6813675
      [4]       15      +          6S94M       100   6813579   6813672
      [5]       15      +           100M       100   6813580   6813679
      ...      ...    ...            ...       ...       ...       ...
  [25415]       15      -           100M       100   6874937   6875036
  [25416]       15      -           100M       100   6874941   6875040
  [25417]       15      -          99M1S       100   6874945   6875043
  [25418]       15      +           100M       100   6874962   6875061
  [25419]       15      -           100M       100   6874966   6875065
              width     njunc
          <integer> <integer>
      [1]     53403         1
      [2]       100         0
      [3]        97         0
      [4]        94         0
      [5]       100         0
      ...       ...       ...
  [25415]       100         0
  [25416]       100         0
  [25417]        99         0
  [25418]       100         0
  [25419]       100         0
  -------
  seqinfo: 66 sequences from an unspecified genome

It is possible to tweak the function to retrieve other potentially-useful information from the bam file, such as the mapping quality and flag.

my.reads <- readGAlignments(file="counts/small_bams/MCL1.DG.15.sm.bam",
                       param=ScanBamParam(which=generegion,
                                          what=c("seq","mapq","flag")))
my.reads

GAlignments object with 25419 alignments and 3 metadata columns:
          seqnames strand          cigar    qwidth     start       end
             <Rle>  <Rle>    <character> <integer> <integer> <integer>
      [1]       15      + 81M53311N11M8S       100   6799340   6852742
      [2]       15      +           100M       100   6813575   6813674
      [3]       15      +          3S97M       100   6813579   6813675
      [4]       15      +          6S94M       100   6813579   6813672
      [5]       15      +           100M       100   6813580   6813679
      ...      ...    ...            ...       ...       ...       ...
  [25415]       15      -           100M       100   6874937   6875036
  [25416]       15      -           100M       100   6874941   6875040
  [25417]       15      -          99M1S       100   6874945   6875043
  [25418]       15      +           100M       100   6874962   6875061
  [25419]       15      -           100M       100   6874966   6875065
              width     njunc |                     seq      mapq      flag
          <integer> <integer> |          <DNAStringSet> <integer> <integer>
      [1]     53403         1 | GTTTGGAAGT...TCTCCTAAAC        60         0
      [2]       100         0 | GAAATGTTTT...ATCAATGTCA        60         0
      [3]        97         0 | TTTTGTTTTA...TCAATGTCAT        60         0
      [4]        94         0 | TTTTTTTGTT...AAATCAATGT        60         0
      [5]       100         0 | GTTTTAATTT...TGTCATTAAC        60         0
      ...       ...       ... .                     ...       ...       ...
  [25415]       100         0 | TCTCTTTATG...TTCCCACCAG        60        16
  [25416]       100         0 | TTTATGGCTG...CACCAGTCGC        60        16
  [25417]        99         0 | TGGCTGCATG...AGTCGCCAGA        60        16
  [25418]       100         0 | GTCCACAGCC...GCCTGGAGAA        60         0
  [25419]       100         0 | ACAGCCACGT...GGAGAACCGC        60        16
  -------
  seqinfo: 66 sequences from an unspecified genome

The flag can represent useful QC information. e.g.

• Read is unmapped
• Read is paired / unpaired
• Read failed QC
• Read is a PCR duplicate (see later)

The combination of any of these properties is used to derive a numeric value, as illustrated in this useful resource

Particular attributes of the reads can be extracted and visualised

hist(mcols(my.reads)$mapq, main="", xlab="MAPQ")

However, there are more-sophisticated visualisation options for aligned reads and range data. We will use the ggbio package, which first requires some discussion of the ggplot2 plotting package.

Composing plots with ggbio

We will now take a brief look at one of the visualisation packages in Bioconductor that takes advantage of the GenomicRanges and GenomicFeatures object-types. In this section we will show a worked example of how to combine several types of genomic data on the same plot. The documentation for ggbio is very extensive and contains lots of examples.

http://www.tengfei.name/ggbio/docs/

The Gviz package is another Bioconductor package that specialising in genomic visualisations, but we will not explore this package in the course.

The Manhattan plot is a common way of visualising genome-wide results, especially when one is concerned with the results of a GWAS study and identifying strongly-associated hits.

The profile is supposed to resemble the Manhattan skyline with particular skyscrapers towering about the lower level buildings.

This type of plot is implemented as the plotGrandLinear function. We have to supply a value to display on the y-axis using the aes function, which is inherited from ggplot2. The positioning of points on the x-axis is handled automatically by ggbio, using the ranges information to get the genomic coordinates of the ranges of interest.

To stop the plots from being too cluttered we will consider the top 200 genes only.

library(ggbio)

Need specific help about ggbio? try mailing 
 the maintainer or visit http://tengfei.github.com/ggbio/

Attaching package: 'ggbio'

The following objects are masked from 'package:ggplot2':

    geom_bar, geom_rect, geom_segment, ggsave, stat_bin,
    stat_identity, xlim

top200 <- sigRegions[order(sigRegions$FDR)[1:200]]
plotGrandLinear(top200 , aes(y = logFC))

using coord:genome to parse x scale

ggbio has alternated the colours of the chromosomes. However, an appealing feature of ggplot2 is the ability to map properties of your plot to variables present in your data. For example, we could create a variable to distinguish between up- and down-regulated genes. The variables used for aesthetic mapping must be present in the mcols section of your ranges object.

mcols(top200)$UpRegulated <- mcols(top200)$logFC > 0
plotGrandLinear(top200, aes(y = logFC, col = UpRegulated))

using coord:genome to parse x scale

plotGrandLinear is a special function in ggbio with preset options for the manhattan style of plot. More often, users will call the autoplot function and ggbio will choose the most appropriate layout. One such layout is the karyogram.

autoplot(top200, layout="karyogram", aes(color=UpRegulated, fill=UpRegulated))

Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.

ggbio is also able to plot the structure of genes according to a particular model represented by a GenomicFeatures object, such as the object we created earlier with the exon coordinates for each gene in the GRCm38 genome.

autoplot(txMm, which=exo[["ENSMUSG00000022146"]])

Parsing transcripts...
Parsing exons...
Parsing cds...
Parsing utrs...
------exons...
------cdss...
------introns...
------utr...
aggregating...
Done
Constructing graphics...

We can even plot the location of sequencing reads if they have been imported using readGAlignments function (or similar).

myreg <- exo[["ENSMUSG00000022146"]] %>% 
    GenomicRanges::reduce() %>% 
    flank(width = 1000, both = T) %>% 
    keepSeqlevels(value = 15, pruning.mode="tidy")
bam <- readGappedReads(file="counts/small_bams/MCL1.DG.15.sm.bam",
                       param=ScanBamParam(which=myreg), use.names = TRUE)
autoplot(bam, geom = "rect") + 
    xlim(GRanges("15", IRanges(6800000, 6900000)))

extracting information...

Like ggplot2, ggbio plots can be saved as objects that can later be modified, or combined together to form more complicated plots. If saved in this way, the plot will only be displayed on a plotting device when we query the object. This strategy is useful when we want to add a common element (such as an ideogram) to a plot composition and don’t want to repeat the code to generate the plot every time.

geneMod <- autoplot(txMm, which = myreg)  + 
    xlim(GRanges("15", IRanges(6810000, 6880000)))
reads.MCL1.DG <- autoplot(bam, stat = "coverage")  + 
    xlim(GRanges("15", IRanges(6810000, 6880000))) +
    labs(title="MCL1.DG")
tracks(GRCm38=geneMod, MCL1.DG=reads.MCL1.DG)

---

References

Gu, Zuguang, Roland Eils, and Matthias Schlesner. 2016. “Complex Heatmaps Reveal Patterns and Correlations in Multidimensional Genomic Data.” Bioinformatics.

Gu, Zuguang, Lei Gu, Roland Eils, Matthias Schlesner, and Benedikt Brors. 2014. “Circlize Implements and Enhances Circular Visualization in R.” Bioinformatics 30 (19): 2811–2.