Version: 0.1.1
Compiled: Sat Oct 24 10:12:10 2015

Objectives

Overview of sequence analysis work flows
Key R and Bioconductor concepts and resources.
Main assumptions: alignment to known reference

Time	Topic
09:15 - 10:15	Sequencing work flows and file types
10:15	Tea/Coffee break
10:30 - 12:30	Introduction to R and Bioconductor
12:30	Lunch
13:30 -14:00	Scalable computing

1 Sequencing work flows

Experimental design
- Keep it simple, e.g., ‘control’ and ‘treatment’ groups
- Replicate within treatments!
Wet-lab sequence preparation (figure from http://rnaseq.uoregon.edu/)
- Record covariates, including processing day – likely ‘batch effects’
(Illumina) Sequencing (Bentley et al., 2008, doi:10.1038/nature07517
- Primary output: FASTQ files of short reads and their quality scores
Alignment
- Choose to match task, e.g., Rsubread, Bowtie2 good for ChIPseq, some forms of RNAseq; BWA, GMAP better for variant calling
- Primary output: BAM files of aligned reads
Reduction
- e.g., RNASeq ‘count table’ (simple spreadsheets), DNASeq called variants (VCF files), ChIPSeq peaks (BED, WIG files)
Analysis
- Differential expression, peak identification, ...
Comprehension
- Biological context

Alt Sequencing Ecosystem

2 Sequence data representations

2.1 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

2.2 Reads: FASTQ files

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

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

Quality scores: ‘phred-like’, encoded. See wikipedia

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

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

Header

@HD VN:1.0 SO:coordinate
@SQ SN:chr1 LN:249250621
@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
...
@SQ SN:chrY LN:59373566
@PG ID:TopHat VN:2.0.8b CL:/home/hpages/tophat-2.0.8b.Linux_x86_64/tophat --mate-inner-dist 150 --solexa-quals --max-multihits 5 --no-discordant --no-mixed --coverage-search --microexon-search --library-type fr-unstranded --num-threads 2 --output-dir tophat2_out/ERR127306 /home/hpages/bowtie2-2.1.0/indexes/hg19 fastq/ERR127306_1.fastq fastq/ERR127306_2.fastq

Alignments: ID, flag, alignment and mate

ERR127306.7941162 403 chr14 19653689 3 72M = 19652348 -1413 ...
ERR127306.22648137 145 chr14 19653692 1 72M = 19650044 -3720 ...
ERR127306.933914 339 chr14 19653707 1 66M120N6M = 19653686 -213 ...
ERR127306.11052450 83 chr14 19653707 3 66M120N6M = 19652348 -1551 ...
ERR127306.24611331 147 chr14 19653708 1 65M120N7M = 19653675 -225 ...
ERR127306.2698854 419 chr14 19653717 0 56M120N16M = 19653935 290 ...
ERR127306.2698854 163 chr14 19653717 0 56M120N16M = 19653935 2019 ...

Alignments: sequence and quality

... GAATTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCC *'%%%%%#&&%''#'&%%%)&&%%$%%'%%'&*****$))$)'')'%)))&)%%%%$'%%%%&"))'')%))
... TTGATCAGTCTCATCTGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAG '**)****)*'*&*********('&)****&***(**')))())%)))&)))*')&***********)****
... TGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAGCAGCCTCTGGTTTCT '******&%)&)))&")')'')'*((******&)&'')'))$))'')&))$)**&&****************
... TGAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAGCAGCCTCTGGTTTCT ##&&(#')$')'%&&#)%$#$%"%###&!%))'%%''%'))&))#)&%((%())))%)%)))%*********
... GAGAGTAACTTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAGCAGCCTCTGGTTTCTT )&$'$'$%!&&%&&#!'%'))%''&%'&))))''$""'%'%&%'#'%'"!'')#&)))))%$)%)&'"')))
... TTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAGCAGCCTCTGGTTTCTTCATGTGGCT ++++++++++++++++++++++++++++++++++++++*++++++**++++**+**''**+*+*'*)))*)#
... TTTGTACCCATCACTGATTCCTTCTGAGACTGCCTCCACTTCCCCAGCAGCCTCTGGTTTCTTCATGTGGCT ++++++++++++++++++++++++++++++++++++++*++++++**++++**+**''**+*+*'*)))*)#

