Contents

Original Authors: Martin Morgan, Sonali Arora
Presenting Authors: Martin Morgan, Lori Shepherd
Date: 22 July, 2019
Back: Monday labs

Objective: Learn the essentials of Bioconductor data structures

Lessons learned:

1 Classes, methods, and packages

This section focuses on classes, methods, and packages, with the goal being to learn to navigate the help system and interactive discovery facilities.

1.1 Motivation

Sequence analysis is specialized

  • Large data needs to be processed in a memory- and time-efficient manner
  • Specific algorithms have been developed for the unique characteristics of sequence data

Additional considerations

  • Re-use of existing, tested code is easier to do and less error-prone than re-inventing the wheel.
  • Interoperability between packages is easier when the packages share similar data structures.

Solution: use well-defined classes to represent complex data; methods operate on the classes to perform useful functions. Classes and methods are placed together and distributed as packages so that we can all benefit from the hard work and tested code of others.

1.2 Objects

Load the Biostrings and GenomicRanges package

library(Biostrings)
library(GenomicRanges)
  • Bioconductor makes extensive use of classes to represent complicated data types
  • Classes foster interoperability – many different packages can work on the same data – but can be a bit intimidating for the user.
  • Formal ‘S4’ object system
    • Often a class is described on a particular home page, e.g., ?GRanges, and in vignettes, e.g., vignette(package="GenomicRanges"), vignette("GenomicRangesIntroduction")
    • Many methods and classes can be discovered interactively , e.g., methods(class="GRanges") to find out what one can do with a GRanges instance, and methods(findOverlaps) for classes that the findOverlaps() function operates on.
    • In more advanced cases, one can look at the actual definition of a class or method using getClass(), getMethod()
  • Interactive help
    • ?findOverlaps,<tab> to select help on a specific method, ?GRanges-class for help on a class.

1.3 Example & short exercise: Biostrings

Example: Biostrings for DNA sequences

library(Biostrings)                     # Biological sequences
data(phiX174Phage)                      # sample data, see ?phiX174Phage
phiX174Phage
##   A DNAStringSet instance of length 6
##     width seq                                            names               
## [1]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA Genbank
## [2]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA RF70s
## [3]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA SS78
## [4]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA Bull
## [5]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA G97
## [6]  5386 GAGTTTTATCGCTTCCATGACG...GATTGGCGTATCCAACCTGCA 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))      # 'DNAStringSet' methods

Exercises

  1. Load the Biostrings package and phiX174Phage data set. What class is phiX174Phage? Find the help page for the class, and identify interesting functions that apply to it.
  2. Discover vignettes in the Biostrings package with vignette(package="Biostrings"). Add another argument to the vignette function to view the ‘BiostringsQuickOverview’ vignette.
  3. If the internet is available, navigate to the Biostrings landing page on http://bioconductor.org. Do this by visiting the biocViews page. Can you find the BiostringsQuickOverview vignette on the web site?
  4. The following code loads some sample data, 6 versions of the phiX174Phage genome as a DNAStringSet object.

    library(Biostrings)
    data(phiX174Phage)

    Explain what the following code does, and how it works

    m <- consensusMatrix(phiX174Phage)[1:4,]
    polymorphic <- which(colSums(m != 0) > 1)
    mapply(substr, polymorphic, polymorphic, MoreArgs=list(x=phiX174Phage))
    ##         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
    ## Genbank "G"  "G"  "A"  "A"  "C"  "C"  "A"  "G"  "C" 
    ## RF70s   "A"  "A"  "A"  "G"  "C"  "T"  "A"  "G"  "C" 
    ## SS78    "A"  "A"  "A"  "G"  "C"  "T"  "A"  "G"  "C" 
    ## Bull    "G"  "A"  "G"  "A"  "C"  "T"  "A"  "A"  "T" 
    ## G97     "A"  "A"  "G"  "A"  "C"  "T"  "G"  "A"  "C" 
    ## NEB03   "A"  "A"  "A"  "G"  "T"  "T"  "A"  "G"  "C"

2 Working with genomic ranges

2.1 IRanges and GRanges

The IRanges package defines an important class for specifying integer ranges, e.g.,

library(IRanges)
ir <- IRanges(start=c(10, 20, 30), width=5)
ir
## IRanges object with 3 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]        10        14         5
##   [2]        20        24         5
##   [3]        30        34         5

