R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Estimate the Dispersion factor in a beta-binomial model.

estimateDispersion

R Documentation

Estimate the Dispersion factor in a beta-binomial model.

Description

This function estimates the dispersion factor in a beta-binomial model of the nucleotide counts. This model assumes that the count for nucleotide j at position i is distributed after a beta-binomial X_ib ~ BB(n_i; alpha, beta_ij), where n_i is the coverage. The base and nucleotide specific parameter beta_ij is estimated from the local mean by the method-of-moments estimate, alpha is a shared overdispersion parameter. It is estimated via a numerical optimization of the likelihood under the null-hypothesis.

Usage

estimateDispersion(test, control, ...)

## S4 method for signature 'deepSNV,missing'
estimateDispersion(test, control, alternative = NULL, interval = c(0,1000))

## S4 method for signature 'matrix,matrix'
estimateDispersion(test, control, alternative = NULL, interval = c(0,1000))

Arguments

`test`	Either a deepSNV object, or a matrix with the test counts.
`control`	Missing if test is a deepSNV object, otherwise missing.
`alternative`	The alternative to be tested. One of "greater", "less", "two-sided" (default). If test is a deepSNV object, automatically taken from the corresponding slot if unspecified.
`interval`	The interval to be screened for the overdispersion factor. Default (0,1000).
`...`	Additional param passed to specific methods

Value

A deepSNV-class object if the input was a deepSNV object. Otherwise the loglikelihood and the estimated parameter.

Author(s)

Moritz Gerstung

Examples

data("RCC", package="deepSNV")
plot(RCC)
summary(RCC)[,1:6]
RCC.bb = estimateDispersion(RCC, alternative = "two.sided")
summary(RCC.bb)

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(deepSNV)
Loading required package: parallel
Loading required package: Rhtslib
Rhtslib htslib version 1.1
Loading required package: IRanges
Loading required package: BiocGenerics

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
Loading required package: stats4

Attaching package: 'S4Vectors'

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

    colMeans, colSums, expand.grid, rowMeans, rowSums

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: Biostrings
Loading required package: XVector
Loading required package: VGAM
Loading required package: splines
Loading required package: VariantAnnotation
Loading required package: Rsamtools

Attaching package: 'VariantAnnotation'

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

    tabulate


Attaching package: 'deepSNV'

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

    dbetabinom, pbetabinom

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

    normalize

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/deepSNV/estimateDispersion-methods.Rd_%03d_medium.png", width=480, height=480)
> ### Name: estimateDispersion
> ### Title: Estimate the Dispersion factor in a beta-binomial model.
> ### Aliases: estimateDispersion estimateDispersion,deepSNV,missing-method
> ###   estimateDispersion,matrix,matrix-method
> 
> ### ** Examples
> 
> data("RCC", package="deepSNV")
> plot(RCC)
> summary(RCC)[,1:6]
     chr      pos ref var         p.val      freq.var
6  chr10 89710231   C   T  1.812318e-06 -0.0522069105
1  chr17  7512879   G   A  4.923709e-03  0.0095680379
10 chr17  7513782   A   G  6.137099e-57  0.0523098637
15  chr3 10158255   A   -  7.414419e-27  0.0374605077
4   chr3 10158274   C   T 6.596999e-152 -0.1420210172
2   chr3 10158337   G   A  8.455105e-09 -0.1458150689
7   chr3 10163012   -   C 2.512916e-211  0.1242571351
16  chr3 10163206   T   -  0.000000e+00  0.2440387988
17  chr3 10163208   G   -  1.908223e-02  0.0007724394
11  chr3 10163428   T   G  1.026389e-79 -0.1111055718
8   chr3 10166219   G   C  6.538408e-74  0.1591560973
3   chr3 10166943   G   A 1.406981e-158 -0.1326056510
12  chr3 10167220   C   G  2.960126e-36  0.0046531880
13  chr3 10167672   A   G  0.000000e+00  0.1399833662
5   chr3 10167709   C   T 6.203881e-304 -0.1440968910
9   chr3 10167762   T   C  0.000000e+00 -0.1323353964
14  chr3 10168683   T   G 2.134656e-304 -0.1324914089
> RCC.bb = estimateDispersion(RCC, alternative = "two.sided")
Note: The initial object used a binomial model. Will be changed to beta-binomial.
Estimated dispersion factor 136.926849574802 
> summary(RCC.bb)
   chr      pos ref var         p.val      freq.var sigma2.freq.var n.tst.fw
2 chr3 10158274   C   T  4.185487e-03 -0.1420210172    1.553566e-05     3458
7 chr3 10163206   T   - 6.777322e-304  0.2440387988    4.701210e-06     4884
8 chr3 10163208   G   -  2.693887e-02  0.0007724394    1.704750e-08       14
1 chr3 10166943   G   A  9.500067e-03 -0.1326056510    2.037310e-05    13302
5 chr3 10167220   C   G  5.461071e-30  0.0046531880    6.646192e-08      108
3 chr3 10167709   C   T  3.466478e-02 -0.1440968910    1.961392e-05     5149
4 chr3 10167762   T   C  8.715720e-03 -0.1323353964    1.549514e-05    10478
6 chr3 10168683   T   G  4.237396e-02 -0.1324914089    1.880836e-05     8747
  cov.tst.fw n.tst.bw cov.tst.bw n.ctrl.fw cov.ctrl.fw n.ctrl.bw cov.ctrl.bw
2      12927    20064      54095      3595        9298     20002       38568
7      22404     8001      30140        15       13853        26       20768
8      19748       21      25563         0       11943         0       18260
1      38634     5043      15262     12414       26061      3500        7585
5      19249      224      51735         0        9809         1       31993
3      28719    14640      30063      4401       15142     12242       19477
4      28367    12397      42504      8947       17750      9984       23847
6      26472    11320      32886      7981       16598      8905       19287
      raw.p.val
2  7.063875e-08
7 1.143813e-308
8  4.546491e-07
1  1.603333e-07
5  9.216686e-35
3  5.850398e-07
4  1.470958e-07
6  7.151483e-07
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>