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:
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.
Load the Biostrings and GenomicRanges package
library(Biostrings)
library(GenomicRanges)
?GRanges
, and in vignettes, e.g.,
vignette(package="GenomicRanges")
,
vignette("GenomicRangesIntroduction")
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.getClass()
, getMethod()
?findOverlaps,<tab>
to select help on a specific method,
?GRanges-class
for help on a class.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
vignette(package="Biostrings")
. Add another argument to the
vignette
function to view the ‘BiostringsQuickOverview’ vignette.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"
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")
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)
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
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.
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.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
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")
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 ]
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]"
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=".")
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
(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
Motivation: reproducible & interoperable
assays()
Addition description about samples, colData()
Additional information about features, rowData()
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
dist()
and cmdscale()
.colData()
is responsible for any
pattern you see?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
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
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
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
Classes – GenomicRanges-like behaivor
Methods
readGAlignments()
, readGAlignmentsList()
summarizeOverlaps()
Exercises
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>
Classes – GenomicRanges-like behavior
Functions and methods
readVcf()
, readGeno()
, readInfo()
,
readGT()
, writeVcf()
, filterVcf()
locateVariants()
(variants overlapping ranges),
predictCoding()
, summarizeVariants()
genotypeToSnpMatrix()
, snpSummary()
Exerises
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
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 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"
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)
##
## Attaching package: 'BSgenome.Hsapiens.UCSC.hg19'
## The following object is masked from 'package:BSgenome.Hsapiens.UCSC.hg38':
##
## Hsapiens
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
sessionInfo()
## R version 3.6.1 Patched (2019-07-16 r76845)
## Platform: x86_64-apple-darwin17.7.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Users/ma38727/bin/R-3-6-branch/lib/libRblas.dylib
## LAPACK: /Users/ma38727/bin/R-3-6-branch/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.0
## [2] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [3] GenomicFeatures_1.37.4
## [4] AnnotationDbi_1.47.0
## [5] VariantAnnotation_1.31.3
## [6] RNAseqData.HNRNPC.bam.chr14_0.23.0
## [7] GenomicAlignments_1.21.4
## [8] Rsamtools_2.1.3
## [9] airway_1.5.0
## [10] SummarizedExperiment_1.15.5
## [11] DelayedArray_0.11.4
## [12] BiocParallel_1.19.0
## [13] matrixStats_0.54.0
## [14] Biobase_2.45.0
## [15] BSgenome.Hsapiens.UCSC.hg38_1.4.1
## [16] BSgenome_1.53.0
## [17] rtracklayer_1.45.1
## [18] GenomicRanges_1.37.14
## [19] GenomeInfoDb_1.21.1
## [20] Biostrings_2.53.2
## [21] XVector_0.25.0
## [22] IRanges_2.19.10
## [23] S4Vectors_0.23.17
## [24] BiocGenerics_0.31.5
## [25] BiocStyle_2.13.2
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 bit64_0.9-7 assertthat_0.2.1
## [4] askpass_1.1 BiocManager_1.30.4 BiocFileCache_1.9.1
## [7] blob_1.2.0 GenomeInfoDbData_1.2.1 yaml_2.2.0
## [10] progress_1.2.2 pillar_1.4.2 RSQLite_2.1.1
## [13] backports_1.1.4 lattice_0.20-38 glue_1.3.1
## [16] digest_0.6.20 htmltools_0.3.6 Matrix_1.2-17
## [19] XML_3.98-1.20 pkgconfig_2.0.2 biomaRt_2.41.6
## [22] bookdown_0.12 zlibbioc_1.31.0 purrr_0.3.2
## [25] tibble_2.1.3 openssl_1.4 magrittr_1.5
## [28] crayon_1.3.4 memoise_1.1.0 evaluate_0.14
## [31] tools_3.6.1 prettyunits_1.0.2 hms_0.5.0
## [34] stringr_1.4.0 compiler_3.6.1 rlang_0.4.0
## [37] grid_3.6.1 RCurl_1.95-4.12 rappdirs_0.3.1
## [40] bitops_1.0-6 rmarkdown_1.14 codetools_0.2-16
## [43] DBI_1.0.0 curl_3.3 R6_2.4.0
## [46] knitr_1.23 dplyr_0.8.3 bit_1.1-14
## [49] zeallot_0.1.0 stringi_1.4.3 Rcpp_1.0.1
## [52] vctrs_0.2.0 dbplyr_1.4.2 tidyselect_0.2.5
## [55] xfun_0.8
Research reported in this tutorial was supported by the National Human Genome Research Institute and the National Cancer Institute of the National Institutes of Health under award numbers U41HG004059 and U24CA180996.
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 633974)