Install from GitHub

The development version is available at guertinlab/bigWig:

if (!requireNamespace("remotes", quietly = TRUE))
    install.packages("remotes")
remotes::install_github("guertinlab/bigWig")

Install from a local directory

After downloading and unpackaging the source files, use the following commands to install the package.

setwd('bigWig-master')
system('R CMD INSTALL bigWig')

bigWig format

bigWig files are genomic intervals with corresponding intensity values that are stored in an indexed binary format. These files are not readily readable by humans, but the format allows for large continuous data to be stored compactly and accessed quickly.

load.bigWig

load.bigWig(filename, udcDir = NULL)

• arguments

  • filename [required] is a string, which is either the the local file directory or URL.

  • udcDir is a string that specifies the location for storing cached copies of remote files locally while they are in use. These are destroyed when you unload the bigWig. If left as the default udcDir = NULL, then it uses /tmp/udcCache.

• value

  • Returns an object of the bigWig class, which contains a C pointer and the following accessible parts:
    • basesCovered: Number of bp covered by bigWig data.
    • chromSizes: Integer vector with number of bases per chromosome (same order as chroms vector.
    • chroms: Character vector with chromosome names for which there is data in the bigWig file.
    • isCompressed: Logical value indicating state of internal bigWig file compression.
    • isSwapped: Logical value indicating if binary representation is byte swapped (big vs. little endian).
    • max: Numeric value containing the maximum value stored in the bigWig file.
    • mean: Numeric value containing the mean value across all covered bases in the bigWig file.
    • min: Numeric value containing the minimum value stored in the bigWig file.
    • primaryDataSize: Integer number of bytes used by data.
    • primaryIndexSize: Integer number of bytes used by index.
    • std: Standard deviation of bigWig values over covered bases.
    • version: Integer version number for bigWig file format.
    • zoomLevels: Integer number of zoom level supported by bigWig file.

load.bigWig creates a bigWig class object in R. This object contains relevant information about the bigWig file and serves as a pointer to the underlying C object of the entire bigWig file. The only parameter required for this is a string of the location and filename. udcDir is only used if you want to keep the downloaded bigWig file locally if filename is a URL.

#load bigWig into variable bw
setwd('./bigWig')

bw=load.bigWig('../inst/extdata/bp.bigWig')

All of the attributes of the object can be accessed using attributes and each individual can be accessed using $

# list all attributes
attributes(bw)

## $handle_ptr
## <pointer: 0x64bc68974090>
## 
## $names
##  [1] "version"          "isCompressed"     "isSwapped"        "primaryDataSize" 
##  [5] "primaryIndexSize" "zoomLevels"       "chroms"           "chromSizes"      
##  [9] "basesCovered"     "mean"             "min"              "max"             
## [13] "std"             
## 
## $class
## [1] "bigWig"

#access an individual attribute
bw$basesCovered

## [1] 15

print.bigWig

print.bigWig(bw) is used to print all of the attributes contained within the object. The function displays information about the bigWig instance in the same format as the UCSC command line tool bigWigInfo.

print.bigWig(bw)

## bigWig
##  version: 4 
##  isCompressed: yes
##  isSwapped: no
##  primaryDataSize: 90 
##  primaryIndexSize: 6,204 
##  zoomLevels: 2 
##  chromCount: 1 
##      chr1 248956422 
##  basesCovered: 15 
##  mean:  2.333333 
##  min:  1 
##  max:  4 
##  std:  1.290994

• arguments
  • bw is the pointer of the underlying C object created in load.bigWig

unload.bigWig

unload.bigWig(bw) destroys the C object and removes it from memory. Use rm() or remove() to clear the R object.

unload.bigWig(bw)

• arguments
  • bw is the pointer of the underlying C object created in load.bigWig

#destroy C object
unload.bigWig(bw)
ls()

## [1] "bw"

#remove variable in R
remove(bw)
ls()

## character(0)

bigWig file structure

query.bigWig allows you to search the bigWig files using chromosome string (chrom='chr1') and genomic window (start=1, end = 115), both are integers and end is inclusive meaning it searches up to and including end. The query results are printed to the command line. Note how the output of query reflects the original structure of the bigWig file (Figure @ref(fig:bigWigs)). Each row that is output from a query.bigWig call is a genomic interval that is referred to as a probe in the relevant functions.

# load the three bigWigs
bw.bp = load.bigWig('../inst/extdata/bp.bigWig')
bw.probes = load.bigWig('../inst/extdata/probes.bigWig')
bw.splitprobes = load.bigWig('../inst/extdata/splitprobes.bigWig')

#note differences in the bigWig structures
query.bigWig(bw.probes, 'chr1', 100, 115)

##      start end value
## [1,]   100 105     1
## [2,]   105 110     4
## [3,]   110 115     2

query.bigWig(bw.splitprobes, 'chr1', 100, 115)

##      start end value
## [1,]   100 102     1
## [2,]   102 105     1
## [3,]   105 109     4
## [4,]   109 110     4
## [5,]   110 112     2
## [6,]   112 115     2

query.bigWig(bw.bp, chrom='chr1',start=100, end=115)

##       start end value
##  [1,]   100 101     1
##  [2,]   101 102     1
##  [3,]   102 103     1
##  [4,]   103 104     1
##  [5,]   104 105     1
##  [6,]   105 106     4
##  [7,]   106 107     4
##  [8,]   107 108     4
##  [9,]   108 109     4
## [10,]   109 110     4
## [11,]   110 111     2
## [12,]   111 112     2
## [13,]   112 113     2
## [14,]   113 114     2
## [15,]   114 115     2

The default behavior is to clip the bigWig intervals to the queried regions. The bw.probes variable and underlying bigWig structure can be used to highlight the clip= option.

query.bigWig(bw.probes, 'chr1', 104, 111, clip=FALSE)

##      start end value
## [1,]   100 105     1
## [2,]   105 110     4
## [3,]   110 115     2

query.bigWig(bw.probes, 'chr1', 104, 111, clip=TRUE)

##      start end value
## [1,]   104 105     1
## [2,]   105 110     4
## [3,]   110 111     2

The query can be set as a variable for storage.

bw.probes.Q = query.bigWig(bw.probes, 'chr1', 100, 115)

Access the array as an indexed array; the following returns the first row.

bw.probes.Q[1,]

## start   end value 
##   100   105     1

Standard [X,Y] indexing returns the specified row and column.

bw.probes.Q[1,2]

## end 
## 105

The genomic coordinate variable strings are keywords that can be used to access the respective columns.

bw.probes.Q[1,'start']

## start 
##   100

region query

The bp and probe query functions takes a region defined by chrom, start and end and returns the result of the specified operation.

region.bpQuery.bigWig(bw, chrom, start, end, strand = NA
                       op = "sum", abs.value = FALSE,
                       gap.value = 0, bwMap = NULL)

region.probeQuery.bigWig(bw, chrom, start, end,
                       op = "wavg", abs.value = FALSE,
                       gap.value = NA)

• arguments
  • bw is the pointer of the underlying C object created in load.bigWig
  • chrom is a string representing the chromosome to which the query interval belongs
  • start is an integer value defining the start of the query interval
  • end is an integer value defining the end of the query interval
  • strand + or - character indicating the strand of the supplied coordinates (bpQuery only)
  • op is a string representing the operation to perform on the interval.
    • sum adds all the counts
    • avg averages the counts
    • min finds the minimum value
    • max finds the maximum value
    • wavg weighted average of the values–only pertains to probeQuery
  • abs.value is a logical argument which determines if the absolute value of the input is performed before the op.
  • gap.value is an integer value that replaces areas that have no overlaps
  • bwMap a bigWig file of coordinates that cannot be uniquely mapped. Note that the sequence read length of the original FASTQ file should determine the k-mer mappability for this file

All bpQuery functions are insensitive to the structure of the original bigWig file, because each base position is evaluated separately. However, probeQuery functions consider each genomic interval as a separate entity, or probe, and evaluates them separately. The following region.probeQuery.bigWig evaluations highlight the different outputs that result from differentially structured bigWig files that have identical values at each genomic position (see Figure @ref(fig:bigWigs)). Note that the output for each command is the sum of the value column output from the first code chunks in Section @ref(struc).

region.probeQuery.bigWig(bw.probes, 'chr1', 100, 115, op = 'sum')

## [1] 7

region.probeQuery.bigWig(bw.splitprobes, 'chr1', 100, 115, op = 'sum')

## [1] 14

region.probeQuery.bigWig(bw.bp, 'chr1', 100, 115, op = 'sum')

## [1] 35

In contrast, the region.bpQuery.bigWig function considers each base position within each genomic interval input separately. These bigWig files have identical values at each position, so the calculations are identical.

region.bpQuery.bigWig(bw.probes, 'chr1', 100, 115, op = 'sum')

## [1] 35

region.bpQuery.bigWig(bw.splitprobes, 'chr1', 100, 115, op = 'sum')

## [1] 35

region.bpQuery.bigWig(bw.bp, 'chr1', 100, 115, op = 'sum')

## [1] 35

operations (op)

region.bpQuery.bigWig(bw.probes,'chr1', 100, 115, op='max')

## [1] 4

region.probeQuery.bigWig(bw.probes,'chr1', 100, 115, op='max')

## [1] 4

region.bpQuery.bigWig(bw.probes, 'chr1', 100, 115, op='min')

## [1] 1

region.probeQuery.bigWig(bw.probes, 'chr1', 100, 115, op='min')

## [1] 1

region.bpQuery.bigWig(bw.splitprobes,'chr1', 100, 115, op='avg')

## [1] 2.333333

region.probeQuery.bigWig(bw.splitprobes,'chr1', 100, 112, op='avg')

## [1] 2.4

region.bpQuery.bigWig(bw.probes,'chr1', 85, 115, op='avg')

## [1] 1.166667

region.probeQuery.bigWig(bw.probes,'chr1', 85, 115, op='avg')

## [1] 2.333333

region.probeQuery.bigWig(bw.splitprobes,'chr1', 102, 109, op='avg')

## [1] 2.5

region.probeQuery.bigWig(bw.splitprobes,'chr1', 102, 109, op='wavg')

## [1] 2.714286

region.probeQuery.bigWig(bw.probes,'chr1', 100, 111, op='wavg')

## [1] 2.454545

query.bigWig(bw.probes,'chr1', 80, 90)

## NULL

region.probeQuery.bigWig(bw.probes, 'chr1', 80, 90, op = 'sum')

## [1] NA

region.bpQuery.bigWig(bw.probes, 'chr1', 80, 90, op='sum')

## [1] 0

region.probeQuery.bigWig(bw.probes, 'chr1', 80, 90, op = 'sum', gap.value=1)

## [1] 1

region.bpQuery.bigWig(bw.probes, 'chr1', 80, 90, op = 'avg', gap.value=100)

## [1] 100

BED format

A standard three column BED file is a tab delimited file that consists of the chromosome name, the start, and end positions of the genomic interval. A BED6 file contains all of the BED columns plus 3 more: name, score, and strand. Only the strand column is considered for bed6 functions described here. Strand defines whether the BED track interval refers to the + or - stand of DNA. More information can be found on UCSC website. BED files are saved with a .bed extension. The bigWig package operates on bed-formatted files that are loaded as data.frames into R.

Load BED file

Load a BED file is to use R’s read.table function, which converts the tab delimited file into an R data.frame. First read in the file string. The header argument refers whether the columns are named in the first row. Here we set header=FALSE. If there are lines prior ot the coordinate information, use the skip= argument to skip the number of lines that precede the genomic intervals.

bed=read.table('../inst/extdata/testBED1.bed',
             header=FALSE, sep='\t', stringsAsFactors=FALSE)
bed

##     V1  V2  V3
## 1 chr1 101 104
## 2 chr1 105 107
## 3 chr1 107 110
## 4 chr1 112 115

To create a BED6 file, you need to define the following columns: chrom, start, end, name, score and strand. bigWig functions don’t use name and score. In the following example we used place holders ‘na’ for name and ‘1’ for score.

bed6=read.table('../inst/extdata/testBED1_strand.bed',
             header=FALSE, sep='\t', stringsAsFactors=FALSE)
bed6

##     V1  V2  V3 V4 V5 V6
## 1 chr1 101 104 na  1  +
## 2 chr1 101 104 na  1  -
## 3 chr1 105 107 na  1  +
## 4 chr1 107 110 na  1  +
## 5 chr1 112 115 na  1  -

center, fiveprime, and threeprime functions

center.bed, fiveprime.bed, and threeprime.bed take an original BED file and transform each row’s start and end columns. The functions are strand-specific and they differ by the anchor point of the window. For a BED6 file, threeprime.bed and fiveprime.bed refers to upstream and downstream are relative to the strand information. The primary usage of theese transformations is to generate constant windows anchored on a genomic feature, such as a transcription start site, intron/exon junction, or sequence motif.

center.bed(bed, upstreamWindow, downstreamWindow)

fiveprime.bed(bed, upstreamWindow, downstreamWindow)

threeprime.bed(bed, upstreamWindow, downstreamWindow)

• arguments
  • bed is the input BED data.frame.
  • upstreamWindow is an integer number of bases to include upstream of the anchor point.
  • downstreamWindow is an integer number of bases to include downstream of the anchor point.

Anchor Point

The anchor point is different for each function.

center.bed uses the center of the original window. The difference between end and start is taken and divided by 2. If the difference is odd, you are left with a X.5, which is rounded down to X. The anchor point is the start + X.
fiveprime.bed uses the start as the anchor point for BED files and BED6 entires with a + in the sixth strand column. The end is the anchor for BED6 entires with a - in the strand column.
threeprime.bed uses the end as the anchor point for BED files and BED6 entires with a + in the sixth strand column. The start is the anchor for BED6 entires with a - in the strand column.

New Window

The new window is calculated by using the anchor point, upstreamWindow and downstreamWindow.
The new start is anchor point - upstreamWindow.
The new end is the anchor point + 1 + downstreamWindow.   Using the previously loaded bed6 file, we can test a few different scenarios:

bed6

##     V1  V2  V3 V4 V5 V6
## 1 chr1 101 104 na  1  +
## 2 chr1 101 104 na  1  -
## 3 chr1 105 107 na  1  +
## 4 chr1 107 110 na  1  +
## 5 chr1 112 115 na  1  -

Using the center.bed function and upstreamWindow = 0 and downstreamWindow = 0 returns the anchor points.

center.bed(bed6, upstreamWindow = 0, downstreamWindow = 0)

##     V1  V2  V3 V4 V5 V6
## 1 chr1 102 103 na  1  +
## 2 chr1 102 103 na  1  -
## 3 chr1 106 107 na  1  +
## 4 chr1 108 109 na  1  +
## 5 chr1 113 114 na  1  -

Equal and Positive Windows

Often when we query signal bigWig tracks around transcription factor binding sites or sequence motifs, the upstreamWindow and downstreamWindow are equal and positive.

center.bed(bed, upstreamWindow = 5, downstreamWindow = 5)

##     V1  V2  V3
## 1 chr1  97 108
## 2 chr1 101 112
## 3 chr1 103 114
## 4 chr1 108 119

The start values are all anchor point - 5 and the end values are all anchor point + 1 + 5.

fiveprime and threeprime functions

The functions fiveprime.bed and threeprime.bed anchor the window as indicated by the strand. By setting upstreamWindow = 0 and downstreamWindow = 0, you can see that the difference between start and end have no influence on the anchor point. These transformations are most typically used when transforming gene annotation files to query bigWig files relative to stranded genomic features.

Changing the strand affects the anchor point:

• If strand = '+' while using fiveprime.bed
  • anchor point = original start
  • start = anchor point - upstreamWindow
  • end = anchor point + 1 + downstreamWindow
• If strand = '-' while using fiveprime.bed
  • anchor point = original end
  • start = anchor point - downstreamWindow
  • end = anchor point + 1 + upstreamWindow
• If strand = '+' while using threeprime.bed
  • anchor point = original end
  • start = anchor point - 1 - upstreamWindow
  • end = anchor point + downstreamWindow
• If strand = '-' while using threeprime.bed
  • anchor point = original start
  • start = anchor point - downstreamWindow
  • end = anchor point + 1 + upstreamWindow

fiveprime.bed(bed, upstreamWindow = 0, downstreamWindow = 0)

##     V1  V2  V3
## 1 chr1 101 102
## 2 chr1 105 106
## 3 chr1 107 108
## 4 chr1 112 113

threeprime.bed(bed, upstreamWindow = 0, downstreamWindow = 0)

##     V1  V2  V3
## 1 chr1 103 104
## 2 chr1 106 107
## 3 chr1 109 110
## 4 chr1 114 115

fiveprime.bed uses the 5\(^\prime\) end or start as the anchor point, while threeprime.bed uses 3\(^\prime\) or end for the anchor point.

Calculating the new window varies slightly. While fiveprime.bed follows center.bed by

• start = anchor point - upstreamWindow
• end = anchor point + 1 + downstreamWindow

threeprime.bed calculates the window by

• start = anchor point - 1 - upstreamWindow
• end = anchor point + downstreamWindow

Both of these function operate like center.bed other than the initial anchor point.

fiveprime.bed(bed, upstreamWindow = 1, downstreamWindow = 5)

##     V1  V2  V3
## 1 chr1 100 107
## 2 chr1 104 111
## 3 chr1 106 113
## 4 chr1 111 118

threeprime.bed(bed, upstreamWindow = 1, downstreamWindow = 5)

##     V1  V2  V3
## 1 chr1 102 109
## 2 chr1 105 112
## 3 chr1 108 115
## 4 chr1 113 120

If using a BED file without a strand column, fiveprime.bed and threeprime.bed assume that the start is the 5\(^\prime\) end of the sequence.

fiveprime.bed(bed6, upstreamWindow=4, downstreamWindow=2)

##     V1  V2  V3 V4 V5 V6
## 1 chr1  97 104 na  1  +
## 2 chr1 102 109 na  1  -
## 3 chr1 101 108 na  1  +
## 4 chr1 103 110 na  1  +
## 5 chr1 113 120 na  1  -

threeprime.bed(bed6, upstreamWindow=4, downstreamWindow=2)

##     V1  V2  V3 V4 V5 V6
## 1 chr1  99 106 na  1  +
## 2 chr1  99 106 na  1  -
## 3 chr1 102 109 na  1  +
## 4 chr1 105 112 na  1  +
## 5 chr1 110 117 na  1  -

downstream, upstream

These two functions transform the BED data.frame by taking the corresponding anchor point and the window.

downstream.bed(bed, downstreamWindow)

upstream.bed(bed, upstreamWindow)

• arguments
  • bed is the input BED data.frame.
  • upstreamWindow integer number of bases to include upstream of the anchor point.
  • downstreamWindow integer number of bases to include downstream of the anchor point.

downstream.bed uses the original start point [5\(^\prime\)] as the anchor point.

• start = anchor point
• end = anchor point + downstreamWindow

upstream.bed uses the original end point [3\(^\prime\)] as the anchor point.

• start = anchor point - upstreamWindow
• end = anchor point

downstream.bed(bed6,5)

##     V1  V2  V3 V4 V5 V6
## 1 chr1 101 106 na  1  +
## 2 chr1  99 104 na  1  -
## 3 chr1 105 110 na  1  +
## 4 chr1 107 112 na  1  +
## 5 chr1 110 115 na  1  -

upstream.bed(bed6,5)

##     V1  V2  V3 V4 V5 V6
## 1 chr1  96 101 na  1  +
## 2 chr1 104 109 na  1  -
## 3 chr1 100 105 na  1  +
## 4 chr1 102 107 na  1  +
## 5 chr1 115 120 na  1  -

For all transformations, negative numbers are allowed for downstreamWindow and upstreamWindow, but beware that downstream operations will fail if the start coordinate is greater than the end coordinate.

foreach

foreach.bed is a way to quickly apply a function across all rows of a BED file. This function is similar to running apply, but faster since it’s tailored to the structure of a BED file and it avoids a lot of boilerplate code.

foreach.bed(bed, func, envir = parent.frame())

• arguments
  • bed is a data.frame structured like a bed file.`
  • func is the function to apply to each entry in bed. Function must have four arguments: index, chrom, start, end and strand. Index will be a one-based integer corresponding to the current BED line. Chrom is a character string with the chromosome name. Start and end are the coordinates for the current entry (remember that BED files are zero-based left-open intervals). Strand is a character string with the entry’s strand ( ‘+’ or ‘-’) or NA if the bed has less than 6 columns
  • envir is the environment where the function is evaluated. Default value is parent.frame() which corresponds to the environment where the foreach.bed was called, giving access (through <<-) to the local variables.

Below is a simple example is to calculate the size of each window:

sizes.bed <- function(bed) {
  N = dim(bed)[1]
  sizes = vector(mode="integer", length=N)

foreach.bed(bed, function(i, chrom, start, end, strand) {
    sizes[i] <<- end - start
  })

  return(sizes)
}
sizes.bed(bed)

## [1] 3 2 3 3

Everything is wrapped into a function sizes.bed. In the example, N returns the number of rows in the BED data.frame, sizes creates a vector of length Nof zeros. This is used in the foreach.bed function as the return. When foreach.bed is called, the BED data.frame and the function are passed as inputs. func iterates through all i’s calculating the window size, end - start, and setting the corresponding place in the vector, sizes[i], equal to the window size. sizes.bed then returns the vector sizes. The result is a vector of length N of window sizes.

Obviously,sizes = bed[,3] - bed[,2] is much faster, but the function can be designed as complicated as necessary.

bed.region

bed.region.bpQuery.bigWig(bw, bed, strand = NA,
                          op = "sum", abs.value = FALSE,
                          gap.value = 0, bwMap = NULL)
bed.region.probeQuery.bigWig(bw, bed, op = "wavg",
                          abs.value = FALSE, gap.value = NA)

• arguments
  • bw is the pointer of the underlying C object created in load.bigWig
  • bed is a dataframe structured like a bed file with columns for chrom, start and end
  • strand + or - character indicating the strand of the supplied coordinates (bpQuery only)
  • op is a string representing the operation to perform on the interval.
    • sum adds all the counts
    • avg averages the counts
    • min finds the minimum value
    • max finds the maximum value
    • wavg weighted average of the values–only pertains to probeQuery
  • abs.value is a logical argument which determines if the absolute value of the input is performed before the op.
  • gap.value is an integer value that replaces areas that have no overlaps
  • bwMap is only for bp functions and is an R object of class bwMap containing the information required to call steps or regions as unmappable.

This function is an extension of region.bpQuery.bigWig and operates identically, except the the region chrom, start, and end is defined by the corresponding column of a BED data.frame. The output is a vector that is equal in length to the number of rows of the BED data.frame. The vectors’s index corresponds to the BED data.frame row index.

bed.region.bpQuery.bigWig(bw.bp, bed6)

## [1]  3  3  8 12  6

bed.region.probeQuery.bigWig(bw.bp, bed6)

## [1] 1 1 4 4 2

bed6.region

bed6.region.bpQuery.bigWig(bw.plus, bw.minus, bed6,
                 op = "sum", abs.value = FALSE, gap.value = 0, bwMap = NULL)

bed6.region.probeQuery.bigWig(bw.plus, bw.minus, bed6,
                 op = "wavg", abs.value = FALSE, gap.value = NA)

• arguments
  • bw.plus is the R pointer created in load.bigWig and refers to the plus strand signal
  • bw.minus is the R pointer created in load.bigWig and refers to the minus strand signal
  • bed6 is a BED6 style data.frame that specifies a strand value in column 6
  • op is a string representing the operation to perform on the interval.
    • sum adds all the counts
    • avg averages the counts
    • min finds the minimum value
    • max finds the maximum value
    • wavg weighted average of the values–only pertains to probeQuery
  • abs.value is a logical argument which determines if the absolute value of the input is performed before the op.
  • gap.value is an integer value that replaces areas that have no overlaps
  • bwMap is an R object of class bwMap containing the information required to call steps or regions as unmappable.

Using the bw.plus and bw.minus strands, we can evaluate a bed6.region function. First, refer to the query.bigWig for each stranded file as a reference.

bwPlus=load.bigWig('../inst/extdata/bpPlus.bigWig')
bwMinus=load.bigWig('../inst/extdata/bpMinus.bigWig')

query.bigWig(bwMinus, chrom='chr1', start=100, end=115)

##       start end value
##  [1,]   100 101     1
##  [2,]   101 102     1
##  [3,]   102 103     1
##  [4,]   103 104     1
##  [5,]   104 105     1
##  [6,]   110 111     2
##  [7,]   111 112     2
##  [8,]   112 113     2
##  [9,]   113 114     2
## [10,]   114 115     2

query.bigWig(bwPlus, chrom='chr1', start=100, end=115)

##      start end value
## [1,]   105 106     4
## [2,]   106 107     4
## [3,]   107 108     4
## [4,]   108 109     4
## [5,]   109 110     4

Note that each interval is specified by the strand in column 6 of the BED data.frame, and each respective stranded bigWig file is queried for that interval.

#bed6 structure
bed6

##     V1  V2  V3 V4 V5 V6
## 1 chr1 101 104 na  1  +
## 2 chr1 101 104 na  1  -
## 3 chr1 105 107 na  1  +
## 4 chr1 107 110 na  1  +
## 5 chr1 112 115 na  1  -

bed6.region.bpQuery.bigWig(bwPlus, bwMinus, bed6)

## [1]  0  3  8 12  6

bed6.region.probeQuery.bigWig(bwPlus, bwMinus, bed6)

## [1] NA  1  4  4  2

bed.step

The bed.step function operates like the bed.region function, but the intervals are specified in the BED data.frame.

bed.step.bpQuery.bigWig(bw, bed, step,
                strand = NA, op = "sum", abs.value = FALSE, gap.value = 0,
                bwMap = NULL, with.attributes = TRUE, as.matrix = FALSE)

bed.step.probeQuery.bigWig(bw, bed, step,
                op = "wavg", abs.value = FALSE, gap.value = NA,
                bwMap = NULL, with.attributes = TRUE, as.matrix = FALSE)            

• arguments
  • bw is the pointer of the underlying C object created in load.bigWig
  • bed is a dataframe structured like a bed file with columns for chrom, start and end
  • op is a string representing the operation to perform on the interval.
    • sum adds all the counts
    • avg averages the counts
    • min finds the minimum value
    • max finds the maximum value
    • wavg weighted average of the values–only pertains to probeQuery
  • abs.value is a logical argument which determines if the absolute value of the input is performed before the op.
  • gap.value is an integer value that replaces areas that have no overlaps
  • with.attributes is a logical argument that determines if the results are returned annotated with their source components and/or step size.
  • bwMap is an R object of class bwMap containing the information required to call steps or regions as unmappable.
  • as.matrix is a logical argument that will return the results in a matrix format. The BED data.frame has to be a fixed window for all entries.

These functions are an extension of step.bpQuery.bigWig and
step.probeQuery.bigWig and they operate identically, except the the region chrom, start, and end is defined by the corresponding column of a BED data.frame. The regions in the BED data.frame need to be exact multiples of the step. The regions defined within the BED file do not need to be the same size, unless as.matrix = TRUE. The output is a list of vectors, one per query BED interval. An additional argument option is to return the data as a matrix. The matrix argument will only return a matrix if the interval widths in the BED data.frame are identical for each row. This is often useful when querying a fixed window around the center of a genomic feature, such as a transcription factor binding site or transcription start site.

#generate a bed
bed.step=data.frame('chr1',100,106)
bed.step[2,] = c('chr1',109,115)
bed.step[,2] = as.numeric(bed.step[,2])
bed.step[,3] = as.numeric(bed.step[,3])
colnames(bed.step)=c('chrom', 'start', 'end')

bed.step.bpQuery.bigWig(bw.splitprobes, bed.step, step = 2,
                       with.attributes = FALSE)

## [[1]]
## [1] 2 2 5
## 
## [[2]]
## [1] 6 4 4

bed.step.bpQuery.bigWig(bw.splitprobes, bed.step, step = 2,
                       as.matrix=TRUE)

##      [,1] [,2] [,3]
## [1,]    2    2    5
## [2,]    6    4    4
## attr(,"step")
## [1] 2

bed6.step

The bed6.step function operates like the bed6.region function, but the intervals and strand information are specified in the BED data.frame. The follow.strand argument is introduced for the bed6 step function. follow.strand reverses the direction of - strand output. This is commonly set to TRUE when the specific genomic feature in the BED6 file has inherent strandedness. It is useful to know how the counts relate to the orientation of the feature, such as a sequence motif or transcription start site.

bed6.step.bpQuery.bigWig(bw.plus, bw.minus, bed6, step,
                       op = "sum", abs.value = FALSE, gap.value = 0,
                       bwMap = NULL, with.attributes = TRUE, 
                       as.matrix = FALSE, follow.strand = FALSE)

bed6.step.probeQuery.bigWig(bw.plus, bw.minus, bed6, step,
                       op = "wavg", abs.value = FALSE, gap.value = NA,
                       with.attributes = TRUE, as.matrix = FALSE,
                       follow.strand = FALSE)      

• arguments
  • bw.plus is the R pointer created in load.bigWig and refers to the plus strand
  • bw.minus is the R pointer created in load.bigWig and refers to the minus strand
  • chrom is a string referring to what chromosome is referenced
  • start is an integer value designation the starting position
  • end is an integer value designation the ending position
  • op is a string representing the operation to perform on the step.
    • sum adds all the counts
    • avg averages the counts
    • min finds the smallest count
    • max finds the largest count
  • abs.value is a logical argument which determines if the absolute value of the input is performed before the op.
  • gap.value is an integer value that replaces areas that have no overlaps
  • bwMap is an R object of class bwMap containing the information required to call steps or regions as unmappable.
  • with.attributes is a logical argument that determines if the results are returned annotated with their source components and/or step size.
  • as.matrix is a logical argument that will return the results in a matrix format. The BED data.frame has to be a fixed window for all entries.
  • follow.strand is a logical value; if TRUE, in ‘BED’ type queries, the result is a matrix, otherwise it’s a list of vectors, one per query ‘BED’ entry.

To highlight the newly introduced as.matrix and follow.strand functionality, we sum the signal of PRO-seq bigWig files around transcription start sites. The as.matrix option returns a matrix and each row corresponds to the respective row of the BED6 transcription start site (TSS) data.frame and each column is a window equivalent to the step size (50) that spans the 500 base pair interval.

plusPRO=load.bigWig('../inst/extdata/plusPRO.bigWig')
minusPRO=load.bigWig('../inst/extdata/minusPRO.bigWig')
bedTSS=read.table(gzfile("../inst/extdata/TSS.bed.gz"))

bedTSSwindow=fiveprime.bed(bedTSS, upstreamWindow = 249,
                  downstreamWindow = 250)

tss.matrix = bed6.step.bpQuery.bigWig(plusPRO, minusPRO, bedTSSwindow,
                  step = 50, as.matrix=TRUE, follow.strand=TRUE)
colnames(tss.matrix) = seq(-225, 274, by = 50)
rownames(tss.matrix) = bedTSSwindow[,5]

#filter out genes that have no signal
tss.matrix=tss.matrix[rowSums(tss.matrix) != 0,]
head(tss.matrix)

##            -225 -175 -125 -75 -25 25 75 125 175 225
## WASH7P        2    0    0   0   0  0  0   0   0   0
## FO538757.1    0    0    0   0   1  1  1   0   0   0
## MIR6859-2     0    0    0   0   0  0  1   0   0   0
## MTND1P23      0    0    0   0   0  0  0  20   6   1
## MTND2P28      0   16    0   0   4  0  0   0   0   1
## MTCO1P12      0    0    0   0   0  0 24  10   0   0

The follow.strand argument orients the BED6 intervals so that downstream and upstream are relative to the strand in column 6. Although the 50 bp window is relatively wide, we can still observe RNA Polymerase peak just downstream of the TSS using PRO-seq data below with follow.strand=TRUE.

plot(seq(-225, 274, by = 50),colMeans(tss.matrix),
                  xlab = 'stranded TSS position',
                  ylab = 'bigWig intensity', type = "l")

Feature-stranded composite profile. The RNA Polymerase peak results from the accumulation of reads that originate from the coding strand orientation. Therefore, it makes sense to orient the reads relative to the TSS position.

However, if follow.strand is not invoked, the - and + strand signals accumulate upstream and downstream, respectively, of the central position.

tss.matrix.notStrand = bed6.step.bpQuery.bigWig(plusPRO, minusPRO,
                  bedTSSwindow, step = 50, as.matrix=TRUE,
                  follow.strand=FALSE)
tss.matrix.notStrand = tss.matrix.notStrand[rowSums(tss.matrix.notStrand) != 0,]

plot(seq(-225, 274, by = 50),colMeans(tss.matrix.notStrand),
                  xlab = 'absolute TSS position',
                  ylab = 'bigWig intensity', type = "l")

Absolute-stranded composite profile. Since a comparable number of genes are oriented in the plus and minus orientatino relative to teh reference genome, the peak downstream of the TSS is split if absolute genomic coordinate is followed.

quantiles

quantiles.metaprofile invokes R’s quantile function for each column of the matrix; passing tss.matrix, to quantiles.metaprofile we observe that the the majority of windows are zero, because the input bigWig was randomly sampled, so it is a low coverage bigWig. The result of quantiles.metaprofile is a list of quantile values for each column.

quantiles.metaprofile(tss.matrix, quantiles = c(0.95, 0.5, 0.05))

## $top
##  [1]  2  2  2  3  4 18 42 30 17 11
## 
## $middle
##  [1] 0 0 0 0 0 0 1 1 1 0
## 
## $bottom
##  [1] 0 0 0 0 0 0 0 0 0 0

subsampling

The subsampled.quantiles.metaprofile function randomly subsamples the matrix, as determined by fraction and n.samples. A matrix is generated that contains n.samples rows and each row is a column mean from a subsample of the data, then quantiles.metaprofile is applied to the matrix.

subsampled.quantiles.metaprofile(tss.matrix,
            quantiles = c(0.875, 0.5, 0.125), fraction = 0.20,
            n.samples = 1000)

## $top
##  [1] 0.5810508 0.6181525 0.6621675 0.7906571 1.0426653 4.0122136 8.4467672
##  [8] 6.6424151 3.6773769 2.4909468
## 
## $middle
##  [1] 0.5071109 0.5596524 0.5903345 0.6912036 0.9702397 3.7993152 8.1021859
##  [8] 6.2533579 3.4990782 2.3680537
## 
## $bottom
##  [1] 0.4569068 0.5043126 0.5322294 0.6244403 0.9043324 3.5868120 7.7596787
##  [8] 5.9490387 3.3370095 2.2543785

confidence intervals

The confinterval.metaprofile function calculates confidence intervals of the matrix columns.

confinterval.metaprofile(tss.matrix, alpha = 0.05)

## $top
##      -225      -175      -125       -75       -25        25        75       125 
## 0.5413434 0.5804930 0.6209649 0.7338425 0.9973144 3.8684260 8.2255507 6.4406883 
##       175       225 
## 3.5702514 2.4136359 
## 
## $middle
##      -225      -175      -125       -75       -25        25        75       125 
## 0.5175646 0.5602254 0.5974088 0.7033760 0.9705588 3.7923843 8.1009638 6.3097382 
##       175       225 
## 3.5110865 2.3730447 
## 
## $bottom
##      -225      -175      -125       -75       -25        25        75       125 
## 0.4937859 0.5399578 0.5738526 0.6729095 0.9438032 3.7163426 7.9763769 6.1787881 
##       175       225 
## 3.4519216 2.3324535

The result is a list of confidence interval values for each step for the given alpha value. There are 3 different levels of confidence intervals: top, middle, and bottom. These are based on 2 values: the population mean, which is the mean of each column in mat and \(\delta\), defined by:

\(\delta=P(1-\alpha/2) * SE\)

whereby SE is the Standard Error of the column.

Using this \(\delta\) and the following means:

\(top = mean + \delta\)
\(middle = mean\)
\(bottom = mean - \delta\)

bootstraped confidence intervals

bootstrapped.confinterval.metaprofile The bootstrap method produces a confidence interval like confinterval.metaprofiles except that it uses multiple samples to form a distribution. From these, the Central Limit Theorem is applied to determine the confidence interval of each column. This tends to be a more robust calculation of the confidence interval and more n.samples provides a better estimation.

bootstrapped.confinterval.metaprofile(tss.matrix, alpha = 0.05,
                                         n.samples = 300)

## $top
##  [1] 0.5703226 0.6036644 0.6522127 0.7754913 1.0331777 3.9773974 8.4002920
##  [8] 6.6352498 3.6534700 2.4674178
## 
## $middle
##  [1] 0.5175646 0.5602254 0.5974088 0.7033760 0.9705588 3.7923843 8.1009638
##  [8] 6.3097382 3.5110865 2.3730447
## 
## $bottom
##  [1] 0.4544926 0.5112279 0.5441887 0.6262264 0.9083384 3.6061727 7.7900403
##  [8] 6.0034166 3.3705247 2.2798732

metaprofile

metaprofile.bigWig creates a class object of the data, which can be passed to plot.profile.bigWig. Set quantiles.metaprofile using the profile.op = quantiles.metaprofile. We will describe matrix.op = NULL in a later section. tss.matrix = bed6.step.bpQuery.bigWig(plusPRO, minusPRO, bedTSSwindow, step = 10, as.matrix=TRUE, follow.strand=TRUE)

metaprofile.bigWig(bedTSSwindow, plusPRO, bw.minus = minusPRO,
                 step = 50, name = "Signal", matrix.op = NULL,
                 profile.op = quantiles.metaprofile)

## $name
## [1] "Signal"
## 
## $X0
## [1] 0
## 
## $step
## [1] 50
## 
## $step
## [1] 50
## 
## $top
##  [1] 0 0 0 0 0 1 2 2 1 1
## 
## $middle
##  [1] 0 0 0 0 0 0 0 0 0 0
## 
## $bottom
##  [1] 0 0 0 0 0 0 0 0 0 0
## 
## attr(,"class")
## [1] "metaprofile"

This function automatically creates the mat variable and will use the default values for the other the inputs. Below we change the profile.op to bootstrapped.confinterval.metaprofile, and the default alpha=0.05 and n.samples=300 were changed to alpha=0.05, and n.samples=1000.

metaprofile.bigWig(bedTSSwindow, plusPRO, bw.minus = minusPRO,
                 step = 50, name = "Signal", matrix.op = NULL,
                 profile.op = bootstrapped.confinterval.metaprofile,
                 alpha=0.05, n.samples=1000)

## $name
## [1] "Signal"
## 
## $X0
## [1] 0
## 
## $step
## [1] 50
## 
## $step
## [1] 50
## 
## $top
##  [1] 0.1939180 0.2058994 0.2205101 0.2616281 0.3490576 1.3469104 2.8375728
##  [8] 2.2347541 1.2332172 0.8339082
## 
## $middle
##  [1] 0.1747550 0.1891594 0.2017143 0.2374940 0.3277079 1.2804937 2.7352800
##  [8] 2.1304750 1.1855138 0.8012555
## 
## $bottom
##  [1] 0.1550780 0.1718419 0.1824529 0.2123020 0.3054410 1.2133448 2.6347477
##  [8] 2.0258953 1.1385519 0.7681351
## 
## attr(,"class")
## [1] "metaprofile"

RPKM

RPKM [Reads Per Kilobase of transcript per Million mapped reads]. This function will scale everything by a factor of

The libSize can be calculated from attributes in the loaded bigWig object. The product bw$basesCovered * bw$mean is the number of raw reads in the bigWig if it is a counts-based bigWig file. Here, the matrix was generated from two stranded bigWig files, so the library size can be calculated from the bigWig objects:
plusPRO$basesCovered * plusPRO$mean + minusPRO$basesCovered * minusPRO$mean

sizeLibrary= (plusPRO$basesCovered * plusPRO$mean) +
              (minusPRO$basesCovered * minusPRO$mean)
rpkm.matrix = rpkm.scale(tss.matrix, step=10, libSize=sizeLibrary)

head(rpkm.matrix)

##                -225     -175 -125 -75        -25       25       75     125
## WASH7P     0.027756 0.000000    0   0 0.00000000 0.000000 0.000000 0.00000
## FO538757.1 0.000000 0.000000    0   0 0.01387800 0.013878 0.013878 0.00000
## MIR6859-2  0.000000 0.000000    0   0 0.00000000 0.000000 0.013878 0.00000
## MTND1P23   0.000000 0.000000    0   0 0.00000000 0.000000 0.000000 0.27756
## MTND2P28   0.000000 0.222048    0   0 0.05551199 0.000000 0.000000 0.00000
## MTCO1P12   0.000000 0.000000    0   0 0.00000000 0.000000 0.333072 0.13878
##                   175      225
## WASH7P     0.00000000 0.000000
## FO538757.1 0.00000000 0.000000
## MIR6859-2  0.00000000 0.000000
## MTND1P23   0.08326799 0.013878
## MTND2P28   0.00000000 0.013878
## MTCO1P12   0.00000000 0.000000

Density to One

densityToOne is a scaling factor that takes each cell in a row of the matrix and divides it by the sum of each row. The total signal in any row equals 1, or NA if the sum of a row is 0. The na.on.zero = TRUE input is used if you want NAs to populate the matrix row when the sum(row)=0. This would happen because dividing by 0 will result in NA. Otherwise if you want 0 to replace NA then na.on.zero=FALSE should be used.

mat1=matrix(c(0,2,23,6,0,9,2,3,0,5,17,68,0,4,24,9), nrow=4, ncol=4)
mat1

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    2    9    5    4
## [3,]   23    2   17   24
## [4,]    6    3   68    9

densityToOne.scale(mat1, na.on.zero = FALSE)

##            [,1]       [,2]      [,3]      [,4]
## [1,] 0.00000000 0.00000000 0.0000000 0.0000000
## [2,] 0.10000000 0.45000000 0.2500000 0.2000000
## [3,] 0.34848485 0.03030303 0.2575758 0.3636364
## [4,] 0.06976744 0.03488372 0.7906977 0.1046512

Max to one

The maxToOne.scale will take the maximum value for each row and set it equal to 1. Every other cell in the row will be divided by the maximum value. Note that if the max=0 then the whole row is set to 0, this avoids NA instances. If the max is equal to the min, the entire row is set to 1.

maxToOne.scale(mat1)

##            [,1]       [,2]      [,3]      [,4]
## [1,] 0.00000000 0.00000000 0.0000000 0.0000000
## [2,] 0.22222222 1.00000000 0.5555556 0.4444444
## [3,] 0.95833333 0.08333333 0.7083333 1.0000000
## [4,] 0.08823529 0.04411765 1.0000000 0.1323529

mat2=matrix(c(0,2,23,6,0,2,2,3,0,2,17,68,0,2,24,9), nrow=4, ncol=4)

maxToOne.scale(mat2)

##            [,1]       [,2]      [,3]      [,4]
## [1,] 0.00000000 0.00000000 0.0000000 0.0000000
## [2,] 1.00000000 1.00000000 1.0000000 1.0000000
## [3,] 0.95833333 0.08333333 0.7083333 1.0000000
## [4,] 0.08823529 0.04411765 1.0000000 0.1323529

Zero to one

zeroToOne.scale sets the minimum value and maximum value of each row to 0 and 1. For every cell the the function applies the formula: \(Xn = \dfrac{Xo − min}{max − min}\)
Like maxToOne.scale, if the max=0 then the whole row is set to 0 and if the max is equal to the min, the entire row is set to 1.

zeroToOne.scale(mat2)

##            [,1] [,2]      [,3]       [,4]
## [1,] 0.00000000    0 0.0000000 0.00000000
## [2,] 1.00000000    1 1.0000000 1.00000000
## [3,] 0.95454545    0 0.6818182 1.00000000
## [4,] 0.04615385    0 1.0000000 0.09230769

metaprofile with a matrix.op

Any of these scaling factors can be passed to the metaprofile.bigWig function using the matrix.op argument:

metaprofile.bigWig(bed6, bw.plus, bw.minus = bw.minus, step = 1,
                  name = "Signal", matrix.op = zeroToOne.scale)

The matrix.op and profile.op options are processed by metaprofile.bigWig and the resultant metaprofile instance can be passed to plot.metaprofile. We chose X0 = 255 to center the plot because the bedTSSwindow was 500 and the step is specified as 10–a quick X0 calculation for a composite profile centered on a feature is the sum of the window and the step divided by 2.

composite.plus = metaprofile.bigWig(bedTSSwindow, plusPRO,
                 bw.minus = minusPRO, step = 10, name = "Signal",
                 matrix.op = zeroToOne.scale,
                 profile.op = bootstrapped.confinterval.metaprofile,
                 alpha = 0.05, n.samples = 200)

composite.min = metaprofile.bigWig(bedTSSwindow, minusPRO,
                 bw.minus = plusPRO, step = 10, name = "Signal",
                 matrix.op = zeroToOne.scale,
                 profile.op = bootstrapped.confinterval.metaprofile,
                 alpha = 0.05, n.samples = 200)

plot.metaprofile(composite.plus, minus.profile = composite.min, X0 = 255,
      draw.error = TRUE, col = c("red", "blue", "lightgrey", "lightgrey"),
      xlab = "Distance (bp)", ylab = "PRO-seq signal")

Composite PRO-seq profiles.