Basic counting

Let’s count how many reads map to the promoter regions of our genes. Using the bamCount() function, this is straightforward.

proms <- GenomicRanges::promoters(genes, upstream=100, downstream=100)
counts <- bamCount(bampath, proms, verbose=FALSE)
str(counts)

##  int [1:49] 141 68 87 199 97 62 149 76 202 88 ...

The object counts is a vector of the same length as the number of ranges that we are analyzing, the i-th count corresponds to the i-th range.

Accounting for fragment length

With the bamCount() function a read is counted in a range if the 5’ end of the read falls in that range. This might be appropriate when analyzing DNase I hypersensitivity tags, however for ChIP-seq data the immunoprecipitated protein is normally located downstream with respect to the 5’ end of the sequenced reads. To correct for that, it is possible to count reads with a strand-specific shift, i.e. reads will be counted in a range if the shifted 5’ end falls in that range. Note that this shift will move reads mapped to the positive strand to the right and reads mapped to the negative strand to the left with respect to the reference orientation. The shift should correspond approximately to half of the average length of the fragments in the sequencing experiment.

counts <- bamCount(bampath, proms, verbose=FALSE, shift=75)
str(counts)

##  int [1:49] 126 81 91 192 105 74 143 79 161 107 ...

Counting on each strand separately

Sometimes it is necessary to consider the two genomic strands separately. This is achieved with the ss option (separate strands, or strand-specific), and depends also on the strand of the GenomicRanges object.

strand(proms)

## factor-Rle of length 49 with 33 runs
##   Lengths: 1 1 1 1 1 1 1 1 2 1 5 1 1 4 1 ... 2 1 1 1 3 1 1 1 2 3 1 1 1 1 2
##   Values : - + - + - + - + - + - + - + - ... - + - + - + - + - + - + - + -
## Levels(3): + - *

counts <- bamCount(bampath, proms, verbose=FALSE, ss=TRUE)
str(counts)

##  int [1:2, 1:49] 87 54 50 18 48 39 162 37 52 45 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:2] "sense" "antisense"
##   ..$ : NULL

Now counts is a matrix with two rows, one for the sense strand, the other for the antisense strand. Note that the sense of a read is decided also by the region it falls into, so if both the region and the read are on the same strand the read is counted as a sense read, otherwise as an antisense read.

Counting on each strand separately

As for the bamCount() function, also with bamProfile() the reads can be counted for each strand separately

sssigs <- bamProfile(bampath, genes, verbose=FALSE, ss=TRUE)
sssigs

## CountSignals object with 49 strand-specific signals
## [1] signal of width 568
## sense      0 0 0 0 0 0 0 0 0 1 ...
## antisense  0 0 0 0 0 0 0 0 0 0 ...
## [2] signal of width 595
## sense      0 0 0 0 0 0 0 0 0 0 ...
## antisense  0 1 1 0 0 0 0 0 1 0 ...
## [3] signal of width 607
## sense      0 0 1 0 0 0 0 0 1 0 ...
## antisense  0 0 0 0 0 0 0 0 0 1 ...
## [4] signal of width 653
## sense      1 0 1 1 4 1 0 3 1 1 ...
## antisense  1 0 0 0 1 0 0 2 1 1 ...
## [5] signal of width 602
## sense      0 0 0 0 1 0 0 2 0 0 ...
## antisense  0 0 1 0 1 0 0 0 0 0 ...
## ....

Now each signal is a matrix with two rows.

str(sssigs[1])

##  int [1:2, 1:568] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:2] "sense" "antisense"
##   ..$ : NULL

#summing up the counts from the two strands is the same as using ss=FALSE
stopifnot(colSums(sssigs[1])==sigs[1])
#the width function takes into account that now the signals are strand-specific
stopifnot(width(sssigs)==width(sigs))
#the length function does not, a strand-specific signal is twice as long
stopifnot(length(sssigs[1])==2*length(sigs[1]))

Let’s summarize this with a plot

xlab <- "offset from start of the region"
ylab <- "counts per base pair (negative means antisense)"
main <- paste0("read profile of the region ", 
    seqnames(genes)[1], ":", start(genes)[1], "-", end(genes)[1])