There are many interesting operations to be performed on ranges, e.g, flank() identifies adjacent ranges

flank(ir, 3)
## IRanges object with 3 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         7         9         3
##   [2]        17        19         3
##   [3]        27        29         3

The IRanges class is part of a class hierarchy. To see this, ask R for the class of ir, and for the class definition of the IRanges class

class(ir)
## [1] "IRanges"
## attr(,"package")
## [1] "IRanges"
getClass(class(ir))
## Class "IRanges" [package "IRanges"]
## 
## Slots:
##                                                             
## Name:              start             width             NAMES
## Class:           integer           integer character_OR_NULL
##                                                             
## Name:        elementType   elementMetadata          metadata
## Class:         character DataTable_OR_NULL              list
## 
## Extends: 
## Class "IPosRanges", directly
## Class "IRanges_OR_IPos", directly
## Class "IntegerRanges", by class "IPosRanges", distance 2
## Class "Ranges", by class "IPosRanges", distance 3
## Class "IntegerRanges_OR_missing", by class "IntegerRanges", distance 3
## Class "List", by class "IPosRanges", distance 4
## Class "Vector", by class "IPosRanges", distance 5
## Class "list_OR_List", by class "IPosRanges", distance 5
## Class "Annotated", by class "IPosRanges", distance 6
## Class "vector_OR_Vector", by class "IPosRanges", distance 6
## 
## Known Subclasses: "NormalIRanges", "GroupingIRanges"

Notice that IRanges extends the Ranges class. Now try entering ?flank (?"flank,<tab>" if not using RStudio, where <tab> means to press the tab key to ask for tab completion). You can see that there are help pages for flank operating on several different classes. Select the completion

?"flank,Ranges-method" 

and verify that you’re at the page that describes the method relevant to an IRanges instance. Explore other range-based operations.

The GenomicRanges package extends the notion of ranges to include features relevant to application of ranges in sequence analysis, particularly the ability to associate a range with a sequence name (e.g., chromosome) and a strand. Create a GRanges instance based on our IRanges instance, as follows

library(GenomicRanges)
gr <- GRanges(c("chr1", "chr1", "chr2"), ir, strand=c("+", "-", "+"))
gr
## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     10-14      +
##   [2]     chr1     20-24      -
##   [3]     chr2     30-34      +
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

The notion of flanking sequence has a more nuanced meaning in biology. In particular we might expect that flanking sequence on the + strand would precede the range, but on the minus strand would follow it. Verify that flank applied to a GRanges object has this behavior.

flank(gr, 3)
## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1       7-9      +
##   [2]     chr1     25-27      -
##   [3]     chr2     27-29      +
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Discover what classes GRanges extends, find the help page documenting the behavior of flank when applied to a GRanges object, and verify that the help page documents the behavior we just observed.

class(gr)
## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"
getClass(class(gr))
## Class "GRanges" [package "GenomicRanges"]
## 
## Slots:
##                                                                       
## Name:         seqnames          ranges          strand         seqinfo
## Class:             Rle IRanges_OR_IPos             Rle         Seqinfo
##                                                       
## Name:  elementMetadata     elementType        metadata
## Class:       DataFrame       character            list
## 
## Extends: 
## Class "GenomicRanges", directly
## Class "Ranges", by class "GenomicRanges", distance 2
## Class "GenomicRanges_OR_missing", by class "GenomicRanges", distance 2
## Class "GenomicRanges_OR_GenomicRangesList", by class "GenomicRanges", distance 2
## Class "GenomicRanges_OR_GRangesList", by class "GenomicRanges", distance 2
## Class "List", by class "GenomicRanges", distance 3
## Class "Vector", by class "GenomicRanges", distance 4
## Class "list_OR_List", by class "GenomicRanges", distance 4
## Class "Annotated", by class "GenomicRanges", distance 5
## Class "vector_OR_Vector", by class "GenomicRanges", distance 5
## 
## Known Subclasses: 
## Class "GPos", directly
## Class "UnstitchedGPos", by class "GPos", distance 2
## Class "StitchedGPos", by class "GPos", distance 2
?"flank,GenomicRanges-method"

Notice that the available flank() methods have been augmented by the methods defined in the GenomicRanges package.

It seems like there might be a number of helpful methods available for working with genomic ranges; we can discover some of these from the command line, indicating that the methods should be on the current search() path

