calculateStats {derfinder} | R Documentation

## Calculate F-statistics at base pair resolution from a loaded BAM files

### Description

After defining the models of interest (see makeModels) and pre-processing the data (see preprocessCoverage), use calculateStats to calculate the F-statistics at base-pair resolution.

### Usage

    calculateStats(coveragePrep, models, lowMemDir = NULL, ...)

### Arguments

`coveragePrep` | A list with `$coverageProcessed`, `$mclapplyIndex`, and `$position` normally generated using preprocessCoverage.
`models` | A list with `$mod` and `$mod0` normally generated using makeModels.
`lowMemDir` | The directory where the processed chunks are saved when using preprocessCoverage with a specified `lowMemDir`.
`...` | Arguments passed to other methods and/or advanced arguments. Advanced arguments: verbose If `TRUE` basic status updates will be printed along the way. scalefac This argument is passed to fstats.apply and should be the same as the one used in preprocessCoverage. Default: 32. method Has to be either 'Matrix' (default), 'Rle' or 'regular'. See details in fstats.apply. adjustF A single value to adjust that is added in the denominator of the F-stat calculation. Useful when the Residual Sum of Squares of the alternative model is very small. Default: 0. Passed to define_cluster.

### Value

A numeric Rle with the F-statistics per base pair that passed the cutoff.

### Author(s)

Leonardo Collado-Torres

### See Also

makeModels, preprocessCoverage

### Examples

    ## Collapse the coverage information
    collapsedFull <- collapseFullCoverage(list(genomeData$coverage),
        verbose = TRUE
    )

    ## Calculate library size adjustments
    sampleDepths <- sampleDepth(collapsedFull, probs = c(0.5), verbose = TRUE)

    ## Build the models
    group <- genomeInfo$pop
    adjustvars <- data.frame(genomeInfo$gender)
    models <- makeModels(sampleDepths, testvars = group, adjustvars = adjustvars)

    ## Preprocess the data
    prep <- preprocessCoverage(genomeData,
        cutoff = 0, scalefac = 32,
        chunksize = 1e3, colsubset = NULL
    )

    ## Run the function
    fstats <- calculateStats(prep, models, verbose = TRUE, method = "regular")
    fstats
    ## Not run:
    ## Compare vs pre-packaged F-statistics
    library("testthat")
    expect_that(fstats, is_equivalent_to(genomeFstats))

    ## End(Not run)

[Package derfinder version 1.28.0 Index]