2 Core principles
2.1 High-throughput analysis needs statistics!
Volume of data
Type of research question
- Designed experiments
- Population samples
- …
Technological artifacts
- Differences in sequencing depth between samples
- Bias in the genomic regions sampled
2.2 Scientific research needs to be reproducible
2.2.1 A motivating case study
Cisplatin-resistant non-small-cell lung cancer gene sets
- Baggerly & Coombes 2009 Ann Appl Stat 3: 1309-1334
Lessons
- Record each step of the analysis
- Coordinated manipulation of feature, sample, and assay data
- Informative labels on visualizations
2.2.2 How to be reproducible?
- Use software ‘objects’ that take care of some of the tedious book-keeping
- Document our analysis in scripts and ‘markdown’ documents
2.2.3 Example: SummarizedExperiment
Underlying data is a matrix
- Regions of interest (e.g., genes) x samples
assay()– e.g., matrix of counts of reads overlapping genes
Include information about rows
rowRanges()– gene identifiers, or genomic ranges describing the coordinates of each gene
Include information about columns
colData()– describing samples, experimental design, …
library(airway) # An 'ExperimentData' package...
data(airway) # ...with a sample data set...
airway # ...that is a SummarizedExperiment## class: RangedSummarizedExperiment
## dim: 64102 8
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowData names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSamplehead(assay(airway)) # contains a matrix of counts## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## ENSG00000000003 679 448 873 408 1138 1047 770
## ENSG00000000005 0 0 0 0 0 0 0
## ENSG00000000419 467 515 621 365 587 799 417
## ENSG00000000457 260 211 263 164 245 331 233
## ENSG00000000460 60 55 40 35 78 63 76
## ENSG00000000938 0 0 2 0 1 0 0
## SRR1039521
## ENSG00000000003 572
## ENSG00000000005 0
## ENSG00000000419 508
## ENSG00000000457 229
## ENSG00000000460 60
## ENSG00000000938 0head(rowRanges(airway)) # information about the genes...## GRangesList object of length 6:
## $ENSG00000000003
## GRanges object with 17 ranges and 2 metadata columns:
## seqnames ranges strand | exon_id exon_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] X [99883667, 99884983] - | 667145 ENSE00001459322
## [2] X [99885756, 99885863] - | 667146 ENSE00000868868
## [3] X [99887482, 99887565] - | 667147 ENSE00000401072
## [4] X [99887538, 99887565] - | 667148 ENSE00001849132
## [5] X [99888402, 99888536] - | 667149 ENSE00003554016
## ... ... ... ... . ... ...
## [13] X [99890555, 99890743] - | 667156 ENSE00003512331
## [14] X [99891188, 99891686] - | 667158 ENSE00001886883
## [15] X [99891605, 99891803] - | 667159 ENSE00001855382
## [16] X [99891790, 99892101] - | 667160 ENSE00001863395
## [17] X [99894942, 99894988] - | 667161 ENSE00001828996
##
## ...
## <5 more elements>
## -------
## seqinfo: 722 sequences (1 circular) from an unspecified genomecolData(airway)[, 1:3] # ...and samples## DataFrame with 8 rows and 3 columns
## SampleName cell dex
## <factor> <factor> <factor>
## SRR1039508 GSM1275862 N61311 untrt
## SRR1039509 GSM1275863 N61311 trt
## SRR1039512 GSM1275866 N052611 untrt
## SRR1039513 GSM1275867 N052611 trt
## SRR1039516 GSM1275870 N080611 untrt
## SRR1039517 GSM1275871 N080611 trt
## SRR1039520 GSM1275874 N061011 untrt
## SRR1039521 GSM1275875 N061011 trt## coordinated subsetting
untrt <- airway[, airway$dex == 'untrt']
head(assay(untrt))## SRR1039508 SRR1039512 SRR1039516 SRR1039520
## ENSG00000000003 679 873 1138 770
## ENSG00000000005 0 0 0 0
## ENSG00000000419 467 621 587 417
## ENSG00000000457 260 263 245 233
## ENSG00000000460 60 40 78 76
## ENSG00000000938 0 2 1 0colData(untrt)[, 1:3]## DataFrame with 4 rows and 3 columns
## SampleName cell dex
## <factor> <factor> <factor>
## SRR1039508 GSM1275862 N61311 untrt
## SRR1039512 GSM1275866 N052611 untrt
## SRR1039516 GSM1275870 N080611 untrt
## SRR1039520 GSM1275874 N061011 untrt2.3 We can ‘stand on the shoulders of giants’
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
- ‘Workflows’ make use of multiple packages for complete end-to-end analysis
- ‘Release’ (every six months) and ‘devel’ branches
- Support site; videos, recent courses
2.4 We should explore our data
Visualization
Inter-operability between packages
- Made easier by using similar data structures
Examples (details later)
SummarizedExperimentDNAStringSetGenomicRanges
2.5 Comprehension is more than statistical analysis
Annotation
- Mapping from technical to user-friendly identifiers
- Assigning genes to pathways
- Placing our results in the context of large-scale analyses
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
2.6 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
?readGAlignmentshelp page andvigentte(package="GenomicAlignments", "summarizeOverlaps") - rtracklayer’s
importandexportfunctions 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?importpage 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, …
2.7 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
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
BSgenomepackages. 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.3362762.8 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
GRangesList
2.8.1 Range operations
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 seqlengthsrange(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 seqlengthsreduce(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 seqlengthscoverage(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 1setdiff(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 seqlengthsIRangesList, 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
2.9 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)
require(Rsamtools)
## our 'region of interest'
roi <- GRanges("chr14", IRanges(19653773, width=1))
## sample data
require('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 genome2.10 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 seqlengthsRelated 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
2.11 Integrated data representations: SummarizedExperiment
- ‘feature’ x ‘sample’
assays() colData()data frame for desciption of samplesrowRanges()GRanges / GRangeList or data frame for description of featuresexptData()to describe the entire objectlibrary(SummarizedExperiment) library(airway) data(airway) airway## class: RangedSummarizedExperiment ## dim: 64102 8 ## metadata(1): '' ## assays(1): counts ## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99 ## rowData names(0): ## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521 ## colData names(9): SampleName cell ... Sample BioSamplecolData(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 SAMN02422677airway[, airway$dex %in% "trt"]## class: RangedSummarizedExperiment ## dim: 64102 4 ## metadata(1): '' ## assays(1): counts ## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99 ## rowData names(0): ## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521 ## colData names(9): SampleName cell ... Sample BioSample
2.12 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
?exonsBypage to retrieve all exons grouped by gene or transcript. - BSgenome packages (e.g., BSgenome.Hsapiens.UCSC.hg19) contain whole genomes of model organisms.
- 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
- 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.
2.13 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.,
yieldSizeargument toRsamtools::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()forlapply()-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.