showMethods(class="GRanges", where=search())

Use help() to list the help pages in the GenomicRanges package, and vignettes() to view and access available vignettes; these are also available in the RStudio ‘Help’ tab.

help(package="GenomicRanges")
vignette(package="GenomicRanges")
vignette(package="GenomicRanges", "GenomicRangesHOWTOs")

2.2 Range-based operations

Alt Ranges Algebra

Alt Ranges Algebra

Ranges - IRanges - start() / end() / width() - Vector-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()

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(transformed_gr, grl)

3 How to use high-throughput sequence data types in Bioconductor

The following sections briefly summarize some of the most important file types in high-throughput sequence analysis. Briefly review these, or those that are most relevant to your research, before starting on the section Data Representation in R / Bioconductor

Alt Files and the Bioconductor packages that input them

Alt Files and the Bioconductor packages that input them

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

3.2 Reads: FASTQ files

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

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

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

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

3.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";
            ...
            ...

4 Exercises: Data representation in R / Bioconductor

This section briefly illustrates how different high-throughput sequence data types are represented in R / Bioconductor. Select relevant data types for your area of interest, and work through the examples. Take time to consult help pages, understand the output of function calls, and the relationship between standard data formats (summarized in the previous section) and the corresponding R / Bioconductor representation.

4.1 Biostrings (DNA or amino acid sequences)

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 BSgenome packages. The following calculates GC content across chr14.

    library(BSgenome.Hsapiens.UCSC.hg38)
    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.3454924

Exercises

  1. Setup

    library(Biostrings)
    url <- "ftp://ftp.ensembl.org/pub/release-92/fasta/mus_musculus/cds/Mus_musculus.GRCm38.cds.all.fa.gz"
    fl <- BiocFileCache::bfcrpath(rnames = url)
    cds <- rtracklayer::import(fl, "fasta")
  2. For simplicity, clean up the data to remove cds with width not a multiple of three. Remove cds that don’t start with a start codon ATG or end with a stop codon c("TAA", "TAG", "TGA")

    pred1 <- width(cds) %% 3 == 0
    table(pred1)
    ## pred1
    ## FALSE  TRUE 
    ##  7219 58251
    pred2 <- narrow(cds, 1, 3) == "ATG"
    stops <- c("TAA", "TAG", "TGA")
    pred3 <- narrow(cds, width(cds) - 2, width(cds)) %in% stops
    table(pred1 & pred2 & pred3)
    ## 
    ## FALSE  TRUE 
    ## 16808 48662
    cds <- cds[ pred1 & pred2 & pred3 ]
  3. What does the distribution of widths of the cds look like? Which cds has maximum width?

    hist(log10(width(cds)))

    cds[ which.max(width(cds)) ]
    ##   A DNAStringSet instance of length 1
    ##      width seq                                           names               
    ## [1] 105642 ATGACTACTCAAGCACCGATG...ATTAATATCCGTTCTATGTAA ENSMUST0000009998...
    names(cds)[ which.max(width(cds)) ]
    ## [1] "ENSMUST00000099981.8 cds chromosome:GRCm38:2:76703980:76982455:-1 gene:ENSMUSG00000051747.14 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:Ttn description:titin [Source:MGI Symbol;Acc:MGI:98864]"
  4. Use letterFrequency() to calculate the GC content of each cds; visualize the distribution of GC content.

    gc <- letterFrequency(cds, "GC", as.prob=TRUE)
    head(gc)
    ##            G|C
    ## [1,] 0.5026624
    ## [2,] 0.5026624
    ## [3,] 0.5026624
    ## [4,] 0.4879075
    ## [5,] 0.5066079
    ## [6,] 0.5714286
    hist(gc)

    plot( log10(width(cds)), gc, pch=".")

  5. Summarize codon usage in each CDS. Which codons are used most frequently over all CDS?

    AMINO_ACID_CODE
    ##     A     R     N     D     C     Q     E     G     H     I     L     K 
    ## "Ala" "Arg" "Asn" "Asp" "Cys" "Gln" "Glu" "Gly" "His" "Ile" "Leu" "Lys" 
    ##     M     F     P     S     T     W     Y     V     U     O     B     J 
    ## "Met" "Phe" "Pro" "Ser" "Thr" "Trp" "Tyr" "Val" "Sec" "Pyl" "Asx" "Xle" 
    ##     Z     X 
    ## "Glx" "Xaa"
    aa <- translate(cds)
    codon_use <- letterFrequency(aa, names(AMINO_ACID_CODE))
    head(codon_use)
    ##       A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V U O B J Z X
    ## [1,] 17 12 10  7 11 12 11 14  5 20 47 12  6 23 17 24 19  3 11 30 0 0 0 0 0 0
    ## [2,] 17 12 10  7 11 12 11 14  5 20 47 12  6 23 17 24 19  3 11 30 0 0 0 0 0 0
    ## [3,] 17 12 10  7 11 12 11 14  5 20 47 12  6 23 17 24 19  3 11 30 0 0 0 0 0 0
    ## [4,] 25 11  8  6  5  5  7 11 10 34 40 11 10 20 14 36 26  2 11 23 0 0 0 0 0 0
    ## [5,] 36 25 22 26 15 38 30 40 20 27 63 32 11 24 40 86 59 11 34 40 0 0 0 0 0 0
    ## [6,]  0  0  0  0 14 10  5  0  0  0  2  8  1  1 29  3  1  0  0  2 0 0 0 0 0 0
  6. (Advanced) – DNAStringSet inherits from Vector and Annotated, which means that each element (sequence) can have additional information, for instance we can associate GC content with each sequence

    mcols(cds) <- DataFrame(
        GC = gc[,"G|C"]
    )
    mcols(cds, use.names = FALSE)
    ## DataFrame with 48662 rows and 1 column
    ##                      GC
    ##               <numeric>
    ## 1     0.502662406815761
    ## 2     0.502662406815761
    ## 3     0.502662406815761
    ## 4     0.487907465825447
    ## 5     0.506607929515419
    ## ...                 ...
    ## 48658 0.627869185456023
    ## 48659 0.625754919150594
    ## 48660 0.659102370146243
    ## 48661 0.649717514124294
    ## 48662 0.596269554753309
    mcols(cds[1:3], use.names = FALSE)
    ## DataFrame with 3 rows and 1 column
    ##                  GC
    ##           <numeric>
    ## 1 0.502662406815761
    ## 2 0.502662406815761
    ## 3 0.502662406815761

