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

Source: R/calculateStats.R

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.

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

Leonardo Collado-Torres

Examples

## Collapse the coverage information
collapsedFull <- collapseFullCoverage(list(genomeData$coverage),
    verbose = TRUE
)
#> 2023-05-07 06:01:09.436288 collapseFullCoverage: Sorting fullCov
#> 2023-05-07 06:01:09.438749 collapseFullCoverage: Collapsing chromosomes information by sample

## Calculate library size adjustments
sampleDepths <- sampleDepth(collapsedFull, probs = c(0.5), verbose = TRUE)
#> 2023-05-07 06:01:09.440814 sampleDepth: Calculating sample quantiles
#> 2023-05-07 06:01:09.44656 sampleDepth: Calculating sample adjustments

## 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")
#> 2023-05-07 06:01:09.708332 calculateStats: calculating the F-statistics
fstats
#> numeric-Rle of length 1434 with 361 runs
#>   Lengths:           1           5           3 ...           1           1
#>   Values : 5.10948e+00 3.85229e+00 5.12448e-01 ... 0.032237427 3.997655116
if (FALSE) {
## Compare vs pre-packaged F-statistics
library("testthat")
expect_that(fstats, is_equivalent_to(genomeFstats))
}