# INTERNAL
Nonpareil.__init_globals <- function(
### Internal Nonpareil function.
soft=TRUE){
if(soft && exists('Nonpareil.LastFactor')) return();
Nonpareil.LastFactor <<- NULL;
Nonpareil.LastXmax <<- NULL;
Nonpareil.OpenColors <<- c();
Nonpareil.OpenNames <<- c();
}
# MODEL FUNCTIONS
#Nonpareil.f <- function(x,a,b){ return( ( 1 - exp(-x/(2^a)) )^b ) }
#Nonpareil.antif <- function(y,a,b){ return( -(2^a)*log(1 - y^(1/b)) ) }
Nonpareil.f <- function(
### Function of the projected model.
x,
### Values of sequencing effort (typically in bp).
a,
### Parameter alpha of the gamma CDF.
b
### Parameter beta of the Gamma CDF.
){ return(pgamma(log(x+1),a, b))
### Predicted values of abundance-weighted average coverage.
}
Nonpareil.antif <- function(
### Complement function of Nonpareil.f.
y,
### Values of abundance-weighted average coverage
a,
### Parameter alpha of the gamma CDF.
b
### Parameter beta of the gamma CDF.
){ return(exp(qgamma(y,a,b))-1)
### Estimated sequencing effort.
}
#Nonpareil.f <- function(x,a,b){ lx<-log(max(1, x)); return(exp((lx^a)/(b + lx^a))) }
#Nonpareil.antif <- function(y,a,b){ ly<-log(max(1e-10, y)); return(exp((b*ly/(1-ly))^(1/a))) }
# NONPAREIL CURVES
Nonpareil.curve.batch <- function(
### Generates a single plot (and a single table of estimates) with multiple Nonpareil
### curves, from different .npo files.
files,
### Vector of characters with the paths to the .npo files.
overlap,
### Value of the '-L' parameter (in Nonpareil, the default is 50). It can be a number (if all the curves
### were generated with the same value) or a vector (in the same order of 'files'). See the 'overlap' value
### of 'Nonpareil.curve()'. Use only with Nonpareil < v2.0.
r=NA,
### Values of 'r' in 'Nonpareil.curve()' (red value).
g=NA,
### Values of 'g' in 'Nonpareil.curve()' (green value).
b=NA,
### Values of 'b' in 'Nonpareil.curve()' (blue value).
libnames=NA,
### Vector of names of the libraries (corresponding to 'libname' in 'Nonpareil.curve()'). It must be characters,
### not factors.
read.lengths=NA,
### A vector of numbers indicating the length of the reads (corresponding to 'read.length' in 'Nonpareil.curve()').
### Use only with Nonpareil < v2.0.
col=NA,
### Values of 'col' in 'Nonpareil.curve()'.
...
### Any other parameter accepted by 'Nonpareil.curve()' is supported.
){
if(!is.vector(files)) files = as.vector(files);
if(missing(overlap)) overlap=rep(NULL, length(files));
new=TRUE;
for(i in 1:length(files)){
o = Nonpareil.curve(as.character(files[i]), overlap, r=r[i], g=g[i], b=b[i], col=col[i],
libname=ifelse(is.na(libnames[1]), NA, as.character(libnames[i])), new=new,
read.length=ifelse(is.na(read.lengths[1]), NA, read.lengths[i]), ...);
if(new){
out.m = matrix(NA, ncol=length(o), nrow=length(files));
colnames(out.m) <- rownames(as.matrix(o));
if(!is.na(libnames[1])) rownames(out.m) <- libnames;
}
out.m[i, ] <- as.numeric(o);
new = FALSE;
}
return(out.m);
### A dataframe containing the values generated by 'Nonpareil.curve()'.
}
Nonpareil.curve <- function(
### Generates a Nonpareil curve from a .npo file.
file,
### Path to the .npo file, containing the read redundancy.
overlap=NULL,
### Value of the '-L' parameter (in Nonpareil, the default is 50). If not set, it tries to find the value in the
### .npo file (supported in Nonpareil >= 2.0), or fails with an error message. Use only with Nonpareil < v2.0.
factor=1,
### Multiplier of the sequencing effort. This can be used to express the sequencing effort in units other than
### base pairs (bp). For example, to express sequencing effort as Gbp, use factor=1e-9. This can also affect
### the fit of the model, and it's considered EXPERIMENTAL.
plotDispersion=NA,
### Indicates if (and how) dispersion of the replicates should be plotted. It requires modelOnly=FALSE to take
### effect. It can be NA, in which case no dispersion is plotted, or any of the following strings: 'sd' (one
### standard deviation around the mean), 'ci95' (95% confidence interval), 'ci90' (90% confidence interval),
### 'ci50' (50% confidence interval), 'iq' (inter-quartile range).
returnModelValues=FALSE,
### If TRUE, returns the coordinates of the model as model.x and model.y.
returnModelParameters=FALSE,
### If TRUE, returns the model itself as model.
xmax=10e12,
### Maximum sequencing effort to plot.
xmin=1e3,
### Minimum sequencing effort to plot
ymax=1,
### Maximum coverage to plot.
ymin=1e-6,
### Minimum coverage to plot.
xlab=NULL,
### Label of the X-axis. If NULL, it's set to sequencing effort and the units (see factor).
ylab=NULL,
### Label of the Y-axis. If NULL, it's set to Estimated average coverage.
col=NA,
### The color of the curve. If passed, it overrides `r`, `g`, and `b`.
r=NA,
### Red component of the curve's color. If NA, it's randomly set. If <=1, it's assumed to be in the range [0,1]; if
### >1, it's assumed to be in the range [0,256].
g=NA,
### Green component of the curve's color. Same as `r`.
b=NA,
### Blue component of the curve's color. Same as `r`.
new=TRUE,
### If FALSE, it attempts to use a previous (active) canvas to plot the curve.
plot=TRUE,
### Determines if the plot should be produced. If FALSE, it still computes the coverage and the model.
libname=NA,
### Name of the library. If NA, it's determined by the file name. This is useful if you plan to call Nonpareil.legend().
modelOnly=FALSE,
### If TRUE, the rarefied data is not presented, only the fitted model.
plotModel=TRUE,
### If FALSE, the model is not plotted (but it's still computed).
plotDiversity=FALSE,
### If TRUE, the diversity estimate is plotted as a small arrow below the Nonpareil curve.
curve.lwd=2,
### Line width of the Nonpareil curve.
curve.alpha=0.4,
### Alpha value (from 0 to 1) of the Nonpareil curve.
model.lwd=1,
### Line width of the model.
model.alpha=1,
### Alpha value (from 0 to 1) of the model.
log='x',
### Axis to plot in logarithmic scale. It can be 'x' (sequencing effort, default), 'y' (coverage), 'xy' (both logarithmic),
### or '' (both linear).
data.consistency=TRUE,
### If TRUE, it checks the consistency of the data before plotting.
useValue='mean',
### Controls how the replicates are to be summarized at each point of sequencing effort. It can be any of: 'mean' (average of
### the replicates), 'median' (median), 'ub' (upper boundary of the 95% confidence interval), 'lb' (lower boundary of the 95%
### confidence interval), 'q1' (quartile 1), 'q3' (quartile 3). Note that the quartile 2 is also 'median'.
star=95,
### Objective coverage (in percentage). By default: 95, which means the sequencing effort required to reach 95% average
### coverage is to be estimated.
read.length=NA,
### Length of the reads. Use only with Nonpareil < v2.0.
weights.exp=NA,
### 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.
...
### Any other parameters accepted by plot().
){
# Create environment
Nonpareil.__init_globals(!new);
# Examine consistency
if(is.null(file)) stop('The file argument is mandatory');
if(!new && (is.null(Nonpareil.LastFactor) || is.null(Nonpareil.LastXmax)))
stop('No previous plot found, please use new=TRUE if this is the first plot.');
# Read metadata (.npo header)
log_sampling <- 0;
meta_data <- gsub('^# @', '', grep("^# @", readLines(file), value=TRUE));
keys <- gsub(': .*', '', meta_data);
vals <- gsub('.*: ', '', meta_data);
if(is.null(overlap) & 'overlap' %in% keys) overlap = as.numeric(vals[keys=='overlap']);
if(is.null(overlap)) stop('The overlap argument is mandatory for Nonpareil versions before 2.0.');
if(is.na(read.length) & 'L' %in% keys) read.length = as.numeric(vals[keys=='L']);
if(is.na(read.length)) read.length=101;
num_reads = NULL;
if('R' %in% keys) num_reads = as.numeric(vals[keys=='R']);
if('divide' %in% keys) log_sampling = as.numeric(vals[keys=='divide']);
if('logsampling' %in% keys) log_sampling = as.numeric(vals[keys=='logsampling']);
# Read input
out <- list(kappa=0, C=0, LRstar=0, LR=0, modelR=0, diversity=0);
a <- read.table(file, sep="\t", h=F);
out$kappa <- a$V2[nrow(a)];
if(is.null(num_reads)) num_reads = max(a$V1);
LR <- exp(log(num_reads) + log(read.length));
out$LR <- LR;
for(i in 2:6) a[, i] <- a[, i]^Nonpareil.coverageFactor(overlap);
if(useValue=='median'){
values = a[, 5]
}else if(useValue=='ub'){
values = pmin(1, a[, 2] + a[, 3]*1.9)
}else if(useValue=='lb'){
values = pmax(0, a[, 2] - a[, 3]*1.9)
}else if(useValue=='q1'){
values = a[, 4]
}else if(useValue=='q3'){
values = a[, 6]
}else{
values = a[, 2]
}
a$V1 = exp(max(log(a$V1)) + max(values^0.27)*(log(a$V1) - max(log(a$V1))));
a$V1 = exp(log(a$V1)*0.61 + 10);
a$V1 = a$V1 * read.length / 101;
horiz.diff = LR / max(a$V1);
a$V1 = a$V1 * horiz.diff;##
horiz.diff = 1;
if(is.na(libname)) {
libname <- basename(file);
if(substr(libname, nchar(libname)-3, nchar(libname))==".npo")
libname <- substr(libname, 0, nchar(libname)-4);
}
# Check data
out$C <- max(values);
if(data.consistency){
twenty.pc = which.max(a$V1[a$V1<=0.2*max(a$V1)]);
if(a[twenty.pc, 5]==0){
warning('Median of the curve is zero at 20% of the reads, check parameters and re-run (e.g., decrease value of -L in nonpareil).');
return(out);
}
if(a[nrow(a), 2]<=1e-5){
warning('Curve too low, check parameters and re-run (e.g., decrease value of -L in nonpareil).');
return(out);
}
if(a[2, 2]>=1-1e-5){
warning('Curve too steep, check parameters and re-run (e.g., decrease value of -i in nonpareil or use -d).');
return(out);
}
#if(min(a[a[,2]>0, 2])<0.01){
# warning('Curve undefined around S* (kappa=05%), check parameters and re-run (e.g., decrease the value of -i in nonpareil).');
# return(0);
#}
}
# For future calls (of this or other Nonpareil.* functions)
Nonpareil.LastFactor <<- factor;
# Draw it
if(plot){
if(new){
if(!is.null(xlab)){ xlab <- as.character(xlab) }
else if(factor==1){ xlab <- 'Sequencing effort (bp)' }
else if(factor==1e-1){ xlab <- 'Sequencing effort (Dbp)' }
else if(factor==1e-2){ xlab <- 'Sequencing effort (Hbp)' }
else if(factor==1e-3){ xlab <- 'Sequencing effort (Kbp)' }
else if(factor==1e-6){ xlab <- 'Sequencing effort (Mbp)' }
else if(factor==1e-9){ xlab <- 'Sequencing effort (Gbp)' }
else if(factor==1e-12){ xlab <- 'Sequencing effort (Tbp)' }
else if(factor==1e-15){ xlab <- 'Sequencing effort (Pbp)' }
else if(factor==1e-18){ xlab <- 'Sequencing effort (Ebp)' }
else if(factor==1e-21){ xlab <- 'Sequencing effort (Zbp)' }
else{ xlab <- paste('Sequencing effort (by ', factor, 'bp)', sep='') }
if(is.null(ylab)) ylab <- 'Estimated average coverage';
plot(1, t='n', xlim=c(xmin, xmax), ylim=c(ymin, ymax), xlab=xlab, ylab=ylab, log=log, ...);
abline(h=c(1,star/100), lty=2, col='red');
abline(v=10^seq(0,15,by=3), lty=2, col='gray80')
}
if(!is.na(col)){
r <- col2rgb(col)['red',1]/256;
g <- col2rgb(col)['green',1]/256;
b <- col2rgb(col)['blue',1]/256;
}
if(is.na(r)) r <- sample(200,1);
if(is.na(g)) g <- sample(200,1);
if(is.na(b)) b <- sample(200,1);
if(r>1) r <- r/256;
if(g>1) g <- g/256;
if(b>1) b <- b/256;
if(!modelOnly){
if(!is.na(plotDispersion)){
if(plotDispersion == 'sd'){
err.y <- c(values+a$V3, rev(values-a$V3));
}else if(plotDispersion == 'ci95'){
err.y <- c(values+a$V3*1.9, rev(values-a$V3*1.9));
}else if(plotDispersion == 'ci90'){
err.y <- c(values+a$V3*1.64, rev(values-a$V3*1.64));
}else if(plotDispersion == 'ci50'){
err.y <- c(values+a$V3*.67, rev(values-a$V3*.67));
}else if(plotDispersion == 'iq'){
err.y <- c(a$V4, rev(a$V6));
}
polygon(horiz.diff*c(a$V1, rev(a$V1))*factor, ifelse(err.y<=ymin*0.1, ymin*0.1, err.y), col=rgb(r,g,b,.2), border=NA);
}
lines(horiz.diff*a$V1*factor, values, col=rgb(r,g,b,curve.alpha), lwd=curve.lwd);
}
# Save some info
Nonpareil.LastXmax <<- xmax;
Nonpareil.OpenColors <<- c(Nonpareil.OpenColors, rgb(r,g,b));
Nonpareil.OpenNames <<- c(Nonpareil.OpenNames, libname);
}
# Model it
sel <- values>0 & values<0.9;
x <- a$V1[sel];
if((length(x)>10) | ! data.consistency){
y <- values[sel];
data <- list(x=x, y=y)
estModel <- TRUE;
if(is.na(weights.exp[1])){
if(log_sampling==0){ weights.exp <- c(-1.1,-1.2,-0.9,-1.3,-1) }
else{ weights.exp <- c(0, 1, -1, 1.3, -1.1, 1.5, -1.5) }
}
weights.i <- 0;
while(estModel){
weights.i <- weights.i+1;
model <- nls(y ~ Nonpareil.f(x, a, b), data=data, weights=(a$V3[sel]^weights.exp[weights.i]), start=list(a=1, b=0.1), lower=c(a=0, b=0), algorithm='port',
control=nls.control(minFactor=1e-25000, tol=1e-15, maxiter=1024, warnOnly=T));
if(!is.na(weights.exp[weights.i+1]) | summary(model)$convInfo$isConv) estModel <- FALSE;
}
if(summary(model)$convInfo$isConv){
model.lty=2;
if(modelOnly) model.lty=1;
if(plot & plotModel){
model.x <- exp(seq(log(xmin), log(xmax), length.out=1e3));
model.y <- predict(model, list(x=model.x));
model.x <- model.x*horiz.diff*factor;
lines(model.x, model.y, col=rgb(r,g,b,model.alpha), lty=model.lty, lwd=model.lwd);
if(modelOnly) points(max(a$V1), predict(model, list(x=max(a$V1))), col=rgb(r,g,b), pch=21, bg='white');
#if(modelOnly) points(min(a$V1[a$V1>0]), predict(model, list(x=min(a$V1[a$V1>0]))), col=rgb(r,g,b), pch=8, bg='white');
}
if(returnModelValues) { out$model.x <- model.x ; out$model.y <- model.y; }
if(returnModelParameters) { out$model=model }
pa <- coef(model)[1];
pb <- coef(model)[2];
if(pa > 1) out$diversity <- (pa-1)/pb;
if(plot & plotDiversity & out$diversity>0)
arrows(x0=exp(out$diversity),
y0=ifelse(log=='y' | log=='xy' | log=='yx', ymin*0.1, ymin-1),
y1=ymin,
length=0.1, col=rgb(r,g,b,model.alpha));
out$LRstar <- Nonpareil.antif(star/100, pa, pb);
out$modelR <- cor(y, predict(model, list(x=x)));
}else{
warning('Model didn\'t converge.');
}
}else{
warning('Insufficient resolution below 90% coverage, skipping model');
}
return(out);
### A list with the following entries:
###
### kappa: "Redundancy" value of the entire dataset.
###
### C: Average coverage of the entire dataset.
###
### LRstar: Estimated sequencing effort required to reach the objective average coverage (star, 95% by default).
###
### LR: Actual sequencing effort of the dataset.
###
### modelR: Pearson's R coefficient betweeen the rarefied data and the projected model.
###
### diversity: Nonpareil-based diversity index. This value's units are the natural logarithm of the units of
### sequencing effort (by default, log-bp), and indicates the inflection point of the fitted model for the
### Nonpareil curve. If the fit doesn't converge, or the model is not estimated, the value is zero (0).
###
### model.x (if retturnModelValues=TRUE): Values of sequencing effort in the projected model.
###
### model.y (if retturnModelValues=TRUE): Values of average coverage in the projected model.
###
### model (if returnModelParameters=TRUE): Fitted model.
###
}
Nonpareil.legend <- function(
### Generates a legend for Nonpareil plots, after calling Nonpareil.curve() or Nonpareil.curve.batch().
x=NULL,
### X coordinate, or any character string accepted by legend (e.g., 'bottomright').
y=.3,
### Y coordinate.
...
### Any other parameters supported by legend().
){
if(is.null(Nonpareil.LastFactor) || is.null(Nonpareil.LastXmax))
stop('There must be at least one active plot to draw a legend. Use Nonpareil.curve().');
if(is.null(x)) x <- 0.75*Nonpareil.LastXmax;
legend(x=x, y=y, legend=Nonpareil.OpenNames, fill=Nonpareil.OpenColors, ...);
}
## PREDICTIONS
Nonpareil.coverageFactor <- function(
### Factor to transform redundancy into coverage (internal function).
overlap
### Value of overlap (-L).
){
if(overlap==25){
return(.845);
}else if(overlap==50){
return(.757);
}else if(overlap==75){
return(.633);
}else if(overlap==100){
return(.137);
}else{
stop('Unsupported overlap. Supported values are: 100, 75, 50, and 25.');
}
### A numeric scalar.
}