Last data update: 2014.03.03

R: arithmetic operation on an interval
cov.intervalR Documentation

arithmetic operation on an interval

Description

This method computes a numeric vector with the result of a certain arithmetic operation with the coverage for two particular BAM or BigWIG files within a specific genomic interval defined by the user

Usage

  ## S4 method for signature 'CoverageBamFile,CoverageBamFile'
cov.interval(tr,ctl,normalization,chr,start,end,bin_width,do)
  ## S4 method for signature 'CoverageBigWigFile,CoverageBigWigFile'
cov.interval(tr,ctl,normalization,chr,start,end,bin_width,do)

Arguments

tr

An instance of a CoverageBamFile or a CoverageBigWigFile object used to compute an arithmetic operation with the coverage values for the genomic interval defined by the chr, start and end parameters

ctl

An instance of a CoverageBamFile or a CoverageBigWigFile object used to compute a certain arithmetic operation with the coverage values for the genomic interval defined by the chr, start and end parameters

normalization

Normalize the coverage values in the resulting numeric vector. Possible normalization types are rpm (Read per million). Default=NULL

chr

chromosome of the genomic interval

start

start coordinate of the genomic interval. If the start and end arguments are not set, then the entire chromosome defined with the chr argument will be used for the calculations

end

end coordinate of the genomic interval

bin_width

calculate the average coverage within bins of size defined by this argument. Default is 1 nucleotide

do

specify the arithmetic operation to perform on the genomic interval. Possible values are 'log2ratio','ratio' or 'substraction'.Default=log2ratio

Details

The cov.interval function receives 2 CoverageBamFiles or CoverageBigWigFiles objects and returns a numeric vector with the result of doing a certain arithmetic operation with the two files using the coverage values for a certain position or genomic bin that is included in the genomic interval defined by the user using the chr, start and end arguments. The size of the bins is controlled by the bin_width argument and for each bin the average coverage value is computed and used for the different arithmetic operations.The arithmetic operation to perform is set with the do argument.

If the normalization parameter is set to 'rpm' then the number of reads aligned into the reference is required to perform the normalization. This number can be provided using the reads_mapped argument of the CoverageBamFile or CoverageBigWigFile objects, or it is automatically calculated if the function is used with a CoverageBamFile object. NOTE: If the reads_mapped argument is set then the automatic calculation will NOT be performed and the provided number will be taken as correct for all downstream calculations.

Author(s)

Ernesto Lowy <ernestolowy@gmail.com>

References

BAM format specification: http://samtools.sourceforge.net/SAMv1.pdf

WIG format specification: https://genome.ucsc.edu/FAQ/FAQformat.html

Examples

##BAM files
  
#get treatment and control test files
treat1file<-system.file("extdata","treat.bam",package="CoverageView")
control1file<-system.file("extdata","ctrl.bam",package="CoverageView")
  
#create two CoverageBamFile objects representing single-end alignments
trm1<-CoverageBamFile(treat1file,run_type="single")
ctl1<-CoverageBamFile(control1file,run_type="single")
  
#calculate the ratio of the coverages for the defined genomic interval using a bin_width equal to 10 nts
cov1=cov.interval(trm1,ctl1,chr="chrI",start=1,end=100,bin_width=10,do="ratio")
  
#create a WIG file with the obtained vector with the ratios
outfolder=system.file("extdata",package="CoverageView")
an_outfile1=paste(outfolder,"out.wig",sep="/")
export.wig(cov1,outfile=an_outfile1)
  
##BigWIG files
  
#get a treatment and control test files
treat2file<-system.file("extdata","treat.bw",package="CoverageView")
control2file<-system.file("extdata","ctrl.bw",package="CoverageView")
  
#create the 'treatment' and 'control' CoverageBigWigFile objects
trm2<-CoverageBigWigFile(path=treat2file,reads_mapped=28564510)
ctl2<-CoverageBigWigFile(path=control2file,reads_mapped=26713667)
  
#calculate the ratio of the coverages for the defined genomic interval
cov2=cov.interval(trm2,ctl2,bin_width=1,chr="chrI",start=1,end=1000,do='ratio')

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(CoverageView)
Loading required package: Rsamtools
Loading required package: GenomeInfoDb
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: S4Vectors

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: IRanges
Loading required package: GenomicRanges
Loading required package: Biostrings
Loading required package: XVector
Loading required package: rtracklayer

Attaching package: 'CoverageView'

The following object is masked from 'package:rtracklayer':

    export.wig

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/CoverageView/cov.interval.Rd_%03d_medium.png", width=480, height=480)
> ### Name: cov.interval
> ### Title: arithmetic operation on an interval
> ### Aliases: cov.interval cov.interval,CoverageBamFile-method
> ###   cov.interval,CoverageBigWigFile-method
> ###   cov.interval,CoverageBamFile,CoverageBamFile-method
> ###   cov.interval,CoverageBigWigFile,CoverageBigWigFile-method
> 
> ### ** Examples
> 
> ##BAM files
>   
> #get treatment and control test files
> treat1file<-system.file("extdata","treat.bam",package="CoverageView")
> control1file<-system.file("extdata","ctrl.bam",package="CoverageView")
>   
> #create two CoverageBamFile objects representing single-end alignments
> trm1<-CoverageBamFile(treat1file,run_type="single")
> ctl1<-CoverageBamFile(control1file,run_type="single")
>   
> #calculate the ratio of the coverages for the defined genomic interval using a bin_width equal to 10 nts
> cov1=cov.interval(trm1,ctl1,chr="chrI",start=1,end=100,bin_width=10,do="ratio")
>   
> #create a WIG file with the obtained vector with the ratios
> outfolder=system.file("extdata",package="CoverageView")
> an_outfile1=paste(outfolder,"out.wig",sep="/")
> export.wig(cov1,outfile=an_outfile1)
>   
> ##BigWIG files
>   
> #get a treatment and control test files
> treat2file<-system.file("extdata","treat.bw",package="CoverageView")
> control2file<-system.file("extdata","ctrl.bw",package="CoverageView")
>   
> #create the 'treatment' and 'control' CoverageBigWigFile objects
> trm2<-CoverageBigWigFile(path=treat2file,reads_mapped=28564510)
> ctl2<-CoverageBigWigFile(path=control2file,reads_mapped=26713667)
>   
> #calculate the ratio of the coverages for the defined genomic interval
> cov2=cov.interval(trm2,ctl2,bin_width=1,chr="chrI",start=1,end=1000,do='ratio')
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>