Version: 0.1.1
 Compiled: Sat Oct 24 10:12:10 2015
Objectives
| Time | Topic | 
|---|---|
| 09:15 - 10:15 | Sequencing work flows and file types | 
| 10:15 | Tea/Coffee break | 
| 10:30 - 12:30 | Introduction to R and Bioconductor | 
| 12:30 | Lunch | 
| 13:30 -14:00 | Scalable computing | 
Wet-lab sequence preparation (figure from http://rnaseq.uoregon.edu/)
(Illumina) Sequencing (Bentley et al., 2008, doi:10.1038/nature07517
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
...Whole genomes: 2bit and .fa formats: rtracklayer, Rsamtools; BSgenome
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.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:1Input 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:3Input: 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";
      ...Language and environment for statistical computing and graphics
factor(), NAVector, class, object
logical, integer, numeric, complex, character, bytematrix – atomic vector with ‘dim’ attributedata.frame – list of equal length atomic vectorslm(), belowFunction, generic, method
rnorm(1000)print().print.factor; methods are invoked indirectly, via the generic.Introspection
class(), str()dim()Help
?print: help on the generic print?print.data.frame: help on print method for objects of class data.frame.Example
x <- rnorm(1000)                   # atomic vectors
y <- x + rnorm(1000, sd=.5)
df <- data.frame(x=x, y=y)         # object of class 'data.frame'
plot(y ~ x, df)                    # generic plot, method plot.formulafit <- lm(y ~x, df)                # object of class 'lm'
methods(class=class(fit))          # introspection##  [1] add1           alias          anova          case.names     coerce         confint       
##  [7] cooks.distance deviance       dfbeta         dfbetas        drop1          dummy.coef    
## [13] effects        extractAIC     family         formula        hatvalues      influence     
## [19] initialize     kappa          labels         logLik         model.frame    model.matrix  
## [25] nobs           plot           predict        print          proj           qr            
## [31] residuals      rstandard      rstudent       show           simulate       slotsFromS3   
## [37] summary        variable.names vcov          
## see '?methods' for accessing help and source codeAnalysis and comprehension of high-throughput genomic data
Packages, vignettes, work flows
Objects
methods(), getClass(), selectMethod()method?"substr,<tab>" to select help on methods, class?D<tab> for help on classesExample
require(Biostrings)                     # Biological sequences
data(phiX174Phage)                      # sample data, see ?phiX174Phage
phiX174Phage##   A DNAStringSet instance of length 6
##     width seq                                                                   names               
## [1]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA Genbank
## [2]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA RF70s
## [3]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA SS78
## [4]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA Bull
## [5]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA G97
## [6]  5386 GAGTTTTATCGCTTCCATGACGCAGAAGTTAAC...TTCGATAAAAATGATTGGCGTATCCAACCTGCA 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))##   [1] !                               !=                             
##   [3] [                               [[                             
##   [5] [[<-                            [<-                            
##   [7] %in%                            <                              
##   [9] <=                              ==                             
##  [11] >                               >=                             
##  [13] $                               $<-                            
##  [15] aggregate                       alphabetFrequency              
##  [17] anyNA                           append                         
##  [19] as.character                    as.complex                     
##  [21] as.data.frame                   as.env                         
##  [23] as.integer                      as.list                        
##  [25] as.logical                      as.matrix                      
##  [27] as.numeric                      as.raw                         
##  [29] as.vector                       c                              
##  [31] chartr                          coerce                         
##  [33] compact                         compare                        
##  [35] compareStrings                  complement                     
##  [37] consensusMatrix                 consensusString                
##  [39] countOverlaps                   countPattern                   
##  [41] countPDict                      dinucleotideFrequencyTest      
##  [43] do.call                         droplevels                     
##  [45] duplicated                      elementLengths                 
##  [47] elementMetadata                 elementMetadata<-              
##  [49] elementType                     endoapply                      
##  [51] eval                            expand                         
##  [53] extractAt                       extractROWS                    
##  [55] Filter                          Find                           
##  [57] findOverlaps                    hasOnlyBaseLetters             
##  [59] head                            high2low                       
##  [61] ifelse                          intersect                      
##  [63] is.na                           is.unsorted                    
##  [65] isEmpty                         isMatchingEndingAt             
##  [67] isMatchingStartingAt            lapply                         
##  [69] length                          lengths                        
##  [71] letterFrequency                 Map                            
##  [73] match                           matchPattern                   
##  [75] matchPDict                      mcols                          
##  [77] mcols<-                         mendoapply                     
##  [79] metadata                        metadata<-                     
##  [81] mstack                          names                          
##  [83] names<-                         narrow                         
##  [85] nchar                           neditEndingAt                  
##  [87] neditStartingAt                 NROW                           
##  [89] nucleotideFrequencyAt           oligonucleotideFrequency       
##  [91] order                           overlapsAny                    
##  [93] PairwiseAlignments              PairwiseAlignmentsSingleSubject
##  [95] parallelSlotNames               PDict                          
##  [97] Position                        PWM                            
##  [99] rank                            Reduce                         
## [101] relist                          relistToClass                  
## [103] rename                          rep                            
## [105] rep.int                         replaceAt                      
## [107] replaceLetterAt                 replaceROWS                    
## [109] rev                             revElements                    
## [111] reverse                         reverseComplement              
## [113] ROWNAMES                        sapply                         
## [115] seqinfo                         seqinfo<-                      
## [117] seqlevelsInUse                  seqtype                        
## [119] seqtype<-                       setdiff                        
## [121] setequal                        shiftApply                     
## [123] show                            showAsCell                     
## [125] sort                            split                          
## [127] split<-                         splitAsList                    
## [129] stack                           stringDist                     
## [131] subseq                          subseq<-                       
## [133] subset                          subsetByOverlaps               
## [135] table                           tail                           
## [137] tapply                          threebands                     
## [139] toString                        translate                      
## [141] trimLRPatterns                  twoWayAlphabetFrequency        
## [143] union                           unique                         
## [145] uniqueLetters                   unlist                         
## [147] unsplit                         unstrsplit                     
## [149] updateObject                    values                         
## [151] values<-                        vcountPattern                  
## [153] vcountPDict                     vmatchPattern                  
## [155] vwhichPDict                     which.isMatchingEndingAt       
## [157] which.isMatchingStartingAt      whichPDict                     
## [159] width                           window                         
## [161] window<-                        with                           
## [163] within                          xtfrm                          
## [165] xvcopy                         
## see '?methods' for accessing help and source codeselectMethod(reverseComplement, class(phiX174Phage))## Method Definition:
## 
## function (x, ...) 
## xvcopy(x, lkup = getDNAComplementLookup(), reverse = TRUE)
## <environment: namespace:Biostrings>
## 
## Signatures:
##         x             
## target  "DNAStringSet"
## defined "DNAStringSet"This very open-ended topic points to some of the most prominent Bioconductor packages for sequence analysis. Use the opportunity in this lab to explore the package vignettes and help pages highlighted below; many of the material will be covered in greater detail in subsequent labs and lectures.
Basics
A package needs to be installed once, using the instructions on the landing page. Once installed, the package can be loaded into an R session
library(GenomicRanges)and the help system queried interactively, as outlined above:
  help(package="GenomicRanges")
  vignette(package="GenomicRanges")
  vignette(package="GenomicRanges", "GenomicRangesHOWTOs")
  ?GRangesDomain-specific analysis – explore the landing pages, vignettes, and reference manuals of two or three of the following packages.
Working with sequences, alignments, common web file formats, and raw data; these packages rely very heavily on the IRanges / GenomicRanges infrastructure that we will encounter later in the course.
?consensusMatrix, for instance. Also check out the BSgenome package for working with whole genome sequences, e.g., ?"getSeq,BSgenome-method"?readGAlignments help page and vigentte(package="GenomicAlignments",   "summarizeOverlaps")import and export functions can read in many common file types, e.g., BED, WIG, GTF, …, in addition to querying and navigating the UCSC genome browser. Check out the ?import page for basic usage.Visualization
Classes
Methods –
reverseComplement()letterFrequency()matchPDict(), matchPWM()Related packages
Example
BSgenome packages. The following calculates GC content across chr14.  require(BSgenome.Hsapiens.UCSC.hg19)
  chr14_range = GRanges("chr14", IRanges(1, seqlengths(Hsapiens)["chr14"]))
  chr14_dna <- getSeq(Hsapiens, chr14_range)
  letterFrequency(chr14_dna, "GC", as.prob=TRUE)##           G|C
## [1,] 0.336276Ranges represent: - Data, e.g., aligned reads, ChIP peaks, SNPs, CpG islands, … - Annotations, e.g., gene models, regulatory elements, methylated regions - Ranges are defined by chromosome, start, end, and strand - Often, metadata is associated with each range, e.g., quality of alignment, strength of ChIP peak
Many common biological questions are range-based - What reads overlap genes? - What genes are ChIP peaks nearest? - …
The GenomicRanges package defines essential classes and methods
GRangesGRangesListRanges - IRanges - start() / end() / width() - List-like – length(), subset, etc. - ‘metadata’, mcols() - GRanges - ‘seqnames’ (chromosome), ‘strand’ - Seqinfo, including seqlevels and seqlengths
Intra-range methods - Independent of other ranges in the same object - GRanges variants strand-aware - shift(), narrow(), flank(), promoters(), resize(), restrict(), trim() - See ?"intra-range-methods"
Inter-range methods - Depends on other ranges in the same object - range(), reduce(), gaps(), disjoin() - coverage() (!) - see ?"inter-range-methods"
Between-range methods - Functions of two (or more) range objects - findOverlaps(), countOverlaps(), …, %over%, %within%, %outside%; union(), intersect(), setdiff(), punion(), pintersect(), psetdiff()
Example
require(GenomicRanges)
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1)                            # 1-based coordinates!## GRanges 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 seqlengthsrange(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 seqlengthsreduce(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 seqlengthscoverage(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 1setdiff(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 seqlengthsIRangesList, GRangesList - List: all elements of the same type - Many *List-aware methods, but a common ‘trick’: apply a vectorized function to the unlisted representaion, then re-list
    grl <- GRangesList(...)
    orig_gr <- unlist(grl)
    transformed_gr <- FUN(orig)
    transformed_grl <- relist(, grl)
    Reference
Classes – GenomicRanges-like behaivor
Methods
readGAlignments(), readGAlignmentsList()summarizeOverlaps()Example
require(GenomicRanges)
require(GenomicAlignments)## Loading required package: GenomicAlignments
## Loading required package: Rsamtoolsrequire(Rsamtools)
## our 'region of interest'
roi <- GRanges("chr14", IRanges(19653773, width=1)) 
## sample data
require('RNAseqData.HNRNPC.bam.chr14')## Loading required package: RNAseqData.HNRNPC.bam.chr14bf <- 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 genomeClasses – GenomicRanges-like behavior
Functions and methods
readVcf(), readGeno(), readInfo(), readGT(), writeVcf(), filterVcf()locateVariants() (variants overlapping ranges), predictCoding(), summarizeVariants()genotypeToSnpMatrix(), snpSummary()Example
  ## input variants
  require(VariantAnnotation)
  fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
  vcf <- readVcf(fl, "hg19")
  seqlevels(vcf) <- "chr22"
  ## known gene model
  require(TxDb.Hsapiens.UCSC.hg19.knownGene)
  coding <- locateVariants(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   QUERYID        TXID
##        <Rle>            <IRanges>  <Rle> | <factor> <integer> <integer> <integer> <character>
##   1    chr22 [50301422, 50301422]      - |   coding       939       939        24       75253
##   2    chr22 [50301476, 50301476]      - |   coding       885       885        25       75253
##   3    chr22 [50301488, 50301488]      - |   coding       873       873        26       75253
##   4    chr22 [50301494, 50301494]      - |   coding       867       867        27       75253
##   5    chr22 [50301584, 50301584]      - |   coding       777       777        28       75253
##   6    chr22 [50302962, 50302962]      - |   coding       698       698        57       75253
##             CDSID      GENEID       PRECEDEID        FOLLOWID
##     <IntegerList> <character> <CharacterList> <CharacterList>
##   1        218562       79087                                
##   2        218562       79087                                
##   3        218562       79087                                
##   4        218562       79087                                
##   5        218562       79087                                
##   6        218563       79087                                
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengthsRelated packages
Reference
assays()colData() data frame for desciption of samplesrowRanges() GRanges / GRangeList or data frame for description of featuresexptData() to describe the entire object
library(SummarizedExperiment)
library(airway)
data(airway)
airway## class: RangedSummarizedExperiment 
## dim: 64102 8 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSamplecolData(airway)## DataFrame with 8 rows and 9 columns
##            SampleName     cell      dex    albut        Run avgLength Experiment    Sample
##              <factor> <factor> <factor> <factor>   <factor> <integer>   <factor>  <factor>
## SRR1039508 GSM1275862   N61311    untrt    untrt SRR1039508       126  SRX384345 SRS508568
## SRR1039509 GSM1275863   N61311      trt    untrt SRR1039509       126  SRX384346 SRS508567
## SRR1039512 GSM1275866  N052611    untrt    untrt SRR1039512       126  SRX384349 SRS508571
## SRR1039513 GSM1275867  N052611      trt    untrt SRR1039513        87  SRX384350 SRS508572
## SRR1039516 GSM1275870  N080611    untrt    untrt SRR1039516       120  SRX384353 SRS508575
## SRR1039517 GSM1275871  N080611      trt    untrt SRR1039517       126  SRX384354 SRS508576
## SRR1039520 GSM1275874  N061011    untrt    untrt SRR1039520       101  SRX384357 SRS508579
## SRR1039521 GSM1275875  N061011      trt    untrt SRR1039521        98  SRX384358 SRS508580
##               BioSample
##                <factor>
## SRR1039508 SAMN02422669
## SRR1039509 SAMN02422675
## SRR1039512 SAMN02422678
## SRR1039513 SAMN02422670
## SRR1039516 SAMN02422682
## SRR1039517 SAMN02422673
## SRR1039520 SAMN02422683
## SRR1039521 SAMN02422677airway[, airway$dex %in% "trt"]## class: RangedSummarizedExperiment 
## dim: 64102 4 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample?select?exonsBy page to retrieve all exons grouped by gene or transcript.open(), read chunk(s), close().yieldSize argument to Rsamtools::BamFile()Rsamtools::ScanBamParam()ShortRead::FastqSampler()lapply()-like operationsParallel evaluation in Bioconductor
bplapply() for lapply()-like functions, increasingly used by package developers to provide easy, standard way of gaining parallel evaluation.R / Bioconductor
Publications (General Bioconductor)
Other