Bioconductor for Sequence Analysis

Martin Morgan
October 29, 2014

Sequence analysis work flows

  1. Experimental design
  2. Wet-lab sequence preparation
  3. (Illumina) Sequencing (Bentley et al., 2008, doi:10.1038/nature07517)
  4. Alignment
  5. Reduction
  6. Analysis
  7. Comprehension

Data movement

Alt Sequencing Ecosystem

Sequence data representations

DNA / amino acid sequences: FASTA files

Input & manipulation: Biostrings

>NM_078863_up_2000_chr2L_16764737_f chr2L:16764737-16766736
gttggtggcccaccagtgccaaaatacacaagaagaagaaacagcatctt
gacactaaaatgcaaaaattgctttgcgtcaatgactcaaaacgaaaatg
...
atgggtatcaagttgccccgtataaaaggcaagtttaccggttgcacggt
>NM_001201794_up_2000_chr2L_8382455_f chr2L:8382455-8384454
ttatttatgtaggcgcccgttcccgcagccaaagcactcagaattccggg
cgtgtagcgcaacgaccatctacaaggcaatattttgatcgcttgttagg
...

Whole genomes: 2bit and .fa formats: rtracklayer, Rsamtools; BSgenome

Reads: FASTQ files

Input & manipulation: ShortRead readFastq(), FastqStreamer(), FastqSampler()

@ERR127302.1703 HWI-EAS350_0441:1:1:1460:19184#0/1
CCTGAGTGAAGCTGATCTTGATCTACGAAGAGAGATAGATCTTGATCGTCGAGGAGATGCTGACCTTGACCT
+
HHGHHGHHHHHHHHDGG<GDGGE@GDGGD<?B8??ADAD<BE@EE8EGDGA3CB85*,77@>>CE?=896=:
@ERR127302.1704 HWI-EAS350_0441:1:1:1460:16861#0/1
GCGGTATGCTGGAAGGTGCTCGAATGGAGAGCGCCAGCGCCCCGGCGCTGAGCCGCAGCCTCAGGTCCGCCC
+
DE?DD>ED4>EEE>DE8EEEDE8B?EB<@3;BA79?,881B?@73;1?########################

Aligned reads: BAM files (e.g., ERR127306_chr14.bam)

Input & manipulation: 'low-level' Rsamtools, scanBam(), BamFile(); 'high-level' GenomicAlignments

Called variants: VCF files

Input and manipulation: VariantAnnotation readVcf(), readInfo(), readGeno() selectively with ScanVcfParam().

Genome annotations: BED, WIG, GTF, etc. files

Input: rtracklayer import()

Sequence data in R / Bioconductor

Role for Bioconductor

Alt Sequencing Ecosystem

Sequences

Biostrings classes for DNA or amino acid sequences

Methods

Related packages

Example 

suppressPackageStartupMessages({  
    library(BSgenome.Hsapiens.UCSC.hg19)
})
chr14_range = GRanges("chr14", IRanges(1, seqlengths(Hsapiens)["chr14"]))
chr14_dna <- getSeq(Hsapiens, chr14_range)
letterFrequency(chr14_dna, "GC", as.prob=TRUE)

##           G|C
## [1,] 0.336276

Ranges

Ranges represent:

Many common biological questions are range-based

The GenomicRanges package defines essential classes and methods

Range operations

Alt Ranges Algebra

Ranges

Intra-range methods

Inter-range methods

Between-range methods

Example

suppressPackageStartupMessages({
    library(GenomicRanges)
})
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1)                            # 1-based coordinates!

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [11, 15]      +
##   [2]        A  [21, 25]      +
##   [3]        A  [23, 27]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

range(gr)                               # intra-range

## GRanges object with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [10, 26]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

reduce(gr)                              # inter-range

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [10, 14]      +
##   [2]        A  [20, 26]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

coverage(gr)

## RleList of length 1
## $A
## integer-Rle of length 26 with 6 runs
##   Lengths: 9 5 5 2 3 2
##   Values : 0 1 0 1 2 1

setdiff(range(gr), gr)                  # 'introns'

## GRanges object with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [15, 19]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

IRangesList, GRangesList

Reference

GenomicAlignments (Aligned reads)

Classes – GenomicRanges-like behaivor

Methods

Example

suppressPackageStartupMessages({
    library(GenomicRanges)
    library(GenomicAlignments)
    library(Rsamtools)
})

## our 'region of interest'
roi <- GRanges("chr14", IRanges(19653773, width=1)) 
## sample data
suppressPackageStartupMessages({
    library('RNAseqData.HNRNPC.bam.chr14')
})
bf <- BamFile(RNAseqData.HNRNPC.bam.chr14_BAMFILES[[1]], asMates=TRUE)
## alignments, junctions, overlapping our roi
paln <- readGAlignmentsList(bf)
j <- summarizeJunctions(paln, with.revmap=TRUE)
j_overlap <- j[j %over% roi]

## supporting reads
paln[j_overlap$revmap[[1]]]