4.2 SummarizedExperiment

Motivation: reproducible & interoperable

  • Matrix of feature x sample measurements, assays()
  • Addition description about samples, colData()

    • Covariates, e.g., age, gender
    • Experimental design, e.g., treatment group
  • Additional information about features, rowData()

    • Gene names, width, GC content, …
    • Genomic ranges(!)
    • Derived values, E.g., log-fold change between treatments, P-value, …
  • Information about the experiment as a whole – metadata()

Example 1: Bulk RNA-seq airway data

  • Attach the airway library and data set

    library(airway)
    ## Loading required package: SummarizedExperiment
    ## Loading required package: Biobase
    ## Welcome to Bioconductor
    ## 
    ##     Vignettes contain introductory material; view with
    ##     'browseVignettes()'. To cite Bioconductor, see
    ##     'citation("Biobase")', and for packages 'citation("pkgname")'.
    ## Loading required package: DelayedArray
    ## Loading required package: matrixStats
    ## 
    ## Attaching package: 'matrixStats'
    ## The following objects are masked from 'package:Biobase':
    ## 
    ##     anyMissing, rowMedians
    ## Loading required package: BiocParallel
    ## 
    ## Attaching package: 'DelayedArray'
    ## The following objects are masked from 'package:matrixStats':
    ## 
    ##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
    ## The following objects are masked from 'package:base':
    ## 
    ##     aperm, apply, rowsum
    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
  • Explore the phenotypic data describing samples. Subset to include just the "untrt" samples.

    colData(airway)
    ## DataFrame with 8 rows and 9 columns
    ##            SampleName     cell      dex    albut        Run avgLength
    ##              <factor> <factor> <factor> <factor>   <factor> <integer>
    ## SRR1039508 GSM1275862   N61311    untrt    untrt SRR1039508       126
    ## SRR1039509 GSM1275863   N61311      trt    untrt SRR1039509       126
    ## SRR1039512 GSM1275866  N052611    untrt    untrt SRR1039512       126
    ## SRR1039513 GSM1275867  N052611      trt    untrt SRR1039513        87
    ## SRR1039516 GSM1275870  N080611    untrt    untrt SRR1039516       120
    ## SRR1039517 GSM1275871  N080611      trt    untrt SRR1039517       126
    ## SRR1039520 GSM1275874  N061011    untrt    untrt SRR1039520       101
    ## SRR1039521 GSM1275875  N061011      trt    untrt SRR1039521        98
    ##            Experiment    Sample    BioSample
    ##              <factor>  <factor>     <factor>
    ## SRR1039508  SRX384345 SRS508568 SAMN02422669
    ## SRR1039509  SRX384346 SRS508567 SAMN02422675
    ## SRR1039512  SRX384349 SRS508571 SAMN02422678
    ## SRR1039513  SRX384350 SRS508572 SAMN02422670
    ## SRR1039516  SRX384353 SRS508575 SAMN02422682
    ## SRR1039517  SRX384354 SRS508576 SAMN02422673
    ## SRR1039520  SRX384357 SRS508579 SAMN02422683
    ## SRR1039521  SRX384358 SRS508580 SAMN02422677
    airway[ , airway$dex == "untrt"]
    ## class: RangedSummarizedExperiment 
    ## dim: 64102 4 
    ## metadata(1): ''
    ## assays(1): counts
    ## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
    ## rowData names(0):
    ## colnames(4): SRR1039508 SRR1039512 SRR1039516 SRR1039520
    ## colData names(9): SampleName cell ... Sample BioSample
  • Calculate library size as the column sums of the assays. Reflect on the relationship between library size and cell / dex column variables and consequences for differential expression analysis.

    colSums(assay(airway))
    ## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 
    ##   20637971   18809481   25348649   15163415   24448408   30818215   19126151 
    ## SRR1039521 
    ##   21164133

