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

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 setting include_deletions to TRUE.

  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 and SEQ field ttA. The cycle position of the A nucleotide will be 3, but the reported position will be 11. If using left_bins or query_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 setting include_insertions to TRUE. 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 for min_nucleotide_depth and min_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 is TRUE 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 integers), 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 with query_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. Because cut 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. Because cut 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:

• BamFiles with a yieldSize value set, once opened and used with pileup, should not be used with other functions that accept a BamFile instance as input. Create a new BamFile 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 to pileup 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 from pileup 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-row data.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-row data.frame pileup will repeatedly process yieldSize 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

• Rsamtools

• BamFile

• ScanBamParam

• scanBam

• cut

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)])