The material in this course requires R version 3.3 and Bioconductor version 3.4
stopifnot(
getRversion() >= '3.3' && getRversion() < '3.4',
BiocInstaller::biocVersion() == "3.4"
)
Version: 0.0.3
Compiled: Wed Jun 22 21:20:26 2016
Physically
Conceptually
Volume of data
Type of research question
Technological artifacts
Cisplatin-resistant non-small-cell lung cancer gene sets
Lessons
SummarizedExperiment
Underlying data is a matrix
assay()
– e.g., matrix of counts of reads overlapping genesInclude information about rows
rowRanges()
– gene identifiers, or genomic ranges describing the coordinates of each geneInclude 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 BioSample
head(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 0
head(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 genome
colData(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 0
colData(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 untrt
Packages!
Visualization
Inter-operability between packages
Examples (details later)
SummarizedExperiment
DNAStringSet
GenomicRanges
Annotation
Objects
methods()
, getClass()
, selectMethod()
method?"substr,<tab>"
to select help on methods, class?D<tab>
for help on classesThis 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
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.
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.
?consensusMatrix
, for instance. Also check out the BSgenome package for working with whole genome sequences, e.g., ?"getSeq,BSgenome-method"
?readGAlignments
help page and vigentte(package="GenomicAlignments", "summarizeOverlaps")
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.Visualization
Classes
Methods –
reverseComplement()
letterFrequency()
matchPDict()
, matchPWM()
Related packages
Example
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
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
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
Classes – GenomicRanges-like behaivor
Methods
readGAlignments()
, readGAlignmentsList()
summarizeOverlaps()
Example
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 genome
Classes – GenomicRanges-like behavior
Functions and methods
readVcf()
, readGeno()
, readInfo()
, readGT()
, writeVcf()
, filterVcf()
locateVariants()
(variants overlapping ranges), predictCoding()
, summarizeVariants()
genotypeToSnpMatrix()
, snpSummary()
Example
## 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
Reference
assays()
colData()
data frame for desciption of samplesrowRanges()
GRanges / GRangeList or data frame for description of featuresexptData()
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
## rowData 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
## rowData names(0):
## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample
?select
?exonsBy
page to retrieve all exons grouped by gene or transcript.open()
, read chunk(s), close()
.yieldSize
argument to Rsamtools::BamFile()
Rsamtools::ScanBamParam()
ShortRead::FastqSampler()
lapply()
-like operationsParallel evaluation in Bioconductor
bplapply()
for lapply()
-like functions, increasingly used by package developers to provide easy, standard way of gaining parallel evaluation.R / Bioconductor
Publications (General Bioconductor)
Other
Acknowledgements
The research reported in this presentation was supported by the National Cancer Institute and the National Human Genome Research Institute of the National Institutes of Health under Award numbers U24CA180996 and U41HG004059, and the National Science Foundation under Award number 1247813. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the National Science Foundation.
sessionInfo()
sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils datasets methods
## [10] base
##
## other attached packages:
## [1] airway_0.107.2 BioC2016Introduction_0.0.3
## [3] Homo.sapiens_1.3.1 GO.db_3.3.0
## [5] OrganismDbi_1.15.1 AnnotationHub_2.5.4
## [7] Gviz_1.17.4 biomaRt_2.29.2
## [9] org.Hs.eg.db_3.3.0 BiocParallel_1.7.4
## [11] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.25.14
## [13] AnnotationDbi_1.35.3 VariantAnnotation_1.19.2
## [15] RNAseqData.HNRNPC.bam.chr14_0.11.0 GenomicAlignments_1.9.4
## [17] Rsamtools_1.25.0 SummarizedExperiment_1.3.5
## [19] Biobase_2.33.0 BSgenome.Hsapiens.UCSC.hg19_1.4.0
## [21] BSgenome_1.41.2 rtracklayer_1.33.7
## [23] GenomicRanges_1.25.8 GenomeInfoDb_1.9.1
## [25] Biostrings_2.41.4 XVector_0.13.2
## [27] IRanges_2.7.11 S4Vectors_0.11.7
## [29] BiocGenerics_0.19.1 ggplot2_2.1.0
## [31] BiocStyle_2.1.10
##
## loaded via a namespace (and not attached):
## [1] httr_1.2.0 splines_3.3.0 Formula_1.2-1
## [4] shiny_0.13.2 interactiveDisplayBase_1.11.3 latticeExtra_0.6-28
## [7] RBGL_1.49.1 yaml_2.1.13 RSQLite_1.0.0
## [10] lattice_0.20-33 biovizBase_1.21.0 chron_2.3-47
## [13] digest_0.6.9 RColorBrewer_1.1-2 colorspace_1.2-6
## [16] htmltools_0.3.5 httpuv_1.3.3 Matrix_1.2-6
## [19] plyr_1.8.4 XML_3.98-1.4 zlibbioc_1.19.0
## [22] xtable_1.8-2 scales_0.4.0 nnet_7.3-12
## [25] survival_2.39-4 magrittr_1.5 mime_0.4
## [28] evaluate_0.9 foreign_0.8-66 graph_1.51.0
## [31] BiocInstaller_1.23.4 tools_3.3.0 data.table_1.9.6
## [34] formatR_1.4 matrixStats_0.50.2 stringr_1.0.0
## [37] munsell_0.4.3 cluster_2.0.4 ensembldb_1.5.8
## [40] RCurl_1.95-4.8 dichromat_2.0-0 bitops_1.0-6
## [43] labeling_0.3 rmarkdown_0.9.6 gtable_0.2.0
## [46] codetools_0.2-14 DBI_0.4-1 reshape2_1.4.1
## [49] R6_2.1.2 gridExtra_2.2.1 knitr_1.13
## [52] Hmisc_3.17-4 stringi_1.1.1 Rcpp_0.12.5
## [55] rpart_4.1-10 acepack_1.3-3.3