Original Authors: Martin Morgan, Sonali Arora
Presenting Author: Martin Morgan (martin.morgan@roswellpark.org)
Date: 11 July, 2016 Back: Monday labs
Objective: Learn the essentials of Bioconductor data structures
Lessons learned:
This section focuses on classes, methods, and packages, with the goal being to learn to navigate the help system and interactive discovery facilities.
Sequence analysis is specialized
Additional considerations
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.
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 elementType
## Class: integer integer characterORNULL character
##
## Name: elementMetadata metadata
## Class: DataTableORNULL list
##
## Extends:
## Class "Ranges", directly
## Class "GRangesOrIRanges", directly
## Class "IntegerList", by class "Ranges", distance 2
## Class "RangesORmissing", by class "Ranges", distance 2
## Class "AtomicList", by class "Ranges", distance 3
## Class "List", by class "Ranges", distance 4
## Class "Vector", by class "Ranges", distance 5
## Class "Annotated", by class "Ranges", 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 elementMetadata
## Class: Rle IRanges Rle DataFrame
##
## Name: seqinfo metadata
## Class: Seqinfo list
##
## Extends:
## Class "GenomicRanges", directly
## Class "GRangesOrIRanges", directly
## Class "Vector", by class "GenomicRanges", distance 2
## Class "GenomicRangesORmissing", by class "GenomicRanges", distance 2
## Class "GenomicRangesORGRangesList", by class "GenomicRanges", distance 2
## Class "GenomicRangesORGenomicRangesList", by class "GenomicRanges", distance 2
## Class "Annotated", by class "GenomicRanges", distance 3
?"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")
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
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
...
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?########################
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
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
Input: rtracklayer import()
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";
...
...
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.
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
library(GenomicRanges)
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1) # 1-based coordinates!
## GRanges object with 3 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] A [11, 15] +
## [2] A [21, 25] +
## [3] A [23, 27] +
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
range(gr) # intra-range
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] A [10, 26] +
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
reduce(gr) # inter-range
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] A [10, 14] +
## [2] A [20, 26] +
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
coverage(gr)
## RleList of length 1
## $A
## integer-Rle of length 26 with 6 runs
## Lengths: 9 5 5 2 3 2
## Values : 0 1 0 1 2 1
setdiff(range(gr), gr) # 'introns'
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] A [15, 19] +
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
IRangesList, GRangesList - List: all elements of the same type - Many *List-aware methods, but a common ‘trick’: apply a vectorized function to the unlisted representaion, then re-list
grl <- GRangesList(...)
orig_gr <- unlist(grl)
transformed_gr <- FUN(orig)
transformed_grl <- relist(transformed_gr, grl)
Reference
Classes
Methods –
reverseComplement()
letterFrequency()
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.hg19)
chr14_range = GRanges("chr14", IRanges(1, seqlengths(Hsapiens)["chr14"]))
chr14_dna <- getSeq(Hsapiens, chr14_range)
letterFrequency(chr14_dna, "GC", as.prob=TRUE)
## G|C
## [1,] 0.336276
Classes – GenomicRanges-like behaivor
Methods
readGAlignments()
, readGAlignmentsList()
summarizeOverlaps()
Example
Find reads supporting the junction identified above, at position 19653707 + 66M = 19653773 of chromosome 14
library(GenomicRanges)
library(GenomicAlignments)
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 njunc
## [1] chr14 - 66M120N6M 72 19653707 19653898 192 1
## [2] chr14 + 7M1270N65M 72 19652348 19653689 1342 1
##
## [[2]]
## GAlignments object with 2 alignments and 0 metadata columns:
## seqnames strand cigar qwidth start end width njunc
## [1] chr14 - 66M120N6M 72 19653707 19653898 192 1
## [2] chr14 + 72M 72 19653686 19653757 72 0
##
## [[3]]
## GAlignments object with 2 alignments and 0 metadata columns:
## seqnames strand cigar qwidth start end width njunc
## [1] chr14 + 72M 72 19653675 19653746 72 0
## [2] chr14 - 65M120N7M 72 19653708 19653899 192 1
##
## ...
## <5 more elements>
## -------
## seqinfo: 93 sequences from an unspecified genome
Classes – GenomicRanges-like behavior
Functions and methods
readVcf()
, readGeno()
, readInfo()
, readGT()
, writeVcf()
, filterVcf()
locateVariants()
(variants overlapping ranges), predictCoding()
, summarizeVariants()
genotypeToSnpMatrix()
, snpSummary()
Example
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())
head(coding)
## GRanges object with 6 ranges and 9 metadata columns:
## seqnames ranges strand | LOCATION LOCSTART LOCEND
## <Rle> <IRanges> <Rle> | <factor> <integer> <integer>
## 1 chr22 [50301422, 50301422] - | coding 939 939
## 2 chr22 [50301476, 50301476] - | coding 885 885
## 3 chr22 [50301488, 50301488] - | coding 873 873
## 4 chr22 [50301494, 50301494] - | coding 867 867
## 5 chr22 [50301584, 50301584] - | coding 777 777
## 6 chr22 [50302962, 50302962] - | coding 698 698
## QUERYID TXID CDSID GENEID PRECEDEID
## <integer> <character> <IntegerList> <character> <CharacterList>
## 1 24 75253 218562 79087
## 2 25 75253 218562 79087
## 3 26 75253 218562 79087
## 4 27 75253 218562 79087
## 5 28 75253 218562 79087
## 6 57 75253 218563 79087
## FOLLOWID
## <CharacterList>
## 1
## 2
## 3
## 4
## 5
## 6
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
Related packages
Reference
Much Bioinformatic data is very large. The discussion so far has assumed that the data can be read into memory. Here we mention two important general strategies for working with large data; we will explore these in greater detail in a later lab, but feel free to ask questions and explore this material now.
Restriction
ScanBamParam()
limits input to desired data at specific genomic rangesIteration
yieldSize
argument of BamFile()
, or FastqStreamer()
allows iteration through large files.Compression
Rle
(run-length encoding) classGRangesList
are efficiently maintain the illusion that vector elements are grouped.Parallel processing
Reference
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.
Identify the regions of interest. We use a ‘TxDb’ package with gene models alreaddy defined
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
## only chromosome 14
seqlevels(exByGn, force=TRUE) = "chr14"
Identify the sample BAM files.
library(RNAseqData.HNRNPC.bam.chr14)
length(RNAseqData.HNRNPC.bam.chr14_BAMFILES)
## [1] 8
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)
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
As an advanced exercise, investigate the relationship between GC content and read count
library(BSgenome.Hsapiens.UCSC.hg19)
sequences <- getSeq(BSgenome.Hsapiens.UCSC.hg19, rowRanges(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 in xy.coords(x, y, xlabel, ylabel, log): 252 y values <= 0 omitted
## from logarithmic plot