## GAlignmentsList object of length 8:
## [[1]] 
## GAlignments object with 2 alignments and 0 metadata columns:
##       seqnames strand      cigar qwidth    start      end width njunc
##   [1]    chr14      -  66M120N6M     72 19653707 19653898   192     1
##   [2]    chr14      + 7M1270N65M     72 19652348 19653689  1342     1
## 
## [[2]] 
## GAlignments object with 2 alignments and 0 metadata columns:
##       seqnames strand     cigar qwidth    start      end width njunc
##   [1]    chr14      - 66M120N6M     72 19653707 19653898   192     1
##   [2]    chr14      +       72M     72 19653686 19653757    72     0
## 
## [[3]] 
## GAlignments object with 2 alignments and 0 metadata columns:
##       seqnames strand     cigar qwidth    start      end width njunc
##   [1]    chr14      +       72M     72 19653675 19653746    72     0
##   [2]    chr14      - 65M120N7M     72 19653708 19653899   192     1
## 
## ...
## <5 more elements>
## -------
## seqinfo: 93 sequences from an unspecified genome

VariantAnnotation (called variants)

Classes – GenomicRanges-like behavior

Functions and methods

Example

## input variants
suppressPackageStartupMessages({
    library(VariantAnnotation)
})
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
seqlevels(vcf) <- "chr22"
## known gene model
suppressPackageStartupMessages({
   library(TxDb.Hsapiens.UCSC.hg19.knownGene)
})
coding <- locateVariants(rowData(vcf),
    TxDb.Hsapiens.UCSC.hg19.knownGene,
    CodingVariants())
head(coding)

## GRanges object with 6 ranges and 9 metadata columns:
##       seqnames               ranges strand | LOCATION  LOCSTART    LOCEND
##          <Rle>            <IRanges>  <Rle> | <factor> <integer> <integer>
##   [1]    chr22 [50301422, 50301422]      - |   coding       939       939
##   [2]    chr22 [50301476, 50301476]      - |   coding       885       885
##   [3]    chr22 [50301488, 50301488]      - |   coding       873       873
##   [4]    chr22 [50301494, 50301494]      - |   coding       867       867
##   [5]    chr22 [50301584, 50301584]      - |   coding       777       777
##   [6]    chr22 [50302962, 50302962]      - |   coding       698       698
##         QUERYID      TXID     CDSID      GENEID       PRECEDEID
##       <integer> <integer> <integer> <character> <CharacterList>
##   [1]        24     75253    218562       79087                
##   [2]        25     75253    218562       79087                
##   [3]        26     75253    218562       79087                
##   [4]        27     75253    218562       79087                
##   [5]        28     75253    218562       79087                
##   [6]        57     75253    218563       79087                
##              FOLLOWID
##       <CharacterList>
##   [1]                
##   [2]                
##   [3]                
##   [4]                
##   [5]                
##   [6]                
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Related packages

Reference

rtracklayer (Genome annotations)

A sequence analysis package tour

This very open-ended topic points to some of the most prominent Bioconductor packages for sequence analysis. Use the opportunity in this lab to explore the package vignettes and help pages highlighted below; many of the material will be covered in greater detail in subsequent labs and lectures.

Basics

Domain-specific analysis – explore the landing pages, vignettes, and reference manuals of two or three of the following packages.

Working with sequences, alignments, common web file formats, and raw data; these packages rely very heavily on the IRanges / GenomicRanges infrastructure that we will encounter later in the course.

Visualization

Lab

Short read quality assessment

