%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{2. Sequence Data Representations}

# Sequence Data Representations

Epigenomics 2014<br>
Author: Martin Morgan (<a href="mailto:mtmorgan@fhcrc.org">mtmorgan@fhcrc.org</a>)<br>
Date: 24 August, 2014

```{r setup, echo=FALSE}
knitr::opts_chunk$set(cache=TRUE)
```

## Bioconductor for sequence analysis

![Alt Sequencing Ecosystem](figures/SequencingEcosystem.png)

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


### 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](http://en.wikipedia.org/wiki/FASTQ_format#Encoding)

Example

```{r Biostrings, message=FALSE}
require(Biostrings)                     # Biological sequences
data(phiX174Phage)                      # sample data, see ?phiX174Phage
phiX174Phage
m <- consensusMatrix(phiX174Phage)[1:4,] # nucl. x position counts
polymorphic <- which(colSums(m != 0) > 1)
m[, polymorphic]
```
```{r showMethods, eval=FALSE}
showMethods(class=class(phiX174Phage), where=search())
```

### Case study: Working with DNA sequence data 

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. 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.
   ```{r phiX}
   library(Biostrings)
   data(phiX174Phage)
   ```
   Explain what the following code does, and how it works
   ```{r consensusMatrix}
   m <- consensusMatrix(phiX174Phage)[1:4,]
   polymorphic <- which(colSums(m != 0) > 1)
   mapply(substr, polymorphic, polymorphic, MoreArgs=list(x=phiX174Phage))
   ```

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

### 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
            
### 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";
              ...
              
## Data representation

### Ranges

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

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


The [GenomicRanges][] package defines essential classes and methods

- `GRanges`

![Alt ](figures/GRanges.png)

- `GRangesList`

![Alt ](figures/GRangesList.png)

### Range operations

![Alt Ranges Algebra](figures/RangeOperations.png)

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

```{r ranges, message=FALSE}
require(GenomicRanges)
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1)                            # 1-based coordinates!
range(gr)                               # intra-range
reduce(gr)                              # inter-range
coverage(gr)
setdiff(range(gr), gr)                  # 'introns'
```

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


### [Biostrings][] (DNA or amino acid sequences)

Classes

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

Methods --

- [Cheat sheat](http://bioconductor.org/packages/release/bioc/vignettes/Biostrings/inst/doc/BiostringsQuickOverview.pdf)
- Manipulation, e.g., `reverseComplement()`
- Summary, e.g., `letterFrequency()`
- Matching, e.g., `matchPDict()`, `matchPWM()`

Related packages

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

Example 

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

```{r BSgenome-require, message=FALSE}
  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)
```
  
### [GenomicAlignments][] (Aligned reads)

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

```{r bam-require}
  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]]]
```
  
### [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()`

Example

- Read variants from a VCF file, and annotate with respect to a known
  gene model
  
```{r vcf, message=FALSE}
  ## 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(rowData(vcf),
      TxDb.Hsapiens.UCSC.hg19.knownGene,
      CodingVariants())
  head(coding)
```

Related packages

- [ensemblVEP][] 
  - Forward variants to Ensembl Variant Effect Predictor
- [VariantTools][], [h5vc][]
  - Call variants

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](http://bioinformatics.oxfordjournals.org/content/early/2014/04/21/bioinformatics.btu168)

### [rtracklayer][] (Genome annotations)

- Import BED, GTF, WIG, etc
- Export GRanges to BED, GTF, WIG, ...
- Access UCSC genome browser

## Big data

Restriction
- e.g., `ScanBamParam()` limits input to desired data at specific
  genomic ranges

Iteration
- e.g., `yieldSize` argument of `BamFile()`, or `FastqStreamer()`
  allows iteration through large files.

Compression
- Genomic vectors represented as `Rle` (run-length encoding) class
- Lists e.g., `GRangesList` are efficiently maintain the illusion that
  vector elements are grouped.

Parallel processing
- e.g., via [BiocParallel][] package

Reference

- Lawrence, M and Morgan, M. Scalable Genomic Computing and
  Visualization with R/Bioconductor. Statistical Science
  [forthcoming](http://www.e-publications.org/ims/submission/STS/user/submissionFile/16331?confirm=d92def1e)

## Exercises

### 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 alreaddy defined
   
```{r summarizeOverlaps-roi, message=FALSE}
   require(TxDb.Hsapiens.UCSC.hg19.knownGene)
   exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
   ## only chromosome 14
   seqlevels(exByGn, force=TRUE) = "chr14"
```

2. Identify the sample BAM files.

```{r summarizeOverlaps-bam, message=FALSE}
   require(RNAseqData.HNRNPC.bam.chr14)
   length(RNAseqData.HNRNPC.bam.chr14_BAMFILES)
```

3. Summarize overlaps, optionally in parallel

```{r summarizeOverlaps}
   ## 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.
   
```{r summarizeOverlaps-explore}
   olaps
   head(assay(olaps))
   colSums(assay(olaps))                # library sizes
   plot(sum(width(olaps)), rowMeans(assay(olaps)), log="xy")
```

5. As an advanced exercise, investigate the relationship between GC
   content and read count
   
```{r summarizeOverlaps-gc}
   require(BSgenome.Hsapiens.UCSC.hg19)
   sequences <- getSeq(BSgenome.Hsapiens.UCSC.hg19, rowData(olaps))
   gcPerExon <- letterFrequency(unlist(sequences), "GC")
   gc <- relist(as.vector(gcPerExon), sequences)
   gc_percent <- sum(gc) / sum(width(olaps))
   plot(gc_percent, rowMeans(assay(olaps)), log="y")
```

[BSgenome]: http://bioconductor.org/packages/release/bioc/html/BSgenome.html
[BiocParallel]: http://bioconductor.org/packages/release/bioc/html/BiocParallel.html
[Biostrings]: http://bioconductor.org/packages/release/bioc/html/Biostrings.html
[DESeq2]: http://bioconductor.org/packages/release/bioc/html/DESeq2.html
[GenomicRanges]: http://bioconductor.org/packages/release/bioc/html/GenomicRanges.html
[GenomicAlignments]: http://bioconductor.org/packages/release/bioc/html/GenomicAlignments.html
[Rsamtools]: http://bioconductor.org/packages/release/bioc/html/Rsamtools.html
[ShortRead]: http://bioconductor.org/packages/release/bioc/html/ShortRead.html
[VariantAnnotation]: http://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html
[VariantTools]: http://bioconductor.org/packages/release/bioc/html/VariantTools.html
[edgeR]: http://bioconductor.org/packages/release/bioc/html/edgeR.html
[ensemblVEP]: http://bioconductor.org/packages/release/bioc/html/ensemblVEP.html 
[h5vc]: http://bioconductor.org/packages/release/bioc/html/h5vc.html
[rtracklayer]: http://bioconductor.org/packages/release/bioc/html/rtracklayer.html
