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, ...)



A list with $coverageProcessed, $mclapplyIndex, and $position normally generated using preprocessCoverage.


A list with $mod and $mod0 normally generated using makeModels.


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:


If TRUE basic status updates will be printed along the way.


This argument is passed to fstats.apply and should be the same as the one used in preprocessCoverage. Default: 32.


Has to be either 'Matrix' (default), 'Rle' or 'regular'. See details in fstats.apply.


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.


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


Leonardo Collado-Torres


## 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
#> 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
expect_that(fstats, is_equivalent_to(genomeFstats))