pileup {Rsamtools} | R Documentation |
Use filters and output formats to calculate pile-up statistics for a BAM file.
Description
pileup
uses PileupParam
and ScanBamParam
objects
to calculate pileup statistics for a BAM file. The result is a
data.frame
with columns summarizing counts of reads overlapping
each genomic position, optionally differentiated on nucleotide,
strand, and position within read.
Usage
pileup(file, index=file, ..., scanBamParam=ScanBamParam(),
pileupParam=PileupParam())
## PileupParam constructor
PileupParam(max_depth=250, min_base_quality=13, min_mapq=0,
min_nucleotide_depth=1, min_minor_allele_depth=0,
distinguish_strands=TRUE, distinguish_nucleotides=TRUE,
ignore_query_Ns=TRUE, include_deletions=TRUE, include_insertions=FALSE,
left_bins=NULL, query_bins=NULL, cycle_bins=NULL)
phred2ASCIIOffset(phred=integer(),
scheme= c("Illumina 1.8+", "Sanger", "Solexa", "Illumina 1.3+",
"Illumina 1.5+"))
Arguments
file |
character(1) or |
index |
When |
... |
Additional arguments, perhaps used by methods. |
scanBamParam |
An instance of |
pileupParam |
An instance of |
max_depth |
integer(1); maximum number of overlapping alignments considered for each position in the pileup. |
min_base_quality |
integer(1); minimum ‘QUAL’ value for
each nucleotide in an alignment. Use |
min_mapq |
integer(1); minimum ‘MAPQ’ value for an alignment to be included in pileup. |
min_nucleotide_depth |
integer(1); minimum count of each nucleotide (independent of other nucleotides) at a given position required for said nucleotide to appear in the result. |
min_minor_allele_depth |
integer(1); minimum count of all nucleotides other than the major allele at a given position required for a particular nucleotide to appear in the result. |
distinguish_strands |
logical(1); |
distinguish_nucleotides |
logical(1); |
ignore_query_Ns |
logical(1); |
include_deletions |
logical(1); |
include_insertions |
logical(1); |
left_bins |
numeric; all same sign; unique positions within a
read to delimit bins. Position within read is based on counting from
the 5' end regardless of strand. Minimum of two values are
required so at least one range can be formed. If you want to delimit bins based on sequencing cycles to, e.g.,
discard later cycles, |
query_bins |
numeric; all same sign; unique positions within a
read to delimit bins. Position within a read is based on counting
from the 5' end when the read aligns to plus strand,
counting from the 3' end when read aligns to minus
strand. Minimum of two values are required so at least one range can
be formed. |
phred |
Either a numeric() (coerced to integer()) PHRED score
(e.g., in the range (0, 41) for the ‘Illumina 1.8+’ scheme)
or character() of printable ASCII characters. When |
scheme |
Encoding scheme, used to translate numeric() PHRED
scores to their ASCII code. Ignored when |
cycle_bins |
DEPRECATED. See |
Details
NB: Use of pileup
assumes familiarity with
ScanBamParam
, and use of left_bins
and
query_bins
assumes familiarity with cut
.
pileup
visits each position in the BAM file, subject to
constraints implied by which
and flag
of
scanBamParam
. For a given position, all reads overlapping the
position that are consistent with constraints dictated by flag
of scanBamParam
and pileupParam
are used for
counting. pileup
also automatically excludes reads with flags
that indicate any of:
unmapped alignment (
isUnmappedQuery
)secondary alignment (
isSecondaryAlignment
)not passing quality controls (
isNotPassingQualityControls
)PCR or optical duplicate (
isDuplicate
)
If no which
argument is supplied to the ScanBamParam
,
behavior reflects that of scanBam
: the entire file is visited
and counted.
Use a yieldSize
value when creating a BamFile
instance to manage memory consumption when using pileup with large BAM
files. There are some differences in how pileup
behavior is
affected when the yieldSize
value is set on the BAM file. See
more details below.
Many of the parameters of the pileupParam
interact. A simple
illustration is ignore_query_Ns
and
distinguish_nucleotides
, as mentioned in the
ignore_query_Ns
section.
Parameters for pileupParam
belong to two categories: parameters
that affect only the filtering criteria (so-called
‘behavior-only’ policies), and parameters that affect
filtering behavior and the schema of the results
(‘behavior+structure’ policies).
distinguish_nucleotides
and distinguish_strands
when set
to TRUE
each add a column (nucleotide
and strand
,
respectively) to the resulting data.frame
. If both are TRUE,
each combination of nucleotide
and strand
at a given
position is counted separately. Setting only one to TRUE
behaves as expected; for example, if only nucleotide
is set to
TRUE
, each nucleotide at a given position is counted
separately, but the distinction of on which strand the nucleotide
appears is ignored.
ignore_query_Ns
determines how ambiguous nucloetides are
treated. By default, ambiguous nucleotides (typically ‘N’ in
BAM files) in alignments are ignored. If ignore_query_Ns
is
FALSE
, ambiguous nucleotides are included in counts; further,
if ignore_query_Ns
is FALSE
and
distinguish_nucleotides
is TRUE
the ‘N’
nucleotide value appears in the nucleotide column when a base at a
given position is ambiguous.
By default, deletions with respect to the reference genome to which
the reads were aligned are included in the counts in a pileup. If
include_deletions
is TRUE
and
distinguish_nucleotides
is TRUE
, the nucleotide column
uses a ‘-’ character to indicate a deletion at a position.
The ‘=’ nucleotide code in the SEQ
field (to mean
‘identical to reference genome’) is supported by pileup, such
that a match with the reference genome is counted separately in the
results if distinguish_nucleotides
is TRUE
.
CIGAR support: pileup
handles the extended CIGAR format by only
heeding operations that contribute to counts (‘M’, ‘D’,
‘I’). If you are confused about how the different CIGAR
operations interact, see the SAM versions of the BAM files used for
pileup
unit testing for a fairly comprehensive human-readable
treatment.
Deletions and clipping:
The extended CIGAR allows a number of operations conceptually similar to a ‘deletion’ with respect to the reference genome, but offer more specialized meanings than a simple deletion. CIGAR ‘N’ operations (not to be confused with ‘N’ used for non-determinate bases in the
SEQ
field) indicate a large skipped region, ‘S’ a soft clip, and ‘H’ a hard clip. ‘N’, ‘S’, and ‘H’ CIGAR operations are never counted: only counts of true deletions (‘D’ in the CIGAR) can be included by settinginclude_deletions
toTRUE
.Soft clip operations contribute to the relative position of nucleotides within a read, whereas hard clip and refskip operations do not. For example, consider a sequence in a bam file that starts at 11, with a CIGAR
2S1M
andSEQ
fieldttA
. The cycle position of theA
nucleotide will be3
, but the reported position will be11
. If usingleft_bins
orquery_bins
it might make sense to first preprocess your files to remove soft clipping.Insertions and padding:
pileup
can include counts of insertion operations by settinginclude_insertions
toTRUE
. Details about insertions are effectively truncated such that each insertion is reduced to a single ‘+’ nucleotide. Information about e.g. the nucleotide code and base quality within the insertion is not included.Because ‘+’ is used as a nucleotide code to represent insertions in
pileup
, counts of the ‘+’ nucleotide participate in voting formin_nucleotide_depth
andmin_minor_allele_depth
functionality.The position of an insertion is the position that precedes it on the reference sequence. Note: this means if
include_insertions
isTRUE
a single read will contribute two nucleotide codes (+
, plus that of the non-insertion base) at a given position if the non-insertion base passes filter criteria.‘P’ CIGAR (padding) operations never affect counts.
The “bin” arguments query_bins
and left_bins
allow users to differentiate pileup counts based on arbitrary
non-overlapping regions within a read. pileup
relies on
cut
to derive bins, but limits input to numeric values
of the same sign (coerced to integer
s), including
+/-Inf
. If a “bin” argument is set pileup
automatically excludes bases outside the implicit outer range. Here
are some important points regarding bin arguments:
query_bins
vs.left_bins
:BAM files store sequence data from the 5' end to the 3' end (regardless of to which strand the read aligns). That means for reads aligned to the minus strand, cycle position within a read is effectively reversed. Take care to use the appropriate bin argument for your use case.
The most common use of a bin argument is to bin later cycles separately from earlier cycles; this is because accuracy typically degrades in later cycles. For this application,
query_bins
yields the correct behavior because bin positions are relative to cycle order (i.e., sensitive to strand).left_bins
(in contrast withquery_bins
) determines bin positions from the 5' end, regardless of strand.Positive or negative bin values can be used to delmit bins based on the effective “start” or “end” of a read. For
left_bin
the effective start is always the 5' end (left-to-right as appear in the BAM file).For
query_bin
the effective start is the first cycle of the read as it came off the sequencer; that means the 5' end for reads aligned to the plus strand and 3' for reads aligned to the minus strand.From effective start of reads: use positive values,
0
, and(+)Inf
. Becausecut
produces ranges in the form (first,last], ‘0’ should be used to create a bin that includes the first position. To account for variable-length reads in BAM files, use(+)Inf
as the upper bound on a bin that extends to the end of all reads.From effective end of reads: use negative values and
-Inf
.-1
denotes the last position in a read. Becausecut
produces ranges in the form (first,last], specify the lower bound of a bin by using one less than the desired value, e.g., a bin that captures only the second nucleotide from the last position:query_bins=c(-3, -2)
. To account for variable-length reads in BAM files, use-Inf
as the lower bound on a bin that extends to the beginning of all reads.
pileup
obeys yieldSize
on BamFile
objects
with some differences from how scanBam
responds to
yieldSize
. Here are some points to keep in mind when using
pileup
in conjunction with yieldSize
:
BamFile
s with ayieldSize
value set, once opened and used withpileup
, should not be used with other functions that accept aBamFile
instance as input. Create a newBamFile
instance instead of trying to reuse.pileup
only returns genomic positions for which all input has been processed.pileup
will hold in reserve positions for which only partial data has been processed. Positions held in reserve will be returned upon subsequent calls topileup
when all the input for a given position has been processed.The correspondence between yieldSize and the number of rows in the
data.frame
returned frompileup
is not 1-to-1.yieldSize
only limits the number of alignments processed from the BAM file each time the file is read. Only a handful of alignments can lead to many distinct records in the result.Like
scanBam
,pileup
uses an empty return object (a zero-rowdata.frame
) to indicate end-of-input.Sometimes reading
yieldSize
records from the BAM file does not result in any completed positions. In order to avoid returning a zero-rowdata.frame
pileup
will repeatedly processyieldSize
additional records until at least one position can be returned to the user.
Value
For pileup
a data.frame
with 1 row per unique
combination of differentiating column values that satisfied filter
criteria, with frequency (count
) of unique combination. Columns
seqnames
, pos
, and count
always appear; optional
strand
, nulceotide
, and left_bin
/
query_bin
columns appear as dictated by arguments to
PileupParam
.
Column details:
seqnames
: factor. SAM ‘RNAME’ field of alignment.pos
: integer(1). Genomic position of base. Derived by offset from SAM ‘POS’ field of alignment.strand
: factor. ‘strand’ to which read aligns.nucleotide
: factor. ‘nucleotide’ of base, extracted from SAM ‘SEQ’ field of alignment.left_bin
/query_bin
: factor. Bin in which base appears.count
: integer(1). Frequency of combination of differentiating fields, as indicated by values of other columns.
See scanBam
for more detailed explanation of SAM fields.
If a which
argument is specified for the scanBamParam
, a
which_label
column (factor
in the form
‘rname:start-end’) is automatically included. The
which_label
column is used to maintain grouping of results,
such that two queries of the same genomic region can be
differentiated.
Order of rows in data.frame
is first by order of
seqnames
column based on the BAM index file, then
non-decreasing order on columns pos
, then nucleotide
,
then strand
, then left_bin
/ query_bin
.
PileupParam
returns an instance of PileupParam class.
phred2ASCIIOffset
returns a named integer vector of the same
length or number of characters in phred
. The names are the
ASCII encoding, and the values are the offsets to be used with
min_base_quality
.
Author(s)
Nathaniel Hayden nhayden@fredhutch.org
See Also
For the relationship between PHRED scores and ASCII encoding, see https://en.wikipedia.org/wiki/FASTQ_format#Encoding.
Examples
## The examples below apply equally to pileup queries for specific
## genomic ranges (specified by the 'which' parameter of 'ScanBamParam')
## and queries across entire files; the entire-file examples are
## included first to make them easy to find. The more detailed examples
## of pileup use queries with a 'which' argument.
library("RNAseqData.HNRNPC.bam.chr14")
fl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
## Minimum base quality to be tallied
p_param <- PileupParam(min_base_quality = 10L)
## no 'which' argument to ScanBamParam: process entire file at once
res <- pileup(fl, pileupParam = p_param)
head(res)
table(res$strand, res$nucleotide)
## no 'which' argument to ScanBamParam with 'yieldSize' set on BamFile
bf <- open(BamFile(fl, yieldSize=5e4))
repeat {
res <- pileup(bf, pileupParam = p_param)
message(nrow(res), " rows in result data.frame")
if(nrow(res) == 0L)
break
}
close(bf)
## pileup for the first half of chr14 with default Pileup parameters
## 'which' argument: process only specified genomic range(s)
sbp <- ScanBamParam(which=GRanges("chr14", IRanges(1, 53674770)))
res <- pileup(fl, scanBamParam=sbp, pileupParam = p_param)
head(res)
table(res$strand, res$nucleotide)
## specify pileup parameters: include ambiguious nucleotides
## (the 'N' level in the nucleotide column of the data.frame)
p_param <- PileupParam(ignore_query_Ns=FALSE, min_base_quality = 10L)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
head(res)
table(res$strand, res$nucleotide)
## Don't distinguish strand, require a minimum frequency of 7 for a
## nucleotide at a genomic position to be included in results.
p_param <- PileupParam(distinguish_strands=TRUE,
min_base_quality = 10L,
min_nucleotide_depth=7)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
head(res)
table(res$nucleotide, res$strand)
## Any combination of the filtering criteria is possible: let's say we
## want a "coverage pileup" that only counts reads with mapping
## quality >= 13, and bases with quality >= 10 that appear at least 4
## times at each genomic position
p_param <- PileupParam(distinguish_nucleotides=FALSE,
distinguish_strands=FALSE,
min_mapq=13,
min_base_quality=10,
min_nucleotide_depth=4)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
head(res)
res <- res[, c("pos", "count")] ## drop seqnames and which_label cols
plot(count ~ pos, res, pch=".")
## ASCII offsets to help specify min_base_quality, e.g., quality of at
## least 10 on the Illumina 1.3+ scale
phred2ASCIIOffset(10, "Illumina 1.3+")
## Well-supported polymorphic positions (minor allele present in at
## least 5 alignments) with high map quality
p_param <- PileupParam(min_minor_allele_depth=5,
min_mapq=40,
min_base_quality = 10,
distinguish_strand=FALSE)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
dim(res) ## reduced to our biologically interesting positions
head(xtabs(count ~ pos + nucleotide, res))
## query_bins
## basic use of positive bins: single pivot; count bases that appear in
## first 10 cycles of a read separately from the rest
p_param <- PileupParam(query_bins=c(0, 10, Inf), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
## basic use of positive bins: simple range; only include bases in
## cycle positions 6-10 within read
p_param <- PileupParam(query_bins=c(5, 10), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
## basic use of negative bins: single pivot; count bases that appear in
## last 3 cycle positions of a read separately from the rest.
p_param <- PileupParam(query_bins=c(-Inf, -4, -1), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
## basic use of negative bins: drop nucleotides in last two cycle
## positions of each read
p_param <- PileupParam(query_bins=c(-Inf, -3), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
## typical use: beginning, middle, and end segments; because of the
## nature of sequencing technology, it is common for bases in the
## beginning and end segments of each read to be less reliable. pileup
## makes it easy to count each segment separately.
## Assume determined ahead of time that for the data 1-7 makes sense for
## beginning, 8-12 middle, >=13 end (actual cycle positions should be
## tailored to data in actual BAM files).
p_param <- PileupParam(query_bins=c(0, 7, 12, Inf), ## note non-symmetric
min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
xt <- xtabs(count ~ nucleotide + query_bin, res)
print(xt)
t(t(xt) / colSums(xt)) ## cheap normalization for illustrative purposes
## 'implicit outer range': in contrast to c(0, 7, 12, Inf), suppose we
## still want to have beginning, middle, and end segements, but know
## that the first three cycles and any bases beyond the 16th cycle are
## irrelevant. Hence, the implicit outer range is (3,16]; all bases
## outside of that are dropped.
p_param <- PileupParam(query_bins=c(3, 7, 12, 16), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
xt <- xtabs(count ~ nucleotide + query_bin, res)
print(xt)
t(t(xt) / colSums(xt))
## single-width bins: count each cycle within a read separately.
p_param <- PileupParam(query_bins=seq(1,20), min_base_quality = 10)
res <- pileup(fl, scanBamParam=sbp, pileupParam=p_param)
xt <- xtabs(count ~ nucleotide + query_bin, res)
print(xt[ , c(1:3, 18:19)]) ## fit on one screen
print(t(t(xt) / colSums(xt))[ , c(1:3, 18:19)])