Contents

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

1 Bioconductor

Physically

Conceptually

2 Core principles

2.1 High-throughput analysis needs statistics!

Volume of data

Type of research question

Technological artifacts

2.2 Scientific research needs to be reproducible

2.2.1 A motivating case study

  • Cisplatin-resistant non-small-cell lung cancer gene sets

  • Hsu et al. 2007 J Clin Oncol 25: 4350-4357 retracted

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

2.3 We can ‘stand on the shoulders of giants’

Packages!

2.4 We should explore our data

Visualization

Inter-operability between packages

Examples (details later)

2.5 Comprehension is more than statistical analysis

Annotation

Objects

2.6 A sequence analysis package tour

Alt Sequencing Ecosystem

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

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

Classes

Methods –

Related packages

Example

  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

2.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

Alt

Alt

2.8.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

2.9 Aligned reads: GenomicAlignments, Rsamtools

Classes – GenomicRanges-like behaivor

Methods

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

2.10 Called variants: VariantAnnotation, VariantFiltering

Classes – GenomicRanges-like behavior

Functions and methods

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

2.11 Integrated data representations: SummarizedExperiment

SummarizedExperiment

2.12 Annotation: org, TxDb, AnnotationHub, biomaRt, …

2.13 Scalable computing

  1. Efficient R code
  1. Iteration
  1. Restriction
  1. Sampling
  1. Parallel evaluation

Parallel evaluation in Bioconductor

3 Resources

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.

3.1 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