plot(sigs[1], ylim=c(-max(sigs[1]), max(sigs[1])), ylab=ylab, xlab=xlab, 
    main=main, type="l")
lines(sssigs[1]["sense",], col="blue")
lines(-sssigs[1]["antisense",], col="red")
legend("topright", c("sense", "antisense", "both"), 
    col=c("blue", "red", "black"), lty=1)

Regions of the same width

In case our ranges have all the same width, a CountSignals object can be immediately converted into a matrix, or an array, with the alignSignals function

#The promoter regions have all the same width
sigs <- bamProfile(bampath, proms, ss=FALSE, verbose=FALSE)
sssigs <- bamProfile(bampath, proms, ss=TRUE, verbose=FALSE)

sigsMat <- alignSignals(sigs)
sigsArr <- alignSignals(sssigs)

The last dimension of the resulting array (or matrix) represents the different ranges, the second-last one represents the base pairs in each region, and in the strand-specific case, the first-one represents the strand of the signal. This can be changed by using the t() function (for matrices) or aperm() (for arrays).

#the dimensions are [base pair, region]
str(sigsMat)

##  int [1:200, 1:49] 0 0 1 0 0 0 0 0 0 1 ...

#the dimensions are [strand, base pair, region]
str(sigsArr)

##  int [1:2, 1:200, 1:49] 0 0 0 0 1 0 0 0 0 0 ...
##  - attr(*, "dimnames")=List of 3
##   ..$ : chr [1:2] "sense" "antisense"
##   ..$ : NULL
##   ..$ : NULL

stopifnot(all(sigsMat == sigsArr["sense",,] + sigsArr["antisense",,]))

Computing the average read profile at promoters in proms is now straightforward

avgSig <- rowMeans(sigsMat)
avgSenseSig <- rowMeans(sigsArr["sense",,])
avgAntisenseSig <- rowMeans(sigsArr["antisense",,])
ylab <- "average counts per base pair"
xlab <- "distance from TSS"
main <- paste0("average profile of ", length(proms), " promoters")
xs <- -99:100
plot(xs, avgSig, ylim=c(0, max(avgSig)), xlab=xlab, ylab=ylab, main=main,
    type="l")
lines(xs, avgSenseSig, col="blue")
lines(xs, avgAntisenseSig, col="red")
legend("topright", c("sense", "antisense", "both"), 
    col=c("blue", "red", "black"), lty=1)

Binning counts

Very often it is better to count reads mapping to small regions instead of single base pairs. Bins are small non-overlapping regions of fixed size tiling a larger region. Instead of splitting your regions of interest into bins, it is easier and much more efficient to provide the binsize option to bamProfile().

binsize <- 20
binnedSigs <- bamProfile(bampath, proms, binsize=binsize)

## Processing /tmp/RtmpGYYvnF/Rinst1c86757385c3/bamsignals/extdata/randomBam.bam: I gEt gOoSeBuMPs wHen I seE yOu pilEuPpiNg...

stopifnot(all(width(binnedSigs)==ceiling(width(sigs)/binsize)))
binnedSigs

## CountSignals object with 49 signals
## [1] signal of width 10
## 4 13 22 8 16 3 10 18 21 26
## [2] signal of width 10
## 6 2 2 9 5 10 8 16 8 2
## [3] signal of width 10
## 9 14 12 6 6 6 3 7 12 12
## [4] signal of width 10
## 7 3 4 5 21 30 37 37 28 27
## [5] signal of width 10
## 7 9 16 8 13 5 7 10 9 13
## ....

In case the ranges’ widths are not multiples of the bin size, a warning will be issued and the last bin in those ranges will be smaller than the others (where “last” depends on the orientation of the region).

Binning means considering a signal at a lower resolution.

avgBinnedSig <- rowMeans(alignSignals(binnedSigs))
#the counts in the bin are the sum of the counts in each base pair
stopifnot(all.equal(colSums(matrix(avgSig, nrow=binsize)),avgBinnedSig))
#let's plot it
ylab <- "average counts per base pair"
plot(xs, avgSig, xlab=xlab, ylab=ylab, main=main, type="l")
lines(xs, rep(avgBinnedSig, each=binsize)/binsize, lty=2)
legend("topright", c("base pair count", "bin count"), lty=c(1, 2))