Example 2 (advanced): single-cell RNA-seq.

  • Hemberg lab scRNA-seq Datasets ‘Manno’ Mouse brain data set (rds)

    url <- "https://scrnaseq-public-datasets.s3.amazonaws.com/scater-objects/manno_mouse.rds"
    fl <- BiocFileCache::bfcrpath(rnames = url)
    sce <- readRDS(fl)

Exercises

  • What’s the frequency of 0 counts in the single cell assay data?
  • What’s the distribution of library sizes in the single cell assay data?
  • Create a random sample of 100 cells and visualize the relationship between samples using dist() and cmdscale().
  • can you identify what column of colData() is responsible for any pattern you see?
  • In exploring the covariates, are the possible problems with confounding?

4.3 GenomicRanges

Example

library(GenomicRanges)
gr <- GRanges(c("chr1:10-14:+", "chr1:20-24:+", "chr1:22-26:+"))
shift(gr, 1)                            # 1-based coordinates!
## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     11-15      +
##   [2]     chr1     21-25      +
##   [3]     chr1     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]     chr1     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]     chr1     10-14      +
##   [2]     chr1     20-26      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
coverage(gr)
## RleList of length 1
## $chr1
## 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]     chr1     15-19      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Exercises

  1. Which of my SNPs overlap genes?

    genes <- GRanges(c("chr1:30-40:+", "chr1:60-70:-"))
    snps <- GRanges(c("chr1:35", "chr1:60", "chr1:45"))
    countOverlaps(snps, genes) > 0
    ## [1]  TRUE  TRUE FALSE
  2. Which gene is ‘nearest’ my regulatory region? Which gene does my regulatory region precede (i.e., upstream of)

    reg <- GRanges(c("chr1:50-55", "chr1:75-80"))
    nearest(reg, genes)
    ## [1] 2 2
    precede(reg, genes)
    ## [1] NA  2
  3. What range do short reads cover? depth of coverage?

    reads <- GRanges(c("chr1:10-19", "chr1:15-24", "chr1:30-41"))
    coverage(reads, width = 100)
    ## RleList of length 1
    ## $chr1
    ## integer-Rle of length 100 with 7 runs
    ##   Lengths:  9  5  5  5  5 12 59
    ##   Values :  0  1  2  1  0  1  0
    as(coverage(reads, width = 100), "GRanges")
    ## GRanges object with 7 ranges and 1 metadata column:
    ##       seqnames    ranges strand |     score
    ##          <Rle> <IRanges>  <Rle> | <integer>
    ##   [1]     chr1       1-9      * |         0
    ##   [2]     chr1     10-14      * |         1
    ##   [3]     chr1     15-19      * |         2
    ##   [4]     chr1     20-24      * |         1
    ##   [5]     chr1     25-29      * |         0
    ##   [6]     chr1     30-41      * |         1
    ##   [7]     chr1    42-100      * |         0
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome

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.4 GenomicAlignments (Aligned reads)

