Lab 1.4: Data Representations
22 July 2019
Contents
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:
- Review of common bioinformatic file formats: FASTA, FASTQ, SAM / BAM, VCF, BED / WIG / GTF
- How to work with Bioconductor objects for DNA sequences, genomic ranges, aligned reads, and called variants.
1 Classes, methods, and packages
This section focuses on classes, methods, and packages, with the goal being to learn to navigate the help system and interactive discovery facilities.
1.1 Motivation
Sequence analysis is specialized
- Large data needs to be processed in a memory- and time-efficient manner
- Specific algorithms have been developed for the unique characteristics of sequence data
Additional considerations
- Re-use of existing, tested code is easier to do and less error-prone than re-inventing the wheel.
- Interoperability between packages is easier when the packages share similar data structures.
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.
1.2 Objects
Load the Biostrings and GenomicRanges package
library(Biostrings)
library(GenomicRanges)- Bioconductor makes extensive use of classes to represent complicated data types
- Classes foster interoperability – many different packages can work on the same data – but can be a bit intimidating for the user.
- Formal ‘S4’ object system
- Often a class is described on a particular home page, e.g.,
?GRanges, and in vignettes, e.g.,vignette(package="GenomicRanges"),vignette("GenomicRangesIntroduction") - Many methods and classes can be discovered interactively , e.g.,
methods(class="GRanges")to find out what one can do with aGRangesinstance, andmethods(findOverlaps)for classes that thefindOverlaps()function operates on. - In more advanced cases, one can look at the actual definition of
a class or method using
getClass(),getMethod()
- Often a class is described on a particular home page, e.g.,
- Interactive help
?findOverlaps,<tab>to select help on a specific method,?GRanges-classfor help on a class.
1.3 Example & short exercise: Biostrings
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 NEB03m <- 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 1methods(class=class(phiX174Phage)) # 'DNAStringSet' methodsExercises
- 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.
- Discover vignettes in the Biostrings package with
vignette(package="Biostrings"). Add another argument to thevignettefunction to view the ‘BiostringsQuickOverview’ vignette. - If the internet is available, 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?
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"
2 Working with genomic ranges
2.1 IRanges and GRanges
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 5There 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 3The 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
##
## 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 seqlengthsThe 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 seqlengthsDiscover 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_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
##
## Known Subclasses: "GPos"?"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")2.2 Range-based operations
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)3 How to use high-throughput sequence data types in Bioconductor
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
3.1 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
...3.2 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
3.3 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.fastqAlignments: 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
3.4 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
3.5 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"; ... ...
4 Exercises: Data representation in R / Bioconductor
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.
4.1 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
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
BSgenomepackages. 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
- Mouse CDS sequence, from Ensembl: https://useast.ensembl.org/info/data/ftp/index.html
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)## Using temporary cache /tmp/RtmpIzVK40/BiocFileCache## adding rname 'ftp://ftp.ensembl.org/pub/release-92/fasta/mus_musculus/cds/Mus_musculus.GRCm38.cds.all.fa.gz'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
ATGor end with a stop codonc("TAA", "TAG", "TGA")pred1 <- width(cds) %% 3 == 0 table(pred1)## pred1 ## FALSE TRUE ## 7219 58251pred2 <- 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 48662cds <- 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.5714286hist(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) –
DNAStringSetinherits fromVectorandAnnotated, which means that each element (sequence) can have additional information, for instance we can associate GC content with each sequencemcols(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.596269554753309mcols(cds[1:3], use.names = FALSE)## DataFrame with 3 rows and 1 column ## GC ## <numeric> ## 1 0.502662406815761 ## 2 0.502662406815761 ## 3 0.502662406815761
4.2 SummarizedExperiment
Motivation: reproducible & interoperable
- Matrix of feature x sample measurements,
assays() Addition description about samples,
colData()- Covariates, e.g., age, gender
- Experimental design, e.g., treatment group
Additional information about features,
rowData()- Gene names, width, GC content, …
- Genomic ranges(!)
- Derived values, E.g., log-fold change between treatments, P-value, …
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 object is masked from 'package:Biostrings': ## ## type## The following objects are masked from 'package:base': ## ## aperm, apply, rowsumdata(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 BioSampleExplore 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 SAMN02422677airway[ , 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 BioSampleCalculate 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
- What’s the frequency of 0 counts in the single cell assay data?
- What’s the distribution of library sizes in the single cell assay data?
- Create a random sample of 100 cells and visualize the relationship
between samples using
dist()andcmdscale(). - can you identify what column of
colData()is responsible for any pattern you see? - In exploring the covariates, are the possible problems with confounding?
4.3 GenomicRanges
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 seqlengthsrange(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 seqlengthsreduce(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 seqlengthscoverage(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 1setdiff(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 seqlengthsExercises
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 FALSEWhich 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 2precede(reg, genes)## [1] NA 2What 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 0as(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
- 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
4.4 GenomicAlignments (Aligned reads)
Classes – GenomicRanges-like behaivor
- GAlignments, GAlignmentPairs, GAlignmentsList
Methods
readGAlignments(),readGAlignmentsList()- Easy to restrict input, iterate in chunks
summarizeOverlaps()
Exercises
Find reads supporting the junction identified above, at position 19653707 + 66M = 19653773 of chromosome 14
library(GenomicRanges) library(GenomicAlignments)## Loading required package: Rsamtoolslibrary(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
4.5 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()
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 QUERYID ## <Rle> <IRanges> <Rle> | <factor> <integer> <integer> <integer> ## [1] chr22 50301422 - | coding 939 939 24 ## [2] chr22 50301476 - | coding 885 885 25 ## [3] chr22 50301488 - | coding 873 873 26 ## [4] chr22 50301494 - | coding 867 867 27 ## [5] chr22 50301584 - | coding 777 777 28 ## [6] chr22 50302962 - | coding 698 698 57 ## TXID CDSID GENEID PRECEDEID FOLLOWID ## <character> <IntegerList> <character> <CharacterList> <CharacterList> ## [1] 75253 218562 79087 <NA> <NA> ## [2] 75253 218562 79087 <NA> <NA> ## [3] 75253 218562 79087 <NA> <NA> ## [4] 75253 218562 79087 <NA> <NA> ## [5] 75253 218562 79087 <NA> <NA> ## [6] 75253 218563 79087 <NA> <NA> ## ------- ## seqinfo: 1 sequence from an unspecified genome; no seqlengths
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
4.6 rtracklayer (Genome annotations)
- Import BED, GTF, WIG, etc
- Export GRanges to BED, GTF, WIG, …
- Access UCSC genome browser
5 Extended exercises
5.1 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.
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] 8Summarize 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 0colSums(assay(olaps)) # library sizes## ERR127306 ERR127307 ERR127308 ERR127309 ERR127302 ERR127303 ERR127304 ## 340646 373268 371639 331518 313800 331135 331606 ## ERR127305 ## 329647plot(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': ## ## Hsapienssequences <- 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
6 End matter
6.1 Session Info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Linux Mint 18.3
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## 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.36.4
## [4] AnnotationDbi_1.46.0
## [5] VariantAnnotation_1.30.1
## [6] RNAseqData.HNRNPC.bam.chr14_0.22.0
## [7] GenomicAlignments_1.20.1
## [8] Rsamtools_2.0.0
## [9] airway_1.4.0
## [10] SummarizedExperiment_1.14.0
## [11] DelayedArray_0.10.0
## [12] BiocParallel_1.18.0
## [13] matrixStats_0.54.0
## [14] Biobase_2.44.0
## [15] BSgenome.Hsapiens.UCSC.hg38_1.4.1
## [16] BSgenome_1.52.0
## [17] rtracklayer_1.44.0
## [18] GenomicRanges_1.36.0
## [19] GenomeInfoDb_1.20.0
## [20] Biostrings_2.52.0
## [21] XVector_0.24.0
## [22] IRanges_2.18.1
## [23] S4Vectors_0.22.0
## [24] BiocGenerics_0.30.0
## [25] BiocStyle_2.12.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 lattice_0.20-38 prettyunits_1.0.2
## [4] assertthat_0.2.1 zeallot_0.1.0 digest_0.6.20
## [7] R6_2.4.0 backports_1.1.4 RSQLite_2.1.1
## [10] evaluate_0.14 httr_1.4.0 pillar_1.4.2
## [13] zlibbioc_1.30.0 rlang_0.4.0 progress_1.2.2
## [16] blob_1.2.0 Matrix_1.2-17 rmarkdown_1.14
## [19] stringr_1.4.0 RCurl_1.95-4.12 bit_1.1-14
## [22] biomaRt_2.40.3 compiler_3.6.1 xfun_0.8
## [25] pkgconfig_2.0.2 htmltools_0.3.6 tibble_2.1.3
## [28] GenomeInfoDbData_1.2.1 bookdown_0.12 codetools_0.2-16
## [31] XML_3.98-1.20 crayon_1.3.4 bitops_1.0-6
## [34] grid_3.6.1 DBI_1.0.0 magrittr_1.5
## [37] stringi_1.4.3 vctrs_0.2.0 tools_3.6.1
## [40] bit64_0.9-7 hms_0.5.0 yaml_2.2.0
## [43] BiocManager_1.30.4 memoise_1.1.0 knitr_1.236.2 Acknowledgements
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)