fastqc is a Java program commonly used for summarizing quality of fastq files. It has a straight-forward graphical user interface. Here we will use the command-line version.

  1. From within Rstudio, choose 'Tools –> Shell…', or log on to your Amazon machine instance using a Mac / linux terminal or on Windows the PuTTY program.

  2. Run fastqc on sample fastq files, sending the output to the ~/fastqc_report directory.

    fastqc fastq/*fastq --threads 8 --outdir=fastqc_reports

  3. Study the quality report and resulting on-line documentation: In the Files tab, click on fastqc_reports. Click on the HTML file there and then click on “View in Web Browser”.

ShortRead provides similar functionality, but from within R. The following shows that R can handle large data, and illustrates some of the basic ways in which one might interact with functionality provided by a Bioconductor package.

## 1. attach ShortRead and BiocParallel
suppressPackageStartupMessages({
    library(ShortRead)
    library(BiocParallel)
})

## 2. create a vector of file paths
fls <- dir("~/fastq", pattern="*fastq", full=TRUE)

## 3. collect statistics
stats0 <- qa(fls)

## 4. generate and browse the report
if (interactive())
    browseURL(report(stats0))

Check out the qa report from all lanes

data(qa_all)
if (interactive())
    browseURL(report(qa_all))

Alignments (and genomic annotations)

This data is from the airway Bioconductor annotation package; see the vignette for details

Integrative Genomics Viewer

  1. Start IGV and select the “hg19” genome.

  2. The sequence names used in the reference genome differ from those used by IGV to represent the identical genome. We need to map between these different sequence names, following the instructions for Creating a Chromosome Name Alias File.

    Copy the file hg19_alias.tab from the location specified in class into the directory <user_home>/igv/genomes/. Restart IGV.

  3. Start igv.

  4. Choose hg19 from the drop-down menu at the top left of the screen

  5. Use File -> Load from URL menu to load a bam file. The URLs will be provided during class.

  6. Zoom in to a particular gene, e.g., SPARCL1, by entering the gene symbol in the box toward the center of the browser window. Adjust the zoom until reads come in to view, and interpret the result.

Bioconductor: we'll explore how to map between different types of identifiers, how to navigate genomic coordinates, and how to query BAM files for aligned reads.

  1. Attach 'Annotation' packages containing information about gene symbols org.Hs.eg.db and genomic coordinates (e.g., genes, exons, cds, transcripts) r Biocannopkg(TxDb.Hsapiens.UCSC.hg19.knownGene). Arrange for the 'seqlevels' (chromosome names) in the TxDb package to match those in the BAM files.

  2. Use an appropriate org.* package to map from gene symbol to Entrez gene id, and the appropriate TxDb.* package to retrieve gene coordinates of the SPARCL1 gene. N.B. – The following uses a single gene symbol, but we could have used 1, 2, or all gene symbols in a vectorized fashion.

  3. Attach the GenomicAlignments package for working with aligned reads. Use range() to get the genomic coordinates spanning the first and last exon of SPARCL1. Input paired reads overlapping SPARCL1.

  4. What questions can you easily answer about these alignments? E.g., how many reads overlap this region of interest?

    ## 1.a 'Annotation' packages
suppressPackageStartupMessages({
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    library(org.Hs.eg.db)
})

## 1.b -- map 'seqlevels' as recorded in the TxDb file to those in the
## BAM file
fl <- "~/igv/genomes/hg19_alias.tab"
map <- with(read.delim(fl, header=FALSE, stringsAsFactors=FALSE),
    setNames(V1, V2))
seqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, force=TRUE) <- map

## 2. Symbol -> Entrez ID -> Gene coordinates
sym2eg <- select(org.Hs.eg.db, "SPARCL1", "ENTREZID", "SYMBOL")
exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
sparcl1exons <- exByGn[[sym2eg$ENTREZID]]

## 3. Aligned reads
suppressPackageStartupMessages({
    library(GenomicAlignments)
})

fl <- "~/bam/SRR1039508_sorted.bam"
sparcl1gene <- range(sparcl1exons)
param <- ScanBamParam(which=sparcl1gene)
aln <- readGAlignmentPairs(fl, param=param)

  5. As another exercise we ask how many of the reads we've input are compatible with the known gene model. We have to find the transcripts that belong to our gene, and then exons grouped by transcript

    ## 5.a. exons-by-transcript for our gene of interest
txids <- select(TxDb.Hsapiens.UCSC.hg19.knownGene, sym2eg$ENTREZID,
    "TXID", "GENEID")$TXID
exByTx <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "tx")[txids]

## 5.b compatible alignments
hits <- findCompatibleOverlaps(query=aln, subject=exByTx)
good <- seq_along(aln) %in% queryHits(hits)
table(good)

    ## good
## FALSE  TRUE 
##    14    55

  6. Finally, let's go from gene model to protein coding sequence. (a) Extract CDS regions grouped by transcript, select just transcripts we're interested in, (b) attach and then extract the coding sequence from the appropriate reference genome. Translating the coding sequences to proteins.

    ## reset seqlevels
restoreSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene)

    ## TxDb object:
## | Db type: TxDb
## | Supporting package: GenomicFeatures
## | Data source: UCSC
## | Genome: hg19
## | Organism: Homo sapiens
## | UCSC Table: knownGene
## | Resource URL: http://genome.ucsc.edu/
## | Type of Gene ID: Entrez Gene ID
## | Full dataset: yes
## | miRBase build ID: GRCh37
## | transcript_nrow: 82960
## | exon_nrow: 289969
## | cds_nrow: 237533
## | Db created by: GenomicFeatures package from Bioconductor
## | Creation time: 2014-09-26 11:16:12 -0700 (Fri, 26 Sep 2014)
## | GenomicFeatures version at creation time: 1.17.17
## | RSQLite version at creation time: 0.11.4
## | DBSCHEMAVERSION: 1.0

    ## a. cds coordinates, grouped by transcript
txids <- select(TxDb.Hsapiens.UCSC.hg19.knownGene, sym2eg$ENTREZID,
    "TXID", "GENEID")$TXID
cdsByTx <- cdsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "tx")[txids]

## b. coding sequence from relevant reference genome
suppressPackageStartupMessages({
    library(BSgenome.Hsapiens.UCSC.hg19)
})

dna <- extractTranscriptSeqs(BSgenome.Hsapiens.UCSC.hg19, cdsByTx)
protein <- translate(dna)

Working with genomic ranges

Visit the “GenomicRanges HOWTOs” vignette. 

browseVignettes("GenomicRanges")

Read section 1, and do exercises 2.2, 2.4, 2.5, 2.8, 2.12, and 2.13. Perhaps select additional topics of particular interest to you.

Resources

R / Bioconductor

Publications and presentations