Nonpareil.curve {Nonpareil}

Generates a Nonpareil curve from an .npo file

Description

Generates a Nonpareil curve from an .npo file

Usage

Nonpareil.curve(
  file,
  plot = TRUE,
  label = NA,
  col = NA,
  enforce.consistency = TRUE,
  star = 95,
  correction.factor = TRUE,
  weights.exp = NA,
  skip.model = FALSE,
  ...
)

Arguments

file

Path to the .npo file, containing the read redundancy.

plot

Determines if the plot should be produced. If FALSE, it computes the coverage and the model wihtout plotting.

label

Name of the dataset. If NA, it is determined by the file name.

col

Color of the curve. If NA, a random color is assigned (even if plot = FALSE).

enforce.consistency

If TRUE, it fails verbosely on insufficient data, otherwise it warns about the inconsistencies and attempts the estimations.

star

Objective coverage in percentage; i.e., coverage value considered near-complete.

correction.factor

Should the overlap-dependent (or kmer-length-dependent) correction factor be applied? If FALSE, redundancy is assumed to equal coverage.

weights.exp

Vector of values to be tested (in order) as exponent of the weights distribution. If the model fails to converge, sometimes manual modifications in this parameter may help. By default (NA), five different values are tested in the following order: For linear sampling, -1.1, -1.2, -0.9, -1.3, -1. For logarithmic sampling (-d option in Nonpareil), 0, 1, -1, 1.3, -1.1, 1.5, -1.5.

skip.model

If set, skips the model estimation altogether.

...

Any additional parameters passed to plot.Nonpareil.Curve.

Value

Returns invisibly a Nonpareil.Curve object

Examples

# Generate a Nonpareil plot
file <- system.file("extdata", "LakeLanier.npo", package = "Nonpareil")
np <- Nonpareil.curve(file)

# Produce the same plot but using powers of 1,000bp as X axis labels
Nonpareil.curve(file, xaxt = "n", xlab = "Sequencing Effort")
axis(
  1L, at = 10L^seq(3L, 12L, by = 3L),
  labels = paste(1L, c("Kbp", "Mbp", "Gbp", "Tbp"))
)

# Show the estimated values
print(np)

# Predict coverage for 20Gbp
predict(np, 20e9)

# Obtain the Nd diversity index
np$diversity