Classes – GenomicRanges-like behaivor

  • GAlignments, GAlignmentPairs, GAlignmentsList

Methods

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

Exercises

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

    library(GenomicRanges)
    library(GenomicAlignments)
    ## Loading required package: Rsamtools
    library(Rsamtools)
    
    ## our 'region of interest'
    roi <- GRanges("chr14", IRanges(19653773, width=1)) 
    ## sample data
    library('RNAseqData.HNRNPC.bam.chr14')
    bf <- BamFile(RNAseqData.HNRNPC.bam.chr14_BAMFILES[[1]], asMates=TRUE)
    ## alignments, junctions, overlapping our roi
    paln <- readGAlignmentsList(bf)
    j <- summarizeJunctions(paln, with.revmap=TRUE)
    j_overlap <- j[j %over% roi]
    
    ## supporting reads
    paln[j_overlap$revmap[[1]]]
    ## GAlignmentsList object of length 8:
    ## [[1]]
    ## GAlignments object with 2 alignments and 0 metadata columns:
    ##       seqnames strand       cigar    qwidth     start       end     width
    ##          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
    ##   [1]    chr14      -   66M120N6M        72  19653707  19653898       192
    ##   [2]    chr14      +  7M1270N65M        72  19652348  19653689      1342
    ##           njunc
    ##       <integer>
    ##   [1]         1
    ##   [2]         1
    ##   -------
    ##   seqinfo: 93 sequences from an unspecified genome
    ## 
    ## [[2]]
    ## GAlignments object with 2 alignments and 0 metadata columns:
    ##       seqnames strand       cigar    qwidth     start       end     width
    ##          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
    ##   [1]    chr14      -   66M120N6M        72  19653707  19653898       192
    ##   [2]    chr14      +         72M        72  19653686  19653757        72
    ##           njunc
    ##       <integer>
    ##   [1]         1
    ##   [2]         0
    ##   -------
    ##   seqinfo: 93 sequences from an unspecified genome
    ## 
    ## [[3]]
    ## GAlignments object with 2 alignments and 0 metadata columns:
    ##       seqnames strand       cigar    qwidth     start       end     width
    ##          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
    ##   [1]    chr14      +         72M        72  19653675  19653746        72
    ##   [2]    chr14      -   65M120N7M        72  19653708  19653899       192
    ##           njunc
    ##       <integer>
    ##   [1]         0
    ##   [2]         1
    ##   -------
    ##   seqinfo: 93 sequences from an unspecified genome
    ## 
    ## ...
    ## <5 more elements>

4.5 VariantAnnotation (Called variants)

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

