← Back to branch summary

~ubuntu-branches/ubuntu/wily/r-bioc-genomicranges/wily-proposed

« back to all changes in this revision

Viewing changes to vignettes/GenomicRangesIntroduction.Rnw

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2014-06-13 15:04:19 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20140613150419-v49mxnlg42rnuks5
Tags: 1.16.3-1
* New upstream version
* New (Build-)Depends: r-bioc-genomeinfodb
* cme fix dpkg-control
* add autopkgtest

Show diffs side-by-side

added added

removed removed

Lines of Context:
vignettes/GenomicRangesIntroduction.Rnw
 
1
%\VignetteIndexEntry{An Introduction to GenomicRanges}
 
2
%\VignetteDepends{Rsamtools}
 
3
%\VignetteKeywords{sequence, sequencing, alignments}
 
4
%\VignettePackage{GenomicRanges}
 
5
\documentclass[10pt]{article}
 
6
 
 
7
\usepackage{times}
 
8
\usepackage{hyperref}
 
9
 
 
10
\textwidth=6.5in
 
11
\textheight=8.5in
 
12
%\parskip=.3cm
 
13
\oddsidemargin=-.1in
 
14
\evensidemargin=-.1in
 
15
\headheight=-.3in
 
16
 
 
17
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
 
18
\newcommand{\Robject}[1]{{\texttt{#1}}}
 
19
\newcommand{\Rpackage}[1]{{\textit{#1}}}
 
20
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
 
21
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
 
22
\newcommand{\Rclass}[1]{{\textit{#1}}}
 
23
\newcommand{\Rcode}[1]{{\texttt{#1}}}
 
24
 
 
25
\newcommand{\software}[1]{\textsf{#1}}
 
26
\newcommand{\R}{\software{R}}
 
27
\newcommand{\Bioconductor}{\software{Bioconductor}}
 
28
\newcommand{\GenomicRanges}{\Rpackage{GenomicRanges}}
 
29
 
 
30
 
 
31
\title{An Introduction to Genomic Ranges Classes}
 
32
\author{Marc Carlson \and Patrick Aboyoun \and Herv\'{e} Pag\`{e}s}
 
33
\date{\today}
 
34
 
 
35
\begin{document}
 
36
 
 
37
\maketitle
 
38
 
 
39
<<options,echo=FALSE>>=
 
40
options(width=72)
 
41
@
 
42
\tableofcontents
 
43
 
 
44
\section{Introduction}
 
45
 
 
46
The \Rpackage{GenomicRanges} package serves as the foundation for
 
47
representing genomic locations within the \software{Bioconductor}
 
48
project.  In the \software{Bioconductor} package hierarchy, it builds
 
49
upon the \Rpackage{IRanges} (infrastructure) package and provides
 
50
support for the \Rpackage{BSgenome} (infrastructure),
 
51
\Rpackage{Rsamtools} (I/O), \Rpackage{ShortRead} (I/O \& QA),
 
52
\Rpackage{rtracklayer} (I/O), and \Rpackage{GenomicFeatures}
 
53
(infrastructure) packages.
 
54
 
 
55
This package lays a foundation for genomic analysis by introducing
 
56
three classes (\Rclass{GRanges}, \Rclass{GRangesList}, and
 
57
\Rclass{GAlignments}), which are used to represent single
 
58
interval range features, multiple interval range features, and genomic
 
59
alignments respectively.  This vignette focuses on these classes and
 
60
their associated methods.
 
61
 
 
62
The \Rpackage{GenomicRanges} package is available at bioconductor.org
 
63
and can be downloaded via \Rfunction{biocLite}:
 
64
 
 
65
<<biocLite, eval=FALSE>>=
 
66
source("http://bioconductor.org/biocLite.R")
 
67
biocLite("GenomicRanges")
 
68
@
 
69
<<initialize, results=hide>>=
 
70
library(GenomicRanges)
 
71
@
 
72
 
 
73
 
 
74
\section{\Rclass{GRanges}: Single Interval Range Features}
 
75
 
 
76
The \Rclass{GRanges} class represents a collection of genomic features
 
77
that each have a single start and end location on the genome. This
 
78
includes features such as contiguous binding sites, transcripts, and
 
79
exons. These objects can be created by using the \Rfunction{GRanges}
 
80
constructor function. For example,
 
81
 
 
82
<<example-GRanges>>=
 
83
gr <-
 
84
  GRanges(seqnames =
 
85
          Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
 
86
          ranges =
 
87
          IRanges(1:10, end = 7:16, names = head(letters, 10)),
 
88
          strand =
 
89
          Rle(strand(c("-", "+", "*", "+", "-")),
 
90
              c(1, 2, 2, 3, 2)),
 
91
          score = 1:10,
 
92
          GC = seq(1, 0, length=10))
 
93
gr
 
94
@
 
95
 
 
96
\noindent
 
97
creates a \Rclass{GRanges} object with 10 single interval features.
 
98
The output of the \Rclass{GRanges} \Rmethod{show} method separates the
 
99
information into a left and right hand region that are separated by
 
100
\Rcode{|} symbols. The genomic coordinates (seqnames, ranges, and strand)
 
101
are located on the left-hand side and the metadata columns (annotation)
 
102
are located on the right. For this example, the metadata is
 
103
comprised of \Rcode{score} and \Rcode{GC} information, but almost
 
104
anything can be stored in the metadata portion of a \Rclass{GRanges}
 
105
object.
 
106
 
 
107
The components of the genomic coordinates within a \Rclass{GRanges}
 
108
object can be extracted using the \Rmethod{seqnames}, \Rmethod{ranges},
 
109
and \Rmethod{strand} accessor functions.
 
110
 
 
111
<<GRanges-location-accessors>>=
 
112
seqnames(gr)
 
113
ranges(gr)
 
114
strand(gr)
 
115
@
 
116
 
 
117
Stored annotations for these coordinates can be extracted as a
 
118
\Rclass{DataFrame} object using the \Rmethod{mcols} accessor.
 
119
 
 
120
<<metadataAccess>>=
 
121
mcols(gr)
 
122
mcols(gr)$score
 
123
@
 
124
 
 
125
Finally, the total lengths of the various sequences that the ranges
 
126
are aligned to can also be stored in the \Rclass{GRanges} object. So
 
127
if this is data from \textit{Homo sapiens}, we can set the values as:
 
128
 
 
129
<<setSeqLengths>>=
 
130
seqlengths(gr) <- c(249250621,243199373,198022430)
 
131
@
 
132
 
 
133
And then retrieves as:
 
134
<<setSeqLengths2>>=
 
135
seqlengths(gr)
 
136
@
 
137
 
 
138
Methods for accessing the \Rmethod{length} and \Rmethod{names} have
 
139
also been defined.
 
140
 
 
141
<<names>>=
 
142
names(gr)
 
143
length(gr)
 
144
@
 
145
 
 
146
 
 
147
\subsection{Splitting and combining \Rclass{GRanges} objects}
 
148
 
 
149
\Rclass{GRanges} objects can be devided into groups using the
 
150
\Rmethod{split} method. This produces a \Rclass{GRangesList} object,
 
151
a class that will be discussed in detail in the next section.
 
152
 
 
153
<<splitAppendGRanges>>=
 
154
sp <- split(gr, rep(1:2, each=5))
 
155
sp
 
156
@
 
157
 
 
158
If you then grab the components of this list, they can also be merged
 
159
by using the \Rmethod{c} and \Rmethod{append} methods.
 
160
 
 
161
<<combine>>=
 
162
c(sp[[1]], sp[[2]])
 
163
@
 
164
 
 
165
\subsection{Subsetting  \Rclass{GRanges} objects}
 
166
 
 
167
The expected subsetting operations are also available for
 
168
\Rclass{GRanges} objects.
 
169
 
 
170
<<subset1>>=
 
171
gr[2:3]
 
172
@
 
173
 
 
174
A second argument to the \Rmethod{[} subset operator can be used
 
175
to specify which metadata columns to extract from the
 
176
\Rclass{GRanges} object. For example,
 
177
 
 
178
<<subset2>>=
 
179
gr[2:3, "GC"]
 
180
@
 
181
 
 
182
You can also assign into elements of the \Rclass{GRanges} object.
 
183
Here is an example where the 2nd row of a \Rclass{GRanges} object
 
184
is replaced with the 1st row of \Robject{gr}.
 
185
 
 
186
<<assign1>>=
 
187
singles <- split(gr, names(gr))
 
188
grMod <- gr
 
189
grMod[2] <- singles[[1]]
 
190
head(grMod, n=3)
 
191
@
 
192
 
 
193
Here is a second example where the metadata for score from the 3rd
 
194
element is replaced with the score from the 2nd row etc.
 
195
 
 
196
<<assign2>>=
 
197
grMod[2,1] <- singles[[3]][,1]
 
198
head(grMod, n=3)
 
199
@
 
200
 
 
201
There are also methods to repeat, reverse, or select specific portions
 
202
of \Rclass{GRanges} objects.
 
203
 
 
204
<<other>>=
 
205
rep(singles[[2]], times = 3)
 
206
rev(gr)
 
207
head(gr,n=2)
 
208
tail(gr,n=2)
 
209
window(gr, start=2,end=4)
 
210
gr[IRanges(start=c(2,7), end=c(3,9))]
 
211
@
 
212
 
 
213
 
 
214
\subsection{Basic interval operations for \Rclass{GRanges} objects}
 
215
 
 
216
Basic interval characteristics of \Rclass{GRanges} objects can
 
217
be extracted using the \Rmethod{start}, \Rmethod{end}, \Rmethod{width},
 
218
and \Rmethod{range} methods.
 
219
 
 
220
<<IRangesStuff>>=
 
221
g <- gr[1:3]
 
222
g <- append(g, singles[[10]])
 
223
start(g)
 
224
end(g)
 
225
width(g)
 
226
range(g)
 
227
@
 
228
 
 
229
The \Rclass{GRanges} class also has many methods for manipulating the
 
230
intervals. For example, the \Rmethod{flank} method can be used to recover
 
231
regions flanking the set of ranges represented by the \Rclass{GRanges}
 
232
object. So to get a \Rclass{GRanges} object containing the ranges that
 
233
include the 10 bases upstream of the ranges:
 
234
 
 
235
<<flank>>=
 
236
flank(g, 10)
 
237
@
 
238
 
 
239
And to include the downstream bases:
 
240
 
 
241
<<flank2>>=
 
242
flank(g, 10, start=FALSE)
 
243
@
 
244
 
 
245
Similar to \Rmethod{flank}, there are also operations to
 
246
\Rmethod{resize} and \Rmethod{shift} our \Rclass{GRanges} object. The
 
247
\Rmethod{shift} method will move the ranges by a specific number of
 
248
base pairs, and the \Rmethod{resize} method will extend the ranges by
 
249
a specified width.
 
250
 
 
251
<<shiftAndResize>>=
 
252
shift(g, 5)
 
253
resize(g, 30)
 
254
@
 
255
 
 
256
The \Rmethod{reduce} will align the ranges and merge overlapping
 
257
ranges to produce a simplified set.
 
258
 
 
259
<<reduce>>=
 
260
reduce(g)
 
261
@
 
262
 
 
263
Sometimes you may be interested in the spaces or the qualities of the
 
264
spaces between the ranges represented by your \Rclass{GRanges} object.
 
265
The \Rmethod{gaps} method will help you calculate the spaces between a
 
266
reduced version of your ranges:
 
267
 
 
268
<<gaps>>=
 
269
gaps(g)
 
270
@
 
271
 
 
272
And sometimes you also may want to know how many quantitatively unique
 
273
fragments your ranges could possibly represent. For this task there
 
274
is the \Rmethod{disjoin} method.
 
275
 
 
276
<<disjoin>>=
 
277
disjoin(g)
 
278
@
 
279
 
 
280
One of the most powerful methods for looking at \Rclass{GRanges}
 
281
objects is the \Rmethod{coverage} method. The \Rmethod{coverage}
 
282
method quantifies the degree of overlap for all the ranges in a
 
283
\Rclass{GRanges} object.
 
284
 
 
285
<<coverage>>=
 
286
coverage(g)
 
287
@
 
288
 
 
289
 
 
290
\subsection{Interval set operations for \Rclass{GRanges} objects}
 
291
 
 
292
There are also operations for calculating relationships between
 
293
different \Rclass{GRanges} objects. Here are a some examples for
 
294
how you can calculate the \Rmethod{union}, the \Rmethod{intersect}
 
295
and the asymmetric difference (using \Rmethod{setdiff}).
 
296
 
 
297
<<intervals1>>=
 
298
g2 <- head(gr, n=2)
 
299
union(g, g2)
 
300
intersect(g, g2)
 
301
setdiff(g, g2)
 
302
@
 
303
 
 
304
In addition, there is similar set of operations that act at the level
 
305
of the individual ranges within each \Rclass{GRanges}. These
 
306
operations all begin with a ``p", which is short for parallel. A
 
307
requirement for this set of operations is that the number of elements
 
308
in each \Rclass{GRanges} object has to be the same, and that both of
 
309
the objects have to have the same seqnames and strand assignments
 
310
throughout.
 
311
 
 
312
<<intervals2>>=
 
313
g3 <- g[1:2]
 
314
ranges(g3[1]) <- IRanges(start=5, end=12)
 
315
punion(g2, g3)
 
316
pintersect(g2, g3)
 
317
psetdiff(g2, g3)
 
318
@
 
319
 
 
320
For even more information on the \Rcode{GRanges} classes be sure to
 
321
consult the manual page.
 
322
 
 
323
<<manPage, eval=FALSE>>=
 
324
?GRanges
 
325
@
 
326
 
 
327
 
 
328
\section{\Rclass{GRangesList}: Multiple Interval Range Features}
 
329
 
 
330
Some important genomic features, such as spliced transcripts that are
 
331
are comprised of exons, are inherently compound structures. Such a
 
332
feature makes much more sense when expressed as a compound object
 
333
such as a \Rclass{GRangesList}. Whenever genomic features consist of
 
334
multiple ranges that are grouped by a parent feature, they can be
 
335
represented as \Rclass{GRangesList} object. Consider the simple
 
336
example of the two transcript \Rfunction{GRangesList} below created
 
337
using the \Rfunction{GRangesList} constructor.
 
338
 
 
339
<<example-GRangesList>>=
 
340
gr1 <-
 
341
  GRanges(seqnames = "chr2", ranges = IRanges(3, 6),
 
342
          strand = "+", score = 5L, GC = 0.45)
 
343
gr2 <-
 
344
  GRanges(seqnames = c("chr1", "chr1"),
 
345
          ranges = IRanges(c(7,13), width = 3),
 
346
          strand = c("+", "-"), score = 3:4, GC = c(0.3, 0.5))
 
347
grl <- GRangesList("txA" = gr1, "txB" = gr2)
 
348
grl
 
349
@
 
350
 
 
351
The \Rmethod{show} method for a \Rclass{GRangesList} object displays
 
352
it as a named list of \Rclass{GRanges} objects, where the names of
 
353
this list are considered to be the names of the grouping feature. In
 
354
the example above, the groups of individual exon ranges are represented
 
355
as separate \Rclass{GRanges} objects which are further organized into a
 
356
list structure where each element name is a transcript name. Many
 
357
other combinations of grouped and labeled \Rclass{GRanges} objects are
 
358
possible of course, but this example is expected to be a common
 
359
arrangement. 
 
360
 
 
361
 
 
362
\subsection{Basic \Rclass{GRangesList} accessors}
 
363
 
 
364
Just as with \Rclass{GRanges} object, the components of the genomic
 
365
coordinates within a \Rclass{GRangesList} object can be extracted
 
366
using simple accessor methods. Not surprisingly, the
 
367
\Rclass{GRangesList} objects have many of the same accessors as
 
368
\Rclass{GRanges} objects.  The difference is that many of these
 
369
methods return a list since the input is now essentially a list of
 
370
\Rclass{GRanges} objects. Here are a few examples:
 
371
 
 
372
<<basicGRLAccessors>>=
 
373
seqnames(grl)
 
374
ranges(grl)
 
375
strand(grl)
 
376
@
 
377
 
 
378
The \Rmethod{length} and \Rmethod{names} methods will return the length
 
379
or names of the list and the \Rmethod{seqlengths} method will return the
 
380
set of sequence lengths.
 
381
 
 
382
<<exceptions>>=
 
383
length(grl)
 
384
names(grl)
 
385
seqlengths(grl)
 
386
@
 
387
 
 
388
The \Rmethod{elementLengths} method returns a list of integers
 
389
corresponding to the result of calling \Rmethod{length} on each
 
390
individual \Rclass{GRanges} object contained by the
 
391
\Rclass{GRangesList}. This is a faster alternative to calling
 
392
\Rmethod{lapply} on the \Rclass{GRangesList}.
 
393
 
 
394
<<elementLengths>>=
 
395
elementLengths(grl)
 
396
@
 
397
 
 
398
You can also use \Rmethod{isEmpty} to test if a \Rclass{GRangesList}
 
399
object contains anything.
 
400
 
 
401
<<isEmpty>>=
 
402
isEmpty(grl)
 
403
@
 
404
 
 
405
Finally, in the context of a \Rclass{GRangesList} object, the
 
406
\Rmethod{mcols} method performs a similar operation to what it
 
407
does on a \Rclass{GRanges} object. However, this metadata now
 
408
refers to information at the list level instead of the level
 
409
of the individual \Rclass{GRanges} objects.
 
410
 
 
411
<<mcolsGRL>>=
 
412
mcols(grl) <- c("Transcript A","Transcript B")
 
413
mcols(grl)
 
414
@
 
415
 
 
416
 
 
417
\subsection{Combining \Rclass{GRangesList} objects}
 
418
 
 
419
\Rclass{GRangesList} objects can be unlisted to combine the separate
 
420
\Rclass{GRanges} objects that they contain as an expanded
 
421
\Rclass{GRanges}.
 
422
 
 
423
<<unlistGRL>>=
 
424
ul <- unlist(grl)
 
425
@
 
426
 
 
427
You can also append values together useing \Rmethod{append} or
 
428
\Rmethod{c}.
 
429
 
 
430
% <<appendGRL>>=
 
431
% grl2 <- shift(grl,10)
 
432
% names(grl2) <- c("shiftTxA","shiftTxB")
 
433
 
 
434
% append(grl, grl2)
 
435
% c(grl, grl2)
 
436
% @
 
437
 
 
438
 
 
439
\subsection{Basic interval operations for \Rclass{GRangesList} objects}
 
440
 
 
441
For interval operations, many of the same methods exist for
 
442
\Rclass{GRangesList} objects that exist for \Rclass{GRanges} objects.
 
443
 
 
444
<<intOpsGRL>>=
 
445
start(grl)
 
446
end(grl)
 
447
width(grl)
 
448
@
 
449
 
 
450
And as with \Rclass{GRanges} objects, you can also shift all the
 
451
\Rclass{GRanges} objects in a \Rclass{GRangesList} object, or
 
452
calculate the coverage. Both of these operations are also carried out
 
453
across each \Rclass{GRanges} list member.
 
454
 
 
455
<<coverageGRL>>=
 
456
shift(grl, 20)
 
457
coverage(grl)
 
458
@
 
459
 
 
460
 
 
461
\subsection{Subsetting \Rclass{GRangesList} objects}
 
462
 
 
463
As you might guess, the subsetting of a \Rclass{GRangesList} object is
 
464
quite different from subsetting on a \Rclass{GRanges} object in that
 
465
it acts as if you are subseting a list. If you try out the following
 
466
you will notice that the standard conventions have been followed.
 
467
 
 
468
<<subsetGRL, eval=FALSE>>=
 
469
grl[1]
 
470
grl[[1]]
 
471
grl["txA"]
 
472
grl$txB
 
473
@
 
474
 
 
475
But in addition to this, when subsetting a \Rclass{GRangesList}, you
 
476
can also pass in a second parameter (as with a \Rclass{GRanges} object)
 
477
to again specify which of the metadata columns you wish to select.
 
478
 
 
479
<<subsetGRL2>>=
 
480
grl[1, "score"]
 
481
grl["txB", "GC"]
 
482
@ 
 
483
 
 
484
The \Rmethod{head}, \Rmethod{tail}, \Rmethod{rep}, \Rmethod{rev},
 
485
and \Rmethod{window} methods all behave as you would expect them
 
486
to for a list object. For example, the elements
 
487
referred to by \Rmethod{window} are now list
 
488
elements instead of \Rclass{GRanges} elements.
 
489
 
 
490
<<otherSubsetGRL>>=
 
491
rep(grl[[1]], times = 3)
 
492
rev(grl)
 
493
head(grl, n=1)
 
494
tail(grl, n=1)
 
495
window(grl, start=1, end=1)
 
496
grl[IRanges(start=2, end=2)]
 
497
@
 
498
 
 
499
 
 
500
\subsection{Looping over \Rclass{GRangesList} objects}
 
501
 
 
502
For \Rclass{GRangesList} objects there is also a family of
 
503
\Rmethod{apply} methods. These include \Rmethod{lapply},
 
504
\Rmethod{sapply}, \Rmethod{mapply}, \Rmethod{endoapply},
 
505
\Rmethod{mendoapply}, \Rmethod{Map}, and \Rmethod{Reduce}.
 
506
 
 
507
The different looping methods defined for \Rclass{GRangesList} objects
 
508
are useful for returning different kinds of results.  The standard
 
509
\Rmethod{lapply} and \Rmethod{sapply} behave according to convention,
 
510
with the \Rmethod{lapply} method returning a list and \Rmethod{sapply}
 
511
returning a more simplified output.
 
512
 
 
513
<<lapply>>=
 
514
lapply(grl, length)
 
515
sapply(grl, length)
 
516
@
 
517
 
 
518
As with \Rclass{IRanges} objects, there is also a multivariate version
 
519
of \Rmethod{sapply}, called \Rmethod{mapply}, defined for
 
520
\Rclass{GRangesList} objects. And, if you don't want the results
 
521
simplified, you can call the \Rmethod{Map} method, which does the same
 
522
things as \Rmethod{mapply} but without simplifying the output.
 
523
 
 
524
<<mapply>>=
 
525
grl2 <- shift(grl, 10)
 
526
names(grl2) <- c("shiftTxA", "shiftTxB")
 
527
 
 
528
mapply(c, grl, grl2)
 
529
Map(c, grl, grl2)
 
530
@
 
531
 
 
532
Sometimes, you may not want to get back a simplified output or a list.
 
533
Sometimes you will want to get back a modified version of the
 
534
\Rclass{GRangesList} that you originally passed in. This is
 
535
conceptually similar to the mathematical notion of an endomorphism.
 
536
This is achieved using the \Rmethod{endoapply} method, which will return
 
537
the results as a \Rclass{GRangesList} object.
 
538
 
 
539
<<endoapply>>=
 
540
endoapply(grl,rev)
 
541
@
 
542
 
 
543
And, there is also a multivariate version of the \Rmethod{endoapply}
 
544
method in the form of the \Rmethod{mendoapply} method.
 
545
 
 
546
<<mendoapply>>=
 
547
mendoapply(c,grl,grl2)
 
548
@
 
549
 
 
550
Finally, the \Rmethod{Reduce} method will allow the \Rclass{GRanges}
 
551
objects to be collapsed across the whole of the \Rclass{GRangesList}
 
552
object.
 
553
 
 
554
% Again, this seems like a sub-optimal example to me.
 
555
 
 
556
<<ReduceGRL>>=
 
557
Reduce(c,grl)
 
558
@
 
559
 
 
560
For even more information on the \Rcode{GRangesList} classes be sure to
 
561
consult the manual page.
 
562
 
 
563
<<manPage2, eval=FALSE>>=
 
564
?GRangesList
 
565
@
 
566
 
 
567
 
 
568
\section{Interval overlaps involving \Rclass{GRanges} and \Rclass{GRangesList} objects}
 
569
Interval overlapping is the process of comparing the ranges in two
 
570
objects to determine if and when they overlap. As such, it is perhaps
 
571
the most common operation performed on \Rclass{GRanges} and
 
572
\Rclass{GRangesList} objects. To this end, the \Rpackage{GenomicRanges}
 
573
package provides a family of interval overlap functions. The most general
 
574
of these functions is \Rfunction{findOverlaps}, which takes a query and a
 
575
subject as inputs and returns a \Rclass{Hits} object containing
 
576
the index pairings for the overlapping elements.
 
577
 
 
578
<<findOverlaps>>=
 
579
mtch <- findOverlaps(gr, grl)
 
580
as.matrix(mtch)
 
581
@
 
582
 
 
583
\noindent
 
584
As suggested in the sections discussing the nature of the
 
585
\Rclass{GRanges} and \Rclass{GRangesList} classes, the index in the
 
586
above matrix of hits for a \Rclass{GRanges} object is a single range
 
587
while for a \Rclass{GRangesList} object it is the set of ranges that
 
588
define a "feature".
 
589
 
 
590
Another function in the overlaps family is \Rfunction{countOverlaps},
 
591
which tabulates the number of overlaps for each element in the query.
 
592
 
 
593
<<countOL>>=
 
594
countOverlaps(gr, grl)
 
595
@
 
596
 
 
597
A third function in this family is \Rfunction{subsetByOverlaps},
 
598
which extracts the elements in the query that overlap at least one
 
599
element in the subject.
 
600
 
 
601
<<subsetByOverlaps>>=
 
602
subsetByOverlaps(gr,grl)
 
603
@
 
604
 
 
605
Finally, you can use the \Rcode{select} argument to get the index
 
606
of the first overlapping element in the subject for each element
 
607
in the query.
 
608
 
 
609
<<select-first>>=
 
610
findOverlaps(gr, grl, select="first")
 
611
findOverlaps(grl, gr, select="first")
 
612
@
 
613
 
 
614
 
 
615
\section{Genomic Alignments} 
 
616
 
 
617
In addition to \Rclass{GRanges} and \Rclass{GRangesList} classes, the
 
618
\GenomicRanges{} package defines the \Rclass{GAlignments} class,
 
619
which is a more specialized container for storing a set of genomic
 
620
alignments. The class is intended to support alignments in general,
 
621
not only those coming from a 'Binary Alignment Map' or 'BAM' files.
 
622
Also alignments with gaps in the reference sequence (a.k.a.
 
623
\emph{gapped alignments}) are supported which, for example, makes
 
624
the class suited for storing junction reads from an RNA-seq experiment.
 
625
 
 
626
More precisely, a \Rclass{GAlignments} object is a vector-like
 
627
object where each element describes an \emph{alignment}, that is,
 
628
how a given sequence (called \emph{query} or \emph{read}, typically
 
629
short) aligns to a reference sequence (typically long).
 
630
 
 
631
As shown later in this document, a \Rclass{GAlignments} object
 
632
can be created from a 'BAM' file. In that case, each element in the
 
633
resulting object will correspond to a record in the file.
 
634
One important thing to note though is that not all the information
 
635
present in the BAM/SAM records is stored in the object. In particular,
 
636
for now, we discard the query sequences (SEQ field), the query ids
 
637
(QNAME field), the query qualities (QUAL), the mapping qualities (MAPQ)
 
638
and any other information that is not needed in order to support the
 
639
basic set of operations described in this document.
 
640
This also means that multi-reads (i.e. reads with multiple hits in the
 
641
reference) don't receive any special treatment i.e. the various SAM/BAM
 
642
records corresponding to a multi-read will show up in the \Rclass{GAlignments}
 
643
object as if they were coming from different/unrelated queries.
 
644
Also paired-end reads will be treated as single-end reads and the
 
645
pairing information will be lost. This might change in the future.
 
646
 
 
647
 
 
648
\subsection{Load a `BAM' file into a \Rclass{GAlignments} object}
 
649
 
 
650
First we use the \Rfunction{readGAlignments} function from the
 
651
\Rpackage{GenomicAlignments} package to load a toy `BAM' file into
 
652
a \Rclass{GAlignments} object:
 
653
<<readGAlignments>>=
 
654
library(GenomicAlignments)
 
655
aln1_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
 
656
aln1 <- readGAlignments(aln1_file)
 
657
aln1
 
658
length(aln1)
 
659
@
 
660
 
 
661
3271 `BAM' records were loaded into the object.
 
662
 
 
663
Note that \Rfunction{readGAlignments} would have discarded
 
664
any `BAM' record describing an unaligned query (see description
 
665
of the <flag> field in the SAM Format Specification
 
666
\footnote{\url{http://samtools.sourceforge.net/SAM1.pdf}}
 
667
for more information).
 
668
The reader interested in tracking down these events can always
 
669
use the \Rfunction{scanBam} function but this goes beyond the scope
 
670
of this document.
 
671
 
 
672
\subsection{Simple accessor methods}
 
673
 
 
674
There is one accessor per field displayed by the \Rmethod{show} method
 
675
and it has the same name as the field. All of them return a vector or
 
676
factor of the same length as the object:
 
677
<<accessors>>=
 
678
head(seqnames(aln1))
 
679
seqlevels(aln1)
 
680
head(strand(aln1))
 
681
head(cigar(aln1))
 
682
head(qwidth(aln1))
 
683
head(start(aln1))
 
684
head(end(aln1))
 
685
head(width(aln1))
 
686
head(ngap(aln1))
 
687
@
 
688
 
 
689
\subsection{More accessor methods}
 
690
 
 
691
 
 
692
 
 
693
\end{document}