Alignments: Tags

... AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:72 YT:Z:UU NH:i:2 CC:Z:chr22 CP:i:16189276 HI:i:0
... AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:72 YT:Z:UU NH:i:3 CC:Z:= CP:i:19921600 HI:i:0
... AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:4 MD:Z:72 YT:Z:UU XS:A:+ NH:i:3 CC:Z:= CP:i:19921465 HI:i:0
... AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:4 MD:Z:72 YT:Z:UU XS:A:+ NH:i:2 CC:Z:chr22 CP:i:16189138 HI:i:0
... AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:5 MD:Z:72 YT:Z:UU XS:A:+ NH:i:3 CC:Z:= CP:i:19921464 HI:i:0
... AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:72 NM:i:0 XS:A:+ NH:i:5 CC:Z:= CP:i:19653717 HI:i:0
... AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:72 NM:i:0 XS:A:+ NH:i:5 CC:Z:= CP:i:19921455 HI:i:1

2.4 Called variants: VCF files

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

Header

 ##fileformat=VCFv4.2
 ##fileDate=20090805
 ##source=myImputationProgramV3.1
 ##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
 ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
 ##phasing=partial
 ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
 ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
 ...
 ##FILTER=<ID=q10,Description="Quality below 10">
 ##FILTER=<ID=s50,Description="Less than 50% of samples have data">
 ...
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
 ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">

Location

 #CHROM POS ID REF ALT QUAL FILTER ...
 20 14370 rs6054257 G A 29 PASS ...
 20 17330 . T A 3 q10 ...
 20 1110696 rs6040355 A G,T 67 PASS ...
 20 1230237 . T . 47 PASS ...
 20 1234567 microsat1 GTC G,GTCT 50 PASS ...

Variant INFO

 #CHROM POS ... INFO ...
 20 14370 ... NS=3;DP=14;AF=0.5;DB;H2 ...
 20 17330 ... NS=3;DP=11;AF=0.017 ...
 20 1110696 ... NS=2;DP=10;AF=0.333,0.667;AA=T;DB ...
 20 1230237 ... NS=3;DP=13;AA=T ...
 20 1234567 ... NS=3;DP=9;AA=G ...

Genotype FORMAT and samples

 ... POS ... FORMAT NA00001 NA00002 NA00003
 ... 14370 ... GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
 ... 17330 ... GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3
 ... 1110696 ... GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4
 ... 1230237 ... GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
 ... 1234567 ... GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3

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

Input: rtracklayer import()

BED: range-based annotation (see http://genome.ucsc.edu/FAQ/FAQformat.html for definition of this and related formats)
WIG / bigWig: dense, continuous-valued data
GTF: gene model

Component coordinates

 7 protein_coding gene 27221129 27224842 . - . ...
 ...
 7 protein_coding transcript 27221134 27224835 . - . ...
 7 protein_coding exon 27224055 27224835 . - . ...
 7 protein_coding CDS 27224055 27224763 . - 0 ...
 7 protein_coding start_codon 27224761 27224763 . - 0 ...
 7 protein_coding exon 27221134 27222647 . - . ...
 7 protein_coding CDS 27222418 27222647 . - 2 ...
 7 protein_coding stop_codon 27222415 27222417 . - 0 ...
 7 protein_coding UTR 27224764 27224835 . - . ...
 7 protein_coding UTR 27221134 27222414 . - . ...

Annotations

 gene_id "ENSG00000005073"; gene_name "HOXA11"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
 ...
 ... transcript_id "ENST00000006015"; transcript_name "HOXA11-001"; transcript_source "ensembl_havana"; tag "CCDS"; ccds_id "CCDS5411";
 ... exon_number "1"; exon_id "ENSE00001147062";
 ... exon_number "1"; protein_id "ENSP00000006015";
 ... exon_number "1";
 ... exon_number "2"; exon_id "ENSE00002099557";
 ... exon_number "2"; protein_id "ENSP00000006015";
 ... exon_number "2";
 ...

3 R

Language and environment for statistical computing and graphics

Full-featured programming language
Interactive and interpretted – convenient and forgiving
Coherent, extensive documentation
Statistical, e.g. factor(), NA
Extensible – CRAN, Bioconductor, github, ...

Vector, class, object

Efficient vectorized calculations on ‘atomic’ vectors logical, integer, numeric, complex, character, byte
Atomic vectors are building blocks for more complicated objects
matrix – atomic vector with ‘dim’ attribute
data.frame – list of equal length atomic vectors
Formal classes represent complicated combinations of vectors, e.g., the return value of lm(), below

Function, generic, method

Functions transform inputs to outputs, perhaps with side effects, e.g., rnorm(1000)
Argument matching first by name, then by position
Functions may define (some) arguments to have default values
Generic functions dispatch to specific methods based on class of argument(s), e.g., print().
Methods are functions that implement specific generics, e.g., print.factor; methods are invoked indirectly, via the generic.

Introspection

General properties, e.g., class(), str()
Class-specific properties, e.g., dim()

Help

?print: help on the generic print
?print.data.frame: help on print method for objects of class data.frame.

Example

x <- rnorm(1000) # atomic vectors
y <- x + rnorm(1000, sd=.5)
df <- data.frame(x=x, y=y) # object of class 'data.frame'
plot(y ~ x, df) # generic plot, method plot.formula

fit <- lm(y ~x, df) # object of class 'lm'
methods(class=class(fit)) # introspection

## [1] add1 alias anova case.names coerce confint 
## [7] cooks.distance deviance dfbeta dfbetas drop1 dummy.coef 
## [13] effects extractAIC family formula hatvalues influence 
## [19] initialize kappa labels logLik model.frame model.matrix 
## [25] nobs plot predict print proj qr 
## [31] residuals rstandard rstudent show simulate slotsFromS3 
## [37] summary variable.names vcov 
## see '?methods' for accessing help and source code

4 Bioconductor

4.1 Overview

Analysis and comprehension of high-throughput genomic data

Statistical analysis: large data, technological artifacts, designed experiments; rigorous
Comprehension: biological context, visualization, reproducibility
High-throughput
Sequencing: RNASeq, ChIPSeq, variants, copy number, ...
Microarrays: expression, SNP, ...
Flow cytometry, proteomics, images, ...

Packages, vignettes, work flows

934 packages
Discover and navigate via biocViews
Package ‘landing page’
Title, author / maintainer, short description, citation, installation instructions, ..., download statistics
All user-visible functions have help pages, most with runnable examples
‘Vignettes’ an important feature in Bioconductor – narrative documents illustrating how to use the package, with integrated code
‘Release’ (every six months) and ‘devel’ branches
Support site; videos, recent courses

Objects

Represent complicated data types
Foster interoperability
S4 object system
Introspection: methods(), getClass(), selectMethod()
‘accessors’ and other documented functions / methods for manipulation, rather than direct access to the object structure
Interactive help
method?"substr,<tab>" to select help on methods, class?D<tab> for help on classes

Example

require(Biostrings) # Biological sequences
data(phiX174Phage) # sample data, see ?phiX174Phage
phiX174Phage

## A DNAStringSet instance of length 6
## width seq names 
## [1] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA Genbank
## [2] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA RF70s
## [3] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA SS78
## [4] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA Bull
## [5] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA G97
## [6] 5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA NEB03

m <- consensusMatrix(phiX174Phage)[1:4,] # nucl. x position counts
polymorphic <- which(colSums(m != 0) > 1)
m[, polymorphic]

## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## A 4 5 4 3 0 0 5 2 0
## C 0 0 0 0 5 1 0 0 5
## G 2 1 2 3 0 0 1 4 0
## T 0 0 0 0 1 5 0 0 1

methods(class=class(phiX174Phage))

## [1] ! != 
## [3] [ [[ 
## [5] [[<- [<- 
## [7] %in% < 
## [9] <= == 
## [11] > >= 
## [13] $ $<- 
## [15] aggregate alphabetFrequency 
## [17] anyNA append 
## [19] as.character as.complex 
## [21] as.data.frame as.env 
## [23] as.integer as.list 
## [25] as.logical as.matrix 
## [27] as.numeric as.raw 
## [29] as.vector c 
## [31] chartr coerce 
## [33] compact compare 
## [35] compareStrings complement 
## [37] consensusMatrix consensusString 
## [39] countOverlaps countPattern 
## [41] countPDict dinucleotideFrequencyTest 
## [43] do.call droplevels 
## [45] duplicated elementLengths 
## [47] elementMetadata elementMetadata<- 
## [49] elementType endoapply 
## [51] eval expand 
## [53] extractAt extractROWS 
## [55] Filter Find 
## [57] findOverlaps hasOnlyBaseLetters 
## [59] head high2low 
## [61] ifelse intersect 
## [63] is.na is.unsorted 
## [65] isEmpty isMatchingEndingAt 
## [67] isMatchingStartingAt lapply 
## [69] length lengths 
## [71] letterFrequency Map 
## [73] match matchPattern 
## [75] matchPDict mcols 
## [77] mcols<- mendoapply 
## [79] metadata metadata<- 
## [81] mstack names 
## [83] names<- narrow 
## [85] nchar neditEndingAt 
## [87] neditStartingAt NROW 
## [89] nucleotideFrequencyAt oligonucleotideFrequency 
## [91] order overlapsAny 
## [93] PairwiseAlignments PairwiseAlignmentsSingleSubject
## [95] parallelSlotNames PDict 
## [97] Position PWM 
## [99] rank Reduce 
## [101] relist relistToClass 
## [103] rename rep 
## [105] rep.int replaceAt 
## [107] replaceLetterAt replaceROWS 
## [109] rev revElements 
## [111] reverse reverseComplement 
## [113] ROWNAMES sapply 
## [115] seqinfo seqinfo<- 
## [117] seqlevelsInUse seqtype 
## [119] seqtype<- setdiff 
## [121] setequal shiftApply 
## [123] show showAsCell 
## [125] sort split 
## [127] split<- splitAsList 
## [129] stack stringDist 
## [131] subseq subseq<- 
## [133] subset subsetByOverlaps 
## [135] table tail 
## [137] tapply threebands 
## [139] toString translate 
## [141] trimLRPatterns twoWayAlphabetFrequency 
## [143] union unique 
## [145] uniqueLetters unlist 
## [147] unsplit unstrsplit 
## [149] updateObject values 
## [151] values<- vcountPattern 
## [153] vcountPDict vmatchPattern 
## [155] vwhichPDict which.isMatchingEndingAt 
## [157] which.isMatchingStartingAt whichPDict 
## [159] width window 
## [161] window<- with 
## [163] within xtfrm 
## [165] xvcopy 
## see '?methods' for accessing help and source code

selectMethod(reverseComplement, class(phiX174Phage))

## Method Definition:
## 
## function (x, ...) 
## xvcopy(x, lkup = getDNAComplementLookup(), reverse = TRUE)
## <environment: namespace:Biostrings>
## 
## Signatures:
## x 
## target "DNAStringSet"
## defined "DNAStringSet"

Alt Sequencing Ecosystem

4.2 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

Bioconductor packages are listed on the biocViews page. Each package has ‘biocViews’ (tags from a controlled vocabulary) associated with it; these can be searched to identify appropriately tagged packages, as can the package title and author.
Each package has a ‘landing page’, e.g., for GenomicRanges. Visit this landing page, and note the description, authors, and installation instructions. Packages are often written up in the scientific literature, and if available the corresponding citation is present on the landing page. Also on the landing page are links to the vignettes and reference manual and, at the bottom, an indication of cross-platform availability and download statistics.
A package needs to be installed once, using the instructions on the landing page. Once installed, the package can be loaded into an R session
```
library(GenomicRanges)
```
and the help system queried interactively, as outlined above:
```
 help(package="GenomicRanges")
 vignette(package="GenomicRanges")
 vignette(package="GenomicRanges", "GenomicRangesHOWTOs")
 ?GRanges
```

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

Important packages for analysis of differential expression include edgeR and DESeq2; both have excellent vignettes for exploration. Additional research methods embodied in Bioconductor packages can be discovered by visiting the biocViews web page, searching for the ‘DifferentialExpression’ view term, and narrowing the selection by searching for ‘RNA seq’ and similar.
Popular ChIP-seq packages include csaw an dDiffBind for comparison of peaks across samples, ChIPQC for quality assessment, and ChIPseeker for annotating results (e.g., discovering nearby genes). What other ChIP-seq packages are listed on the biocViews page?
Working with called variants (VCF files) is facilitated by packages such as VariantAnnotation, VariantFiltering, ensemblVEP, and SomaticSignatures; packages for calling variants include, e.g., h5vc and VariantTools.
Several packages identify copy number variants from sequence data, including cn.mops; from the biocViews page, what other copy number packages are available? The CNTools package provides some useful facilities for comparison of segments across samples.
Microbiome and metagenomic analysis is facilitated by packages such as phyloseq and metagenomeSeq.
Metabolomics, chemoinformatics, image analysis, and many other high-throughput analysis domains are also represented in Bioconductor; explore these via biocViews and title searches.

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.

The Biostrings package is used to represent DNA and other sequences, with many convenient sequence-related functions. Check out the functions documented on the help page ?consensusMatrix, for instance. Also check out the BSgenome package for working with whole genome sequences, e.g., ?"getSeq,BSgenome-method"
The GenomicAlignments package is used to input reads aligned to a reference genome. See for instance the ?readGAlignments help page and vigentte(package="GenomicAlignments", "summarizeOverlaps")
rtracklayer’s import and export functions can read in many common file types, e.g., BED, WIG, GTF, ..., in addition to querying and navigating the UCSC genome browser. Check out the ?import page for basic usage.
The ShortRead and Rsamtools packages can be used for lower-level access to FASTQ and BAM files, respectively. Explore the ShortRead vignette and Scalable Genomics labs to see approaches to effectively processing the large files.

Visualization

The Gviz package provides great tools for visualizing local genomic coordinates and associated data.
epivizr drives the epiviz genome browser from within R; rtracklayer provides easy ways to transfer data to and manipulate UCSC browser sessions.
Additionl packages include ggbio, OmicCircos, ...

4.3 DNA or amino acid sequences: Biostrings, ShortRead, BSgenome

Classes

XString, XStringSet, e.g., DNAString (genomes), DNAStringSet (reads)

Methods –

Cheat sheat
Manipulation, e.g., reverseComplement()
Summary, e.g., letterFrequency()
Matching, e.g., matchPDict(), matchPWM()

Related packages

BSgenome
- Whole-genome representations
- Model and custom
ShortRead
- FASTQ files

Example

Whole-genome sequences are distrubuted by ENSEMBL, NCBI, and others as FASTA files; model organism whole genome sequences are packaged into more user-friendly BSgenome packages. The following calculates GC content across chr14.

 require(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

4.4 Ranges: GenomicRanges, IRanges

Ranges represent: - Data, e.g., aligned reads, ChIP peaks, SNPs, CpG islands, ... - Annotations, e.g., gene models, regulatory elements, methylated regions - Ranges are defined by chromosome, start, end, and strand - Often, metadata is associated with each range, e.g., quality of alignment, strength of ChIP peak

Many common biological questions are range-based - What reads overlap genes? - What genes are ChIP peaks nearest? - ...

The GenomicRanges package defines essential classes and methods

GRanges

Alt

GRangesList

Alt

4.4.1 Range operations

Alt Ranges Algebra

Ranges - IRanges - start() / end() / width() - List-like – length(), subset, etc. - ‘metadata’, mcols() - GRanges - ‘seqnames’ (chromosome), ‘strand’ - Seqinfo, including seqlevels and seqlengths

Intra-range methods - Independent of other ranges in the same object - GRanges variants strand-aware - shift(), narrow(), flank(), promoters(), resize(), restrict(), trim() - See ?"intra-range-methods"

Inter-range methods - Depends on other ranges in the same object - range(), reduce(), gaps(), disjoin() - coverage() (!) - see ?"inter-range-methods"

Between-range methods - Functions of two (or more) range objects - findOverlaps(), countOverlaps(), ..., %over%, %within%, %outside%; union(), intersect(), setdiff(), punion(), pintersect(), psetdiff()

Example

require(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 - List: all elements of the same type - Many *List-aware methods, but a common ‘trick’: apply a vectorized function to the unlisted representaion, then re-list

 grl <- GRangesList(...)
 orig_gr <- unlist(grl)
 transformed_gr <- FUN(orig)
 transformed_grl <- relist(, grl)

Reference

Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, et al. (2013) Software for Computing and Annotating Genomic Ranges. PLoS Comput Biol 9(8): e1003118. doi:10.1371/journal.pcbi.1003118

4.5 Aligned reads: GenomicAlignments, Rsamtools

Classes – GenomicRanges-like behaivor

GAlignments, GAlignmentPairs, GAlignmentsList
SummarizedExperiment
Matrix where rows are indexed by genomic ranges, columns by a DataFrame.

Methods

readGAlignments(), readGAlignmentsList()
Easy to restrict input, iterate in chunks
summarizeOverlaps()

Example

Find reads supporting the junction identified above, at position 19653707 + 66M = 19653773 of chromosome 14

require(GenomicRanges)
require(GenomicAlignments)

## Loading required package: GenomicAlignments
## Loading required package: Rsamtools

require(Rsamtools)
## our 'region of interest'
roi <- GRanges("chr14", IRanges(19653773, width=1)) 
## sample data
require('RNAseqData.HNRNPC.bam.chr14')

## Loading required package: 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

4.6 Called variants: VariantAnnotation, VariantFiltering

Classes – GenomicRanges-like behavior

VCF – ‘wide’
VRanges – ‘tall’

Functions and methods

I/O and filtering: readVcf(), readGeno(), readInfo(), readGT(), writeVcf(), filterVcf()
Annotation: locateVariants() (variants overlapping ranges), predictCoding(), summarizeVariants()
SNPs: genotypeToSnpMatrix(), snpSummary()

Example

Read variants from a VCF file, and annotate with respect to a known gene model

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

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

Related packages

ensemblVEP
- Forward variants to Ensembl Variant Effect Predictor
VariantTools, h5vc
- Call variants
VariantFiltering
- Filter variants using criteria such as coding consequence, MAF, ..., inheritance model

Reference

Obenchain, V, Lawrence, M, Carey, V, Gogarten, S, Shannon, P, and Morgan, M. VariantAnnotation: a Bioconductor package for exploration and annotation of genetic variants. Bioinformatics, first published online March 28, 2014 doi:10.1093/bioinformatics/btu168

4.7 Integrated data representations: SummarizedExperiment

SummarizedExperiment

‘feature’ x ‘sample’ assays()
colData() data frame for desciption of samples
rowRanges() GRanges / GRangeList or data frame for description of features

exptData() to describe the entire object

library(SummarizedExperiment)
library(airway)
data(airway)
airway

## class: RangedSummarizedExperiment 
## dim: 64102 8 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

colData(airway)

## DataFrame with 8 rows and 9 columns
## SampleName cell dex albut Run avgLength Experiment Sample
## <factor> <factor> <factor> <factor> <factor> <integer> <factor> <factor>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345 SRS508568
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346 SRS508567
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349 SRS508571
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350 SRS508572
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353 SRS508575
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354 SRS508576
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357 SRS508579
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358 SRS508580
## BioSample
## <factor>
## SRR1039508 SAMN02422669
## SRR1039509 SAMN02422675
## SRR1039512 SAMN02422678
## SRR1039513 SAMN02422670
## SRR1039516 SAMN02422682
## SRR1039517 SAMN02422673
## SRR1039520 SAMN02422683
## SRR1039521 SAMN02422677

airway[, airway$dex %in% "trt"]

## class: RangedSummarizedExperiment 
## dim: 64102 4 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

4.8 Annotation: org, TxDb, AnnotationHub, biomaRt, ...

Bioconductor provides extensive access to ‘annotation’ resources (see the AnnotationData biocViews hierarchy); some interesting examples to explore during this lab include:
biomaRt, PSICQUIC, KEGGREST and other packages for querying on-line resources; each of these have informative vignettes.
AnnotationDbi is a cornerstone of the Annotation Data packages provided by Bioconductor.
- org packages (e.g., org.Hs.eg.db) contain maps between different gene identifiers, e.g., ENTREZ and SYMBOL. The basic interface to these packages is described on the help page ?select
- TxDb packages (e.g., TxDb.Hsapiens.UCSC.hg19.knownGene) contain gene models (exon coordinates, exon / transcript relationships, etc) derived from common sources such as the hg19 knownGene track of the UCSC genome browser. These packages can be queried, e.g., as described on the ?exonsBy page to retrieve all exons grouped by gene or transcript.
- BSgenome packages (e.g., BSgenome.Hsapiens.UCSC.hg19) contain whole genomes of model organisms.
VariantAnnotation and ensemblVEP provide access to sequence annotation facilities, e.g., to identify coding variants; see the Introduction to VariantAnnotation vignette for a brief introduction.
Take a quick look at the annotation work flow on the Bioconductor web site.

4.9 Scalable computing

Efficient R code

Vectorize!
Reuse others’ work Know – DESeq2, GenomicRanges, Biostrings, dplyr, data.table, Rcpp

Iteration

Chunk-wise
open(), read chunk(s), close().
e.g., yieldSize argument to Rsamtools::BamFile()

Restriction

Limit to columns and / or rows of interest
Exploit domain-specific formats, e.g., BAM files and Rsamtools::ScanBamParam()
Use a data base

Sampling

Iterate through large data, retaining a manageable sample, e.g., ShortRead::FastqSampler()

Parallel evaluation

After writing efficient code
Typically, lapply()-like operations
Cores on a single machine (‘easy’); clusters (more tedious); clouds

Parallel evaluation in Bioconductor

BiocParallel – bplapply() for lapply()-like functions, increasingly used by package developers to provide easy, standard way of gaining parallel evaluation.
GenomicFiles – Framework for working on groups of files, ranges, or ranges x files
Bioconductor AMI (Amazon Machine Instance) including pre-configured StarCluster, and docker containers.

5 Resources

R / Bioconductor

Web site – install, learn, use, develop R / Bioconductor packages
Support – seek help and guidance; also StackOverflow for R programming questions
biocViews – discover packages
Package landing pages, e.g., GenomicRanges, including title, description, authors, installation instructions, vignettes (e.g., GenomicRanges ‘How To’), etc.
Course and other help material (e.g., videos, EdX course, community blogs, ...)

Publications (General Bioconductor)

Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, et al. (2013) Software for Computing and Annotating Genomic Ranges. PLoS Comput Biol 9(8): e1003118. doi: 10.1371/journal.pcbi.1003118
Lawrence, M, and Morgan, M. 2014. Scalable Genomics with R and Bioconductor. Statistical Science 2014, Vol. 29, No. 2, 214-226. http://arxiv.org/abs/1409.2864v1

Other

Lawrence, M. 2014. Software for Enabling Genomic Data Analysis. Bioc2014 conference slides.

Introduction to Sequence Analysis, R, and Bioconductor

Martin Morgan (martin.morgan@roswellpark.org)
Roswell Park Cancer Institute, Buffalo, NY
19 October, 2015

Contents

1 Sequencing work flows

2 Sequence data representations

2.1 DNA / amino acid sequences: FASTA files

2.2 Reads: FASTQ files

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

2.4 Called variants: VCF files

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

3 R

4 Bioconductor

4.1 Overview

4.2 A sequence analysis package tour

4.3 DNA or amino acid sequences: Biostrings, ShortRead, BSgenome

4.4 Ranges: GenomicRanges, IRanges

4.4.1 Range operations

4.5 Aligned reads: GenomicAlignments, Rsamtools

4.6 Called variants: VariantAnnotation, VariantFiltering

4.7 Integrated data representations: SummarizedExperiment

4.8 Annotation: org, TxDb, AnnotationHub, biomaRt, ...

4.9 Scalable computing

5 Resources

Introduction to Sequence Analysis, R, and Bioconductor

Martin Morgan (martin.morgan@roswellpark.org) Roswell Park Cancer Institute, Buffalo, NY 19 October, 2015

Contents

1 Sequencing work flows

2 Sequence data representations

2.1 DNA / amino acid sequences: FASTA files

2.2 Reads: FASTQ files

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

2.4 Called variants: VCF files

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

3 R

4 Bioconductor

4.1 Overview

4.2 A sequence analysis package tour

4.3 DNA or amino acid sequences: Biostrings, ShortRead, BSgenome

4.4 Ranges: GenomicRanges, IRanges

4.4.1 Range operations

4.5 Aligned reads: GenomicAlignments, Rsamtools

4.6 Called variants: VariantAnnotation, VariantFiltering

4.7 Integrated data representations: SummarizedExperiment

4.8 Annotation: org, TxDb, AnnotationHub, biomaRt, ...

4.9 Scalable computing

5 Resources

Martin Morgan (martin.morgan@roswellpark.org)
Roswell Park Cancer Institute, Buffalo, NY
19 October, 2015