Compute the log odd ratio of RIP over background.

Description

The RIPScore is computed as the log odd ratio of the posterior for the RIP state (z_i = 2) over the posterior for the background state (z_i = 1)

Usage

computeLogOdd(nbhGR)

Arguments

Details

To assess the statistical significance of the RIP predictions, we assign each bin a RIPScore defined as the log odd ratio of the posterior for the RIP state (z_i = 2) over the posterior for the background state (z_i = 1). When control is available, the RIPScore is updated as the difference between the RIPScores evaluated separately for RIP and control libraries. The scoring system captures the model confidence for the RIP state of each bin in the RIP library penalized by the false confidence for the RIP state of the same bin in the control library. In addition, RIPScore obviates scaling of read counts. Since sequencing depth usually differs between RIP and control libraries, scaling is necessary if the statistical score were derived from the read count differences. On the other hand, simplistic linear scaling may distort the data.

Value

A vector of log odd scores for each bin in nbhGR.

Author(s)

Yue Li

See Also

seekRIP, scoreMergedBins, logScoreWithoutControl, logScoreWithControl

Examples

# Retrieve system files
extdata.dir <- system.file("extdata", package="RIPSeeker") 

bamFiles <- list.files(extdata.dir, ".bam$", recursive=TRUE, full.names=TRUE)

bamFiles <- grep("PRC2", bamFiles, value=TRUE)

# Retrieve system files
extdata.dir <- system.file("extdata", package="RIPSeeker") 

bamFiles <- list.files(extdata.dir, ".bam$", recursive=TRUE, full.names=TRUE)

bamFiles <- grep("PRC2", bamFiles, value=TRUE)

alignGal <- getAlignGal(bamFiles[1], reverseComplement=TRUE, genomeBuild="mm9")

alignGR <- as(alignGal, "GRanges")

alignGRList <- GRangesList(as.list(split(alignGR, seqnames(alignGR))))

################ run main function for HMM inference on a single chromosome ################
nbhGR <- mainSeekSingleChrom(alignGR=alignGRList$chrX, K = 2, binSize=1e5)

ripscore <- computeLogOdd(nbhGR)

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(RIPSeeker)
Loading required package: S4Vectors
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


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: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: Rsamtools
Loading required package: Biostrings
Loading required package: XVector
Loading required package: GenomicAlignments
Loading required package: rtracklayer
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/RIPSeeker/computeLogOdd.Rd_%03d_medium.png", width=480, height=480)
> ### Name: computeLogOdd
> ### Title: Compute the log odd ratio of RIP over background.
> ### Aliases: computeLogOdd
> 
> ### ** Examples
> 
> # Retrieve system files
> extdata.dir <- system.file("extdata", package="RIPSeeker") 
> 
> bamFiles <- list.files(extdata.dir, ".bam$", recursive=TRUE, full.names=TRUE)
> 
> bamFiles <- grep("PRC2", bamFiles, value=TRUE)
> 
> # Retrieve system files
> extdata.dir <- system.file("extdata", package="RIPSeeker") 
> 
> bamFiles <- list.files(extdata.dir, ".bam$", recursive=TRUE, full.names=TRUE)
> 
> bamFiles <- grep("PRC2", bamFiles, value=TRUE)
> 
> alignGal <- getAlignGal(bamFiles[1], reverseComplement=TRUE, genomeBuild="mm9")
Processing /home/ddbj/local/lib64/R/library/RIPSeeker/extdata/PRC2/SRR039210_processed_tophat/accepted_hits_noDup_sel_chrX.bam ... All hits are returned with flags.
> 
> alignGR <- as(alignGal, "GRanges")
> 
> alignGRList <- GRangesList(as.list(split(alignGR, seqnames(alignGR))))
> 
> ################ run main function for HMM inference on a single chromosome ################
> nbhGR <- mainSeekSingleChrom(alignGR=alignGRList$chrX, K = 2, binSize=1e5)

**********
chrX:
**********

***0. Using predefined binSize of 100000 bp.
***1. Traning NB HMM to derive posterior (and Viterbi state sequence:)

***1. Initializing negative binomial HMM (nbh) with 2 states:


Starting NB mixture model (nbm_em) for K=2 clusters:
Iteration 0:	-9843.544
Iteration 1:	-7325.446
Iteration 2:	-6986.980
Iteration 3:	-6747.130
Iteration 4:	-6590.919
Iteration 5:	-6509.158
Iteration 6:	-6455.276

***2. Traininig nbh with forward-backward algorithm:

Iteration 0:	-6410.144
Iteration 1:	-6324.948
Iteration 2:	-6292.593
Iteration 3:	-6265.340
Iteration 4:	-6241.037
Iteration 5:	-6219.898
Iteration 6:	-6204.007
Iteration 7:	-6192.605
Iteration 8:	-6184.278
Iteration 9:	-6178.123
> 
> ripscore <- computeLogOdd(nbhGR)
> 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>