M3D (version 1.6.2)

pvals: Computes p-values

Description

Returns p-values for each region reflecting the probability of observing the mean test-statistic of the between group comparisons among the inter-replicate comparisons.

Usage

pvals(rrbs, CpGs, MMD, group1, group2, smaller = FALSE, comparison = "allReps", method = "empirical", closePara = 0.005)

Arguments

rrbs

An rrbs object containing methylation and coverage data as created using the BiSeq pacakge

CpGs

A GRanges object with each row being a testing region

MMD

A matrix containing the M3D test-statistic, the difference the full and methylation blind metrics, for each region in the CpGs object. Each column is a comparison between two samples, which are described in the column names.

group1

The name of the first group for the comparison. This is stored in colData(rrbs)

group2

The name of the second group for the comparison. This is stored in colData(rrbs)

smaller

Determines whether the p-value is computed whether the test-statistic is greater or lesser than inter-replicate values. For our purposes, it should be set to FALSE.

comparison

Details which groups we are using to define our empirical testing distribution. The default is to use all of them, however, should the user find one group contains unusually high variability, then that group can be selected. Values are either 'allReps', 'Group1' or 'Group2'.

method

Determines which method is used to calculate p-values. 'empirical' uses the empirical distribution directly, without modelling. This is the default. 'model', fits an exponential distribution to the tail of our null distribution.

closePara

Sets a threshold for how close the exponential curve should fit the empirical distribution in the 'model' method. If the method produces errors, consider raising this parameter.

Value

Returns a list P, with 2 entries. 'FDRmean' is the Benjamini-Hochberg adjusted p-values. The unadjusted p-values are stored in 'Pmean'.

Examples

data(rrbsDemo)
data(CpGsDemo)
data(MMDlistDemo)
M3Dstat <- MMDlistDemo$Full-MMDlistDemo$Coverage
group1 <- unique(colData(rrbsDemo)$group)[1]
group2 <-unique(colData(rrbsDemo)$group)[2]
PDemo <- pvals(rrbsDemo, CpGsDemo, M3Dstat,
            group1, group2, smaller=FALSE,comparison='allReps')