Sequence Data Representations

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

Bioconductor for sequence analysis

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

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

• 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

• GRangesList

Range operations

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

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

  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

• 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

  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

• 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

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

• 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

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

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

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