Exerises

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

    ## input variants
    library(VariantAnnotation)
    fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
    vcf <- readVcf(fl, "hg19")
    seqlevels(vcf) <- "chr22"
    ## known gene model
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    coding <- locateVariants(rowRanges(vcf),
        TxDb.Hsapiens.UCSC.hg19.knownGene,
        CodingVariants())
    ## Warning in valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE): GRanges object contains 2405 out-of-bound ranges located on
    ##   sequences 75253, 74357, 74359, 74360, 74361, 74362, 74363, 74358,
    ##   74364, 74365, 75254, 75259, 74368, 74369, 74366, 74367, 74370,
    ##   74372, 74373, 74374, 74375, 74378, 74377, 74380, 74381, 75262,
    ##   75263, 75265, 75266, 75268, 75269, 75271, 75273, 75276, 75281,
    ##   75282, 75283, 74389, 74383, 74384, 74385, 74386, 74387, 75287,
    ##   75288, 75286, 75289, 74390, 74391, 74392, 74393, 74394, 75291,
    ##   74395, 74396, 74397, 74398, 75302, 75304, 75305, and 75306. Note
    ##   that ranges located on a sequence whose length is unknown (NA) or on
    ##   a circular sequence are not considered out-of-bound (use
    ##   seqlengths() and isCircular() to get the lengths and circularity
    ##   flags of the underlying sequences). You can use trim() to trim these
    ##   ranges. See ?`trim,GenomicRanges-method` for more information.
    head(coding)
    ## GRanges object with 6 ranges and 9 metadata columns:
    ##                   seqnames    ranges strand | LOCATION  LOCSTART    LOCEND
    ##                      <Rle> <IRanges>  <Rle> | <factor> <integer> <integer>
    ##       rs114335781    chr22  50301422      - |   coding       939       939
    ##         rs8135963    chr22  50301476      - |   coding       885       885
    ##   22:50301488_C/T    chr22  50301488      - |   coding       873       873
    ##   22:50301494_G/A    chr22  50301494      - |   coding       867       867
    ##   22:50301584_C/T    chr22  50301584      - |   coding       777       777
    ##       rs114264124    chr22  50302962      - |   coding       698       698
    ##                     QUERYID        TXID         CDSID      GENEID
    ##                   <integer> <character> <IntegerList> <character>
    ##       rs114335781        24       75253        218562       79087
    ##         rs8135963        25       75253        218562       79087
    ##   22:50301488_C/T        26       75253        218562       79087
    ##   22:50301494_G/A        27       75253        218562       79087
    ##   22:50301584_C/T        28       75253        218562       79087
    ##       rs114264124        57       75253        218563       79087
    ##                         PRECEDEID        FOLLOWID
    ##                   <CharacterList> <CharacterList>
    ##       rs114335781            <NA>            <NA>
    ##         rs8135963            <NA>            <NA>
    ##   22:50301488_C/T            <NA>            <NA>
    ##   22:50301494_G/A            <NA>            <NA>
    ##   22:50301584_C/T            <NA>            <NA>
    ##       rs114264124            <NA>            <NA>
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Related packages

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.6 rtracklayer (Genome annotations)

  • Import BED, GTF, WIG, etc
  • Export GRanges to BED, GTF, WIG, …
  • Access UCSC genome browser

5 Extended exercises

5.1 Summarize overlaps

The goal is to count the number of reads overlapping exons grouped into genes. This type of count data is the basic input for RNASeq differential expression analysis, e.g., through DESeq2 and edgeR.

  1. Identify the regions of interest. We use a ‘TxDb’ package with gene models already defined; the genome (hg19) is determined by the genome used for read alignment in the sample BAM files.

    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
    ## only chromosome 14
    seqlevels(exByGn, pruning.mode="coarse") = "chr14"
  2. Identify the sample BAM files.

    library(RNAseqData.HNRNPC.bam.chr14)
    length(RNAseqData.HNRNPC.bam.chr14_BAMFILES)
    ## [1] 8
  3. Summarize overlaps, optionally in parallel

    ## next 2 lines optional; non-Windows
    library(BiocParallel)
    register(MulticoreParam(workers=detectCores()))
    olaps <- summarizeOverlaps(exByGn, RNAseqData.HNRNPC.bam.chr14_BAMFILES)
  4. Explore our handiwork, e.g., library sizes (column sums), relationship between gene length and number of mapped reads, etc.

    olaps
    ## class: RangedSummarizedExperiment 
    ## dim: 779 8 
    ## metadata(0):
    ## assays(1): counts
    ## rownames(779): 10001 100113389 ... 9950 9985
    ## rowData names(0):
    ## colnames(8): ERR127306 ERR127307 ... ERR127304 ERR127305
    ## colData names(0):
    head(assay(olaps))
    ##           ERR127306 ERR127307 ERR127308 ERR127309 ERR127302 ERR127303
    ## 10001           103       139       109       125       152       168
    ## 100113389         0         0         0         0         0         0
    ## 100113391         0         0         0         0         0         0
    ## 100124539         0         0         0         0         0         0
    ## 100126297         0         0         0         0         0         0
    ## 100126308         0         0         0         0         0         0
    ##           ERR127304 ERR127305
    ## 10001           181       150
    ## 100113389         0         0
    ## 100113391         0         0
    ## 100124539         0         0
    ## 100126297         0         0
    ## 100126308         0         0
    colSums(assay(olaps))                # library sizes
    ## ERR127306 ERR127307 ERR127308 ERR127309 ERR127302 ERR127303 ERR127304 
    ##    340646    373268    371639    331518    313800    331135    331606 
    ## ERR127305 
    ##    329647
    plot(sum(width(olaps)), rowMeans(assay(olaps)), log="xy")
    ## Warning in xy.coords(x, y, xlabel, ylabel, log): 252 y values <= 0 omitted
    ## from logarithmic plot