Sequence Data Representations

Epigenomics 2014
Author: Martin Morgan (mtmorgan@fhcrc.org)
Date: 24 August, 2014

Bioconductor for sequence analysis

Alt Sequencing Ecosystem

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

Example

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

    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"
    

Aligned reads: BAM files (e.g., ERR127306_chr14.bam)

Input & manipulation: 'low-level' Rsamtools, scanBam(), BamFile(); 'high-level' GenomicAlignments

Called variants: VCF files

Input and manipulation: VariantAnnotation readVcf(), readInfo(), readGeno() selectively with ScanVcfParam().

Genome annotations: BED, WIG, GTF, etc. files

Input: rtracklayer import()

Data representation

Ranges

Ranges represent:

Many common biological questions are range-based

The GenomicRanges package defines essential classes and methods

Alt

Alt

Range operations

Alt Ranges Algebra

Ranges

Intra-range methods

Inter-range methods

Between-range methods

Example

require(GenomicRanges)
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1)                            # 1-based coordinates!
## GRanges 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]      +
##   ---
##   seqlengths:
##     A
##    NA
range(gr)                               # intra-range
## GRanges with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [10, 26]      +
##   ---
##   seqlengths:
##     A
##    NA
reduce(gr)                              # inter-range
## GRanges with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [10, 14]      +
##   [2]        A  [20, 26]      +
##   ---
##   seqlengths:
##     A
##    NA
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 with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]        A  [15, 19]      +
##   ---
##   seqlengths:
##     A
##    NA

IRangesList, GRangesList

Reference

Biostrings (DNA or amino acid sequences)

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

GenomicAlignments (Aligned reads)

Classes – GenomicRanges-like behaivor

Methods

Example

  require(GenomicRanges)
  require(GenomicAlignments)
## Loading required package: GenomicAlignments
## Loading required package: Rsamtools
## 
## Attaching package: 'GenomicAlignments'
## 
## The following objects are masked from 'package:locfit':
## 
##     left, right
  require(Rsamtools)

  ## our 'region of interest'
  roi <- GRanges("chr14", IRanges(19653773, width=1)) 
  ## sample data
  require('RNAseqData.HNRNPC.bam.chr14')
## Loading required package: 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 of length 8:
## [[1]] 
## GAlignments 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 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 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>
## ---
## seqlengths:
##                   chr1                 chr10 ...                  chrY
##              249250621             135534747 ...              59373566

VariantAnnotation (Called variants)

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(rowData(vcf),
      TxDb.Hsapiens.UCSC.hg19.knownGene,
      CodingVariants())
  head(coding)
## GRanges with 6 ranges and 7 metadata columns:
##       seqnames               ranges strand | LOCATION   QUERYID      TXID
##          <Rle>            <IRanges>  <Rle> | <factor> <integer> <integer>
##   [1]    chr22 [50301422, 50301422]      - |   coding        24     75253
##   [2]    chr22 [50301476, 50301476]      - |   coding        25     75253
##   [3]    chr22 [50301488, 50301488]      - |   coding        26     75253
##   [4]    chr22 [50301494, 50301494]      - |   coding        27     75253
##   [5]    chr22 [50301584, 50301584]      - |   coding        28     75253
##   [6]    chr22 [50302962, 50302962]      - |   coding        57     75253
##           CDSID      GENEID       PRECEDEID        FOLLOWID
##       <integer> <character> <CharacterList> <CharacterList>
##   [1]    218562       79087                                
##   [2]    218562       79087                                
##   [3]    218562       79087                                
##   [4]    218562       79087                                
##   [5]    218562       79087                                
##   [6]    218563       79087                                
##   ---
##   seqlengths:
##    chr22
##       NA

Related packages

Reference

rtracklayer (Genome annotations)

Big data

Restriction

Iteration

Compression

Parallel processing

Reference

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
   require(TxDb.Hsapiens.UCSC.hg19.knownGene)
   exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
   ## only chromosome 14
   seqlevels(exByGn, force=TRUE) = "chr14"
  1. Identify the sample BAM files.
   require(RNAseqData.HNRNPC.bam.chr14)
   length(RNAseqData.HNRNPC.bam.chr14_BAMFILES)
## [1] 8
  1. 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)
  1. Explore our handiwork, e.g., library sizes (column sums), relationship between gene length and number of mapped reads, etc.
   olaps
## class: SummarizedExperiment 
## dim: 779 8 
## exptData(0):
## assays(1): counts
## rownames(779): 10001 100113389 ... 9950 9985
## rowData metadata column 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: 252 y values <= 0 omitted from logarithmic plot

plot of chunk summarizeOverlaps-explore

  1. As an advanced exercise, investigate the relationship between GC content and read count
   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")
## Warning: 252 y values <= 0 omitted from logarithmic plot

plot of chunk summarizeOverlaps-gc