Fixing some indentation

Author: andrewdavidsmith

R/compared_methods.R

@@ -1,26 +1,25 @@
-#  Copyright (C) 2016 University of Southern California and
-# Chao Deng and Andrew D. Smith and Timothy Daley
+# Copyright (C) 2016-2022 University of Southern California and
+#  Chao Deng and Andrew D. Smith and Timothy Daley
 #
-#  Authors: Chao Deng
+# Authors: Chao Deng
 #
-#  This program is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
 #
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
 #
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 ## Zero Truncated Poisson (ZTP) estimator
-## Cohen, A. Clifford. "Estimating the parameter in a conditional Poisson 
+## Cohen, A. Clifford. "Estimating the parameter in a conditional Poisson
 ## distribution." Biometrics 16, no. 2 (1960): 203-211.
-ztp.rSAC <- function(n, r=1) 
+ztp.rSAC <- function(n, r=1)
 {
  ## mle based on zero-truncated Poisson
  n[, 2] <- as.numeric(n[, 2])
@@ -31,22 +30,22 @@ ztp.rSAC <- function(n, r=1)
  lambda <- result$root
  L <- sum(n[, 2]) / (1 - ppois(0, lambda))
 
- ## estimator 
+ ## estimator
  function(x) {
  L * ppois(q=r - 1, lambda=lambda * x, lower.tail=FALSE)}
 }
 
 
 ## BBC estimator
 ## Boneh, Shahar, Arnon Boneh, and Richard J. Caron. "Estimating the prediction
-## function and the number of unseen species in sampling with replacement." 
+## function and the number of unseen species in sampling with replacement."
 ## Journal of the American Statistical Association 93, no. 441 (1998): 372-379.
 bbc.rSAC <- function(n, r=1) {
  n[, 2] <- as.numeric(n[, 2])
  S <- sum(n[, 2])
  ## BBC estimator without bias correction
  tmp <- function(t) {sapply(t, function(x) {
-  n[, 2] %*% (exp(-n[, 1]) - ppois(r-1, n[, 1] * x)) + S})}
+ n[, 2] %*% (exp(-n[, 1]) - ppois(r-1, n[, 1] * x)) + S})}
  ## bias correction
  index.f1 <- which(n[, 1] == 1)
  f1 <- n[index.f1, 2]
@@ -66,8 +65,8 @@ bbc.rSAC <- function(n, r=1) {
 
 
 ## CS estimator
-## Chao, Anne, and Tsung-Jen Shen. "Nonparametric prediction in species 
-## sampling." Journal of agricultural, biological, and environmental 
+## Chao, Anne, and Tsung-Jen Shen. "Nonparametric prediction in species
+## sampling." Journal of agricultural, biological, and environmental
 ## statistics 9, no. 3 (2004): 253-269.
 cs.rSAC <- function(n, r=1, k=10) {
  n[, 2] <- as.numeric(n[, 2])
@@ -82,9 +81,9 @@ cs.rSAC <- function(n, r=1, k=10) {
  S.rare <- sum(n[index.rare, 2])
  C.rare <- 1 - f1 / (n[index.rare, 1] %*% n[index.rare, 2])
  ## estimate parameters
- gamma.rare <- max(S.rare / C.rare * 
+ gamma.rare <- max(S.rare / C.rare *
  ((n[index.rare, 1] * (n[index.rare, 1] - 1)) %*% n[index.rare, 2]) /
- (n[index.rare, 1] %*% n[index.rare, 2])^2 - 1, 0)
+  (n[index.rare, 1] %*% n[index.rare, 2])^2 - 1, 0)
  f0 <- S.rare / C.rare + f1 / C.rare * gamma.rare - S.rare
  ## estimator
  # consistent type for estimator arithmetic to avoid warnings
@@ -104,7 +103,7 @@ fisher.alpha <- function(n) {
  n[, 2] <- as.numeric(n[, 2])
  N <- n[, 1] %*% n[, 2]
  S <- sum(n[, 2])
- result <- uniroot(function(x) (exp(x) - 1) / x - N / S, c(0.001, 1e9), 
+ result <- uniroot(function(x) (exp(x) - 1) / x - N / S, c(0.001, 1e9),
  tol=0.0001, extendInt="upX")
  alpha <- S / result$root
  return(alpha)
@@ -115,20 +114,22 @@ fisher.rSAC <- function(n, r=1) {
  n[, 2] <- as.numeric(n[, 2])
  alpha <- fisher.alpha(n)
  N <- n[, 1] %*% n[, 2]
- f.rSAC <- function(t) {sapply(t, function(x) alpha *
- integrate(function(z) (z^(r-1) / (1 - z)), lower=0, 
- upper=N*x / (N*x + alpha))$value)}
+ f.rSAC <- function(t) {
+ sapply(t, function(x) alpha *
+ integrate(function(z) (z^(r-1) / (1 - z)), lower=0,
+ upper=N*x / (N*x + alpha))$value)
+ }
  return(f.rSAC)
 }
 
 
 ## fit a Poisson distribution
 pois.rSAC <- function(n, L, r=1) {
-  lambda <- n[, 1] %*% n[, 2] / L
-  f.rSAC <- function(t) {
-  L * ppois(q=r - 1, lambda=lambda * t, lower.tail=FALSE)
-  }
-  return(f.rSAC)
+ lambda <- n[, 1] %*% n[, 2] / L
+ f.rSAC <- function(t) {
+ L * ppois(q=r - 1, lambda=lambda * t, lower.tail=FALSE)
+ }
+ return(f.rSAC)
 }
 
 ## estimate the parameter for negative binomial distribuiton
@@ -145,14 +146,14 @@ nb.fitting <- function(n, L, size=SIZE.INIT)
 
  ## target function f
  f <- function(x) {
-  return( -nb.loglikelihood(n, zero.counts, size = x, mu = m)/L )
+ return( -nb.loglikelihood(n, zero.counts, size = x, mu = m)/L )
  }
 
  ## derivative of f
  gr <- function(x)
  {
  first.term <- ( digamma(x) * zero.counts +
- digamma(n[, 1] + x) %*% n[, 2] )/L
+  digamma(n[, 1] + x) %*% n[, 2] )/L
  second.term <- digamma(x)
  third.term <- log(x) - log(x + m)
  result <- first.term - second.term + third.term
@@ -163,10 +164,10 @@ nb.fitting <- function(n, L, size=SIZE.INIT)
  ## optimization
  if (v > m) {
  res <- optim(m^2 / (v - m), f, gr, method = "L-BFGS-B",
- lower = 0.00001, upper = 100000)
+  lower = 0.00001, upper = 100000)
  } else {
  res <- optim(size, f, gr, method = "L-BFGS-B",
- lower = 0.00001, upper = 100000)
+  lower = 0.00001, upper = 100000)
  }
 
  loglikelihood <- nb.loglikelihood(n, zero.counts, size=res$par, mu=m)

R/kmer.R

@@ -1,30 +1,28 @@
-#  Copyright (C) 2016 University of Southern California and
-# Chao Deng and Andrew D. Smith and Timothy Daley
+# Copyright (C) 2016-2022 University of Southern California and
+#  Chao Deng and Andrew D. Smith and Timothy Daley
 #
-#  Authors: Chao Deng
+# Authors: Chao Deng
 #
-#  This program is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
 #
-#  This program is distributed in the hope that it will be useful,
-#  but WITHOUT ANY WARRANTY; without even the implied warranty of
-#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-#  GNU General Public License for more details.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
 #
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-#
-
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 ## predict the fraction of k-mers represented at least r times in the sample
 kmer.frac <- function(n, r=2, mt=20) {
  return(preseqR.sample.cov(n=n, r=r-1, mt=mt))
 }
 
 
-## the fraction of k-mers represented at least r times as a function of 
+## the fraction of k-mers represented at least r times as a function of
 ## sample sizes
 kmer.frac.curve <- function(n, k, read.len, seq, r=2, mt=20) {
  f <- kmer.frac(n, r=r, mt=mt)
@@ -50,7 +48,7 @@ kmer.frac.bootstrap <- function(n, r=2, mt=20, times=30, conf=0.95) {
 }
 
 
-## the fraction of k-mers represented at least r times as a function of 
+## the fraction of k-mers represented at least r times as a function of
 ## sample sizes
 kmer.frac.curve.bootstrap <- function(n, k, read.len, seq, r=2, mt=20,
  times=30, conf=0.95)
@@ -66,9 +64,9 @@ kmer.frac.curve.bootstrap <- function(n, k, read.len, seq, r=2, mt=20,
  # consistent vector-vector arithmetic
  unit <- as.numeric(unit)
  seq.effort <- seq / unit
- result <- matrix(c(seq, f$f(seq.effort), f$lb(seq.effort), 
+ result <- matrix(c(seq, f$f(seq.effort), f$lb(seq.effort),
  f$ub(seq.effort)), ncol=4, byrow=FALSE)
- colnames(result) <- c("bases", paste("frac(X>=", r, ")", sep=""), 
+ colnames(result) <- c("bases", paste("frac(X>=", r, ")", sep=""),
  "lb", "ub")
  return(result)
 }

R/rSAC.R

@@ -1,23 +1,22 @@
-#  Copyright (C) 2016 University of Southern California and
+# Copyright (C) 2016-2022 University of Southern California and
 # Chao Deng and Andrew D. Smith and Timothy Daley
 #
-#  Authors: Chao Deng
+# Authors: Chao Deng
 #
-#  This program is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
 #
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
 #
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>
 
-## the pair, a pole and a root, is a defect if the distance is less than 
+## the pair, a pole and a root, is a defect if the distance is less than
 ## the value defined by the variable PRECISION
 PRECISION <- 1e-3
 
@@ -62,8 +61,8 @@ preseqR.interpolate.rSAC <- function(n, ss, r=1)
  if (x <= r - 1) {
  return(1)
  } else {
-# consistent vector-vector arithmetic
-denom <- as.numeric(denom)
+ # consistent vector-vector arithmetic
+ denom <- as.numeric(denom)
  logp = lchoose(N - x, size - 0:(r-1)) + lchoose(x, 0:(r-1)) - denom
  return(sum(exp(logp)))
  }})
@@ -93,7 +92,7 @@ preseqR.interpolate.rSAC <- function(n, ss, r=1)
 # frequencies <- lapply(1:i, function(x) rbinom(n[x, 2], n[x, 1], p))
 # f <- unlist(frequencies)
 # h <- hist(f, breaks=-1:max(f), plot=FALSE)$counts[-1]
-# matrix(c(which(h != 0), h[which(h != 0)]), byrow = FALSE, ncol=2) 
+# matrix(c(which(h != 0), h[which(h != 0)]), byrow = FALSE, ncol=2)
 #}
 
 
@@ -113,7 +112,7 @@ discoveryrate.ps <- function(n, mt)
  change.sign <- 0
  for (i in 1:(min(mt-1, length(hist.count)))) {
  PS.coeffs <- c(
- PS.coeffs, 
+ PS.coeffs,
  (-1)^change.sign * hist.count[i] - PS.coeffs[length(PS.coeffs)])
  change.sign <- change.sign + 1
  }
@@ -250,7 +249,7 @@ ds.rSAC.bootstrap <- function(n, r=1, mt=20, times=30, conf=0.95)
  return(apply(result, FUN=var, MARGIN=1))
  }
  }
- 
+
  se <- function(x) sqrt(variance(x))
 
  ## prevent later binding!!!
@@ -298,7 +297,7 @@ preseqR.rSAC <- function(n, r=1, mt=20, size=SIZE.INIT, mu=MU.INIT)
 
 ## the bootstrap version of preseqR.rSAC
 ## with confidence interval
-preseqR.rSAC.bootstrap <- function(n, r=1, mt=20, size=SIZE.INIT, 
+preseqR.rSAC.bootstrap <- function(n, r=1, mt=20, size=SIZE.INIT,
  mu=MU.INIT, times=30, conf=0.95)
 {
  checking.hist(n)