diff --git a/NAMESPACE b/NAMESPACE
index c56243b..b30de66 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,20 +1,31 @@
 # Generated by roxygen2: do not edit by hand
 
 export(aldex2propr)
-export(basis_shrinkage)
 export(calculate_theta)
 export(compute_iqlr)
+export(getAdjacencyFDR)
+export(getAdjacencyFDR.propd)
+export(getAdjacencyFDR.propr)
+export(getAdjacencyFstat)
+export(getCutoffFDR)
+export(getCutoffFDR.propd)
+export(getCutoffFDR.propr)
+export(getCutoffFstat)
 export(getMatrix)
 export(getRatios)
 export(getResults)
+export(getSignificantResultsFDR)
+export(getSignificantResultsFDR.propd)
+export(getSignificantResultsFDR.propr)
+export(getSignificantResultsFstat)
 export(index_reference)
 export(logratio)
 export(logratio_with_alpha)
 export(logratio_without_alpha)
+export(pcor.bshrink)
 export(propd)
 export(propr)
 export(ratios)
-export(runCutoff)
 export(runGraflex)
 export(runNormalization)
 export(runPostHoc)
diff --git a/NEWS.md b/NEWS.md
index a5efd80..95bf39a 100755
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,16 @@
+## propr 5.1.0
+---------------------
+* Allowed `updateCutoffs` function to compute the FDR values for negative cutoffs
+* Added `getCutoffFDR` to get a significant cutoff based on the permuted FDR
+* Changed `runCutoff` into `getCutoffFstat`
+* Added `getSignificantResultsFDR` and `getSignificantResultsFstat` to get the significant pairs as results table
+* Added `getAdjacencyFDR` and `getAdjacencyFstat` functions to get the adjacency matrix with the significant pairs
+* Optimized `graflex` related functions to Rcpp C++
+* Removed fixseed, since setting `set.seed()` before calling any function is enough to ensure reproducibility
+* Fix bug: when `ivar=NA`, the already preprocessed data given by the user are replaced when there are zeros, which should not
+* Replaced `ppcor::pcor` by `corpcor::cor2pcor` for coherence with the shrunk partial correlations
+* Parallelized `propd` when ncores > 1
+
 ## propr 5.0.4
 ---------------------
 * Added fixseed option for runGraflex
diff --git a/R/1-propr-OOP.R b/R/1-propr-OOP.R
index 6f15110..2bc642a 100644
--- a/R/1-propr-OOP.R
+++ b/R/1-propr-OOP.R
@@ -6,7 +6,10 @@ setClass(
     counts = "data.frame",
     alpha = "numeric",
     metric = "character",
+    direct = "logical",
+    has_meaningful_negative_values = "logical",
     ivar = "ANY",
+    lambda = "ANY",
     logratio = "data.frame",
     matrix = "matrix",
     pairs = "numeric",
diff --git a/R/1-propr.R b/R/1-propr.R
index 915ecae..bb283ad 100755
--- a/R/1-propr.R
+++ b/R/1-propr.R
@@ -18,8 +18,9 @@
 #'   - "clr" (default): Centered log-ratio transformation.
 #'   - "alr": Additive log-ratio transformation ("pcor.bshrink" metric only).
 #'   - "iqlr": Inter-quartile log-ratio transformation from ALDEx2.
-#'   - The explicit name(s) of variable(s) to use as a reference.
-#'  Use NA to skip log-ratio transformation.
+#'   - The explicit name(s) or index(es) of variable(s) to use as a reference.
+#'   - Use NA to skip log-ratio transformation and any other pre-processing, like
+#'  zero replacement. This is useful when the input data is already pre-processed.
 #' @param select A numeric vector representing the indices of features to be
 #'  used for computing the Propr matrix. This argument is optional. If
 #'  provided, it reduces the data size by using only the selected features.
@@ -29,9 +30,6 @@
 #' @param alpha The alpha parameter used in the alpha log-ratio transformation.
 #' @param p The number of permutations to perform for calculating the false
 #'  discovery rate (FDR). The default is 0.
-#' @param fixseed A logical value indicating whether to fix the random seed 
-#' for generating the pemurations. If `TRUE`, the function will fix the random
-#' The default is `FALSE`.
 #' @param ... Additional arguments passed to \code{corpcor::pcor.shrink},
 #'  if "pcor.shrink" metric is selected.
 #'
@@ -62,6 +60,7 @@ propr <- function(counts,
                              "phs",
                              "cor",
                              "vlr",
+                             "ppcor",
                              "pcor",
                              "pcor.shrink",
                              "pcor.bshrink"),
@@ -70,7 +69,6 @@ propr <- function(counts,
                   symmetrize = FALSE,
                   alpha = NA,
                   p = 0,
-                  fixseed = FALSE,
                   ...) {
   ##############################################################################
   ### CLEAN UP ARGS
@@ -82,24 +80,37 @@ propr <- function(counts,
 
   # Special handling for 'metric'
   metric <- metric[1]
+  ivar_pcor <- NA
   if (metric == "pcor.bshrink") {
+    if (length(ivar) != 1) {
+      stop("For 'pcor.bshrink', the 'ivar' argument must be a single value.")
+    }
     if (!ivar %in% c("clr", "alr")) {
       stop("For 'pcor.bshrink', the 'ivar' argument must be 'clr' or 'alr'.")
     }
     message("Alert: Log-ratio transform will be handled by 'bShrink'.")
     ivar_pcor <- ivar
-    ivar <-
-      NA # skips log-ratio transform as performed for all other metrics
+    ivar <- NA # skips log-ratio transform
+  } else {
+    if (length(ivar) == 1 && !is.na(ivar) && ivar == "alr") {
+      stop("Please give the index or name of the reference gene instead of 'alr'.")
+    } else if (length(ivar) > 1 && any(c("alr","clr","iqlr") %in% ivar)) {
+      stop("Please check the ivar argument is correct.")
+    }
   }
 
   ##############################################################################
   ### PERFORM ZERO REPLACEMENT AND LOG-RATIO TRANSFORM
   ##############################################################################
 
-  # NOTE: ct will never have zeros; counts may or may not have zeros!
+  # NOTE: counts are the original counts, while ct may have zeros replaced
   counts <- as_safe_matrix(counts)
-  ct <- simple_zero_replacement(counts)
-  lr <- logratio(counts, ivar, alpha)
+  if (length(ivar) == 1 && is.na(ivar) && is.na(ivar_pcor)) {
+    ct <- counts
+  } else {
+    ct <- simple_zero_replacement(counts)
+  }
+  lr <- logratio(ct, ivar, alpha)
 
   ##############################################################################
   ### OPTIONALLY REDUCE DATA SIZE BEFORE COMPUTING MATRIX
@@ -117,28 +128,50 @@ propr <- function(counts,
   ##############################################################################
 
   lrv <- lr2vlr(lr)
+  lambda <- NULL
+
   if (metric == "rho") {
     mat <- lr2rho(lr)
+
   } else if (metric == "phi") {
     mat <- lr2phi(lr)
     if (symmetrize)
       symRcpp(mat) # optionally force symmetry
+
   } else if (metric == "phs") {
     mat <- lr2phs(lr)
+
   } else if (metric == "cor") {
     mat <- stats::cor(lr)
+
   } else if (metric == "vlr") {
     mat <- lrv
-  } else if (metric == "pcor") {
+
+  } else if (metric == "ppcor") {
     packageCheck("ppcor")
     mat <- ppcor::pcor(lr)$estimate
+
+  } else if (metric == "pcor") {
+    packageCheck("corpcor")
+    cov <- cov(lr)
+    mat <- corpcor::cor2pcor(cov)
+    class(mat) <- "matrix"
+
   } else if (metric == "pcor.shrink") {
     packageCheck("corpcor")
-    mat <- corpcor::pcor.shrink(lr, ...)
+    cov <- corpcor::cov.shrink(lr)
+    mat <- corpcor::cor2pcor(cov)
+    mat <- matrix(mat, ncol=ncol(lr), nrow=ncol(lr))
     class(mat) <- "matrix"
+    lambda <- attr(cov, "lambda")
+
   } else if (metric == "pcor.bshrink") {
-    mat <- basis_shrinkage(ct, outtype = ivar_pcor)
-  } else{
+    with(pcor.bshrink(ct, outtype = ivar_pcor), {
+      mat <<- matrix
+      lambda <<- lambda
+    })
+
+  } else {
     stop("Provided 'metric' not recognized.")
   }
 
@@ -155,8 +188,11 @@ propr <- function(counts,
   result@logratio <- as.data.frame(lr)
   result@pairs <- vector("numeric")
   result@permutes <- list(NULL)
+  result@lambda <- lambda
+  result@direct <- ifelse(metric[1] %in% c("rho", "cor", "pcor", "pcor.shrink", "pcor.bshrink"), TRUE, FALSE)
+  result@has_meaningful_negative_values <- ifelse(metric[1] %in% c("cor", "pcor", "pcor.shrink", "pcor.bshrink"), TRUE, FALSE)  # metrics like proportionality has negative values that are difficult to interpret, whereas correlation metrics have a clear interpretation
 
-  # ivar should not be NA, otherwise updateCutoffs does not work
+  # ivar should not be NA for pcor.bshrink, otherwise updateCutoffs does not work
   if (metric == 'pcor.bshrink') result@ivar <- ivar_pcor
 
   # Clean row and column names
@@ -178,7 +214,7 @@ propr <- function(counts,
     )
 
   # permute data
-  if (p > 0) result <- updatePermutes(result, p, fixseed)
+  if (p > 0) result <- updatePermutes(result, p)
 
   ##############################################################################
   ### GIVE HELPFUL MESSAGES TO USER
diff --git a/R/1a-propr-backend.R b/R/1a-propr-backend.R
index 9c776da..c80a342 100644
--- a/R/1a-propr-backend.R
+++ b/R/1a-propr-backend.R
@@ -44,6 +44,12 @@ index_reference <- function(counts, ivar) {
     use <- sort(ivar) # use features given by number
   }
 
+  if (length(use) == 0) {
+      stop("No reference detected. Please provide a valid 'ivar' argument.")
+  } else if(length(use) > ncol(counts)) {
+    stop("Detected more reference variables than existing in the data. Please provide a valid 'ivar' argument.")
+  }
+
   return(use)
 }
 
@@ -76,30 +82,6 @@ compute_iqlr <- function(counts) {
   return(use)
 }
 
-#' Simple Zero Replacement in a Count Matrix
-#'
-#' This function replaces zeros with the next smallest non-zero value in the
-#'  input count matrix. If the matrix contains no zeros, it produces an
-#'  informational message indicating that no replacements were made.
-#'
-#' @inheritParams logratio_with_alpha
-#' @return A matrix with zero values replaced by the next smallest non-zero value.
-#' @examples
-#' # Sample input count data with zeros
-#' data <- matrix(c(0, 2, 3, 4, 5, 0), nrow = 2, byrow = TRUE)
-#'
-#' @export
-simple_zero_replacement <- function(ct) {
-  if (any(ct == 0)) {
-    zeros <- ct == 0
-    ct[zeros] <- min(ct[!zeros])
-  } else{
-    message("Alert: No 0s found that need replacement.")
-  }
-
-  return(ct)
-}
-
 #' Perform log-ratio transformation without alpha parameter
 #'
 #' This function applies a log-ratio transformation to a given data matrix
@@ -112,11 +94,14 @@ simple_zero_replacement <- function(ct) {
 #' # Sample input data
 #' data <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, byrow = TRUE)
 #'
-#' # Applying log-ratio transformation to rows using columns 2 and 3
+#' # Applying log-ratio transformation to rows using columns 2 and 3 as reference
 #' result <- logratio_without_alpha(data, use = c(2, 3))
 #'
 #' @export
 logratio_without_alpha <- function(ct, use) {
+  if (any(ct == 0)){
+    stop("Error: please first replace the zeros before running logratio.")
+  }
   # Use g(x) = Mean[log(x)] to log-ratio transform data
   logX <- log(ct)
   logSet <- logX[, use, drop = FALSE]
@@ -186,15 +171,13 @@ logratio_with_alpha <- function(ct, use, alpha) {
 #' @export
 logratio <- function(counts, ivar, alpha) {
   counts <- as_safe_matrix(counts)
-  if (is.na(ivar)) {
+  if (length(ivar) == 1 && is.na(ivar)) {
     message("Alert: Skipping built-in log-ratio transformation.")
-    ct <- simple_zero_replacement(counts)
-    lr <- ct
+    lr <- counts
   } else {
     use <- index_reference(counts, ivar) # Index reference
     if (is.na(alpha)) {
-      ct <- simple_zero_replacement(counts)
-      lr <- logratio_without_alpha(ct, use)
+      lr <- logratio_without_alpha(counts, use)
     } else {
       lr <- logratio_with_alpha(counts, use, alpha)
     }
@@ -202,59 +185,66 @@ logratio <- function(counts, ivar, alpha) {
   return(lr)
 }
 
-#' Covariance Shrinkage and Partial Correlation Calculation
+#' Basis Covariance Shrinkage and Partial Correlation Calculation
 #'
-#' This function performs covariance shrinkage and calculates the partial
-#'  correlation matrix based on the input count data matrix. The function can
+#' This function performs covariance shrinkage on the basis matrix
+#'  and calculates the partial correlation matrix. The function can
 #'  output the results in two formats: centered log-ratio (clr) or
 #'  additive log-ratio (alr).
 #'
-#' @param M A data matrix representing the counts. Each row should represent
+#' @param ct A data matrix representing the counts. Each row should represent
 #'  observations, and each column should represent different variables.
 #' @param outtype A character vector specifying the output type. It can take
-#'  either "clr" (centered log-ratio) or "alr" (additive log-ratio).
+#'  either "clr" (centered log-ratio) or "alr" (additive log-ratio). Since
+#'  the reference does not affect the logratio partial correlation coefficients,
+#'  the index of the reference is not needed for "alr". Also, "clr" is recommended
+#'  because of the same reason, at the same time avoiding losing one dimension.
 #'
-#' @return A matrix representing the partial correlation matrix.
+#' @return A matrix representing the shrunk partial correlation matrix.
 #'
 #' @examples
 #' # Sample input count data
 #' data <- iris[,1:4]
 #'
 #' # Calculate partial correlation matrix using clr transformation
-#' result_clr <- basis_shrinkage(data, outtype = "clr")
+#' result_clr <- pcor.bshrink(data, outtype = "clr")
 #'
 #' # Calculate partial correlation matrix using alr transformation
-#' result_alr <- basis_shrinkage(data, outtype = "alr")
+#' result_alr <- pcor.bshrink(data, outtype = "alr")
 #'
 #' @export
-basis_shrinkage <- function(M, outtype = c("clr", "alr")) {
+pcor.bshrink <- function(ct, outtype = c("clr", "alr")) {
   packageCheck("corpcor")
   outtype <- match.arg(outtype)
 
   # transform counts to log proportions
-  P <- M / rowSums(M)
-  B <- log(P)
+  logP <- log( ct / rowSums(ct) )
 
   # covariance shrinkage
-  D  <- ncol(M)
-  Cb <- corpcor::cov.shrink(B, verbose = FALSE)
+  covB <- corpcor::cov.shrink(logP, verbose = FALSE)
+  lambda <- attr(covB, "lambda")
+
+  # convert basis covariance matrix to clr/alr covariance matrix
+  D  <- ncol(ct)
   if (outtype == "alr") {
     F   <- cbind(diag(rep(1, D - 1)), rep(-1, D - 1))
-    Cov <- F %*% Cb %*% t(F)
+    cov <- F %*% covB %*% t(F)
   } else if (outtype == "clr") {
     G   <- diag(rep(1, D)) - matrix(1 / D, D, D)
-    Cov <- G %*% Cb %*% G
+    cov <- G %*% covB %*% G
   }
 
   # partial correlation
-  PC <- corpcor::cor2pcor(Cov)
+  pcor <- corpcor::cor2pcor(cov)
 
   # make output to have same dimensions as input
+  # alr partial correlation has one less dimension,
+  # so here we add a row and a column of 0s
   if (outtype == "alr") {
-    PC <- cbind(PC, replicate(nrow(PC), 0))
-    PC <- rbind(PC, replicate(ncol(PC), 0))
-    PC[nrow(PC), ncol(PC)] <- 1
+    pcor <- cbind(pcor, 0)
+    pcor <- rbind(pcor, 0)
+    pcor[ncol(pcor), ncol(pcor)] <- 1
   }
 
-  return(PC)
+  return(list(matrix = pcor, lambda = lambda))
 }
diff --git a/R/2-propd.R b/R/2-propd.R
index cb75129..3aaeab4 100644
--- a/R/2-propd.R
+++ b/R/2-propd.R
@@ -5,9 +5,6 @@
 #' @param alpha The alpha parameter used in the alpha log-ratio transformation.
 #' @param p The number of permutations to perform for calculating the false
 #'  discovery rate (FDR). The default is 0.
-#' @param fixseed A logical value indicating whether to fix the random seed 
-#' for generating the pemurations. If `TRUE`, the function will fix the random
-#' The default is `FALSE`.
 #' @param weighted A logical value indicating whether weighted calculations
 #'  should be performed. If \code{TRUE}, the function will use limma-based
 #'  weights for the calculations.
@@ -36,7 +33,6 @@ propd <- function(counts,
                   group,
                   alpha = NA,
                   p = 0,
-                  fixseed = FALSE,
                   weighted = FALSE) {
   ##############################################################################
   ### CLEAN UP ARGS
@@ -108,7 +104,7 @@ propd <- function(counts,
     round(result@results$theta, 14) # round floats to 1
 
   # permute data
-  if (p > 0) result <- updatePermutes(result, p, fixseed)
+  if (p > 0) result <- updatePermutes(result, p)
 
   ##############################################################################
   ### GIVE HELPFUL MESSAGES TO USER
diff --git a/R/2b-propd-frontend.R b/R/2b-propd-frontend.R
index 5df15e3..b0e9f5f 100644
--- a/R/2b-propd-frontend.R
+++ b/R/2b-propd-frontend.R
@@ -84,12 +84,11 @@ updateF <- function(propd,
     # Establish data with regard to a reference Z
     if (any(propd@counts == 0)) {
       message("Alert: Building reference set with ivar and counts offset by 1.")
-      X <- as.matrix(propd@counts + 1)
+      X <- as.matrix(propd@counts + 1)  # TODO this should be checked, if we still want to do this
     } else{
       message("Alert: Building reference set with ivar and counts.")
       X <- as.matrix(propd@counts)
     }
-
     logX <- log(X)
     z.set <- logX[, use, drop = FALSE]
     z.geo <- rowMeans(z.set)
diff --git a/R/2c-propd-experimental.R b/R/2c-propd-experimental.R
index a403ac7..301eaed 100644
--- a/R/2c-propd-experimental.R
+++ b/R/2c-propd-experimental.R
@@ -1,54 +1,3 @@
-#' Calculate a theta Cutoff
-#'
-#' This function uses the F distribution to calculate a cutoff of
-#'  theta for a p-value given by the \code{pval} argument.
-#'
-#' If the argument \code{fdr = TRUE}, this function returns the
-#'  empiric cutoff that corresponds to the FDR-adjusted p-value
-#'  stored in the \code{@@results$FDR} slot.
-#'
-#' @param object A \code{\link{propd}} object.
-#' @param pval A p-value at which to calculate a theta cutoff.
-#' @param fdr A boolean. Toggles whether to calculate the theta
-#'  cutoff for an FDR-adjusted p-value.
-#' @return A cutoff of theta from [0, 1].
-#' @export
-runCutoff <- function(object, pval = 0.05, fdr = FALSE) {
-  if (!"Fstat" %in% colnames(object@results)) {
-    stop("Please run updateF() on propd object before calling qtheta.")
-  }
-
-  if (pval < 0 | pval > 1) {
-    stop("Provide a p-value cutoff from [0, 1].")
-  }
-
-  if (fdr) {
-    message("Alert: Returning an empiric cutoff based on the $FDR slot.")
-    index <- object@results$FDR < pval
-    if (any(index)) {
-      cutoff <- max(object@results$theta[index])
-    } else{
-      stop("No pairs below p-value.")
-    }
-
-  } else{
-    # Compute based on theory
-    K <- length(unique(object@group))
-    N <-
-      length(object@group) + object@dfz # population-level metric (i.e., N)
-    Q <- stats::qf(pval, K - 1, N - K, lower.tail = FALSE)
-    # # Fstat <- (N - 2) * (1 - object@theta$theta) / object@theta$theta
-    # # Q = Fstat
-    # # Q = (N-2) * (1-theta) / theta
-    # # Q / (N-2) = (1/theta) - 1
-    # # 1/theta = Q / (N-2) + 1 = Q(N-2)/(N-2)
-    # # theta = (N-2)/(Q+(N-2))
-    cutoff <- (N - 2) / (Q + (N - 2))
-  }
-
-  return(cutoff)
-}
-
 #' Get Per-Feature Theta
 #'
 #' This function calculates the differential proportionality
diff --git a/R/2d-propd-graflex.R b/R/2d-propd-graflex.R
deleted file mode 100644
index 68d6a67..0000000
--- a/R/2d-propd-graflex.R
+++ /dev/null
@@ -1,134 +0,0 @@
-#' Permute Odds Ratio
-#'
-#' This function permutes \code{p} odds ratios for the
-#'  edge overlap between two non-random graphs.
-#'  It does this by randomly shuffling the rows and
-#'  columns of \code{A} (jointly, thus preserving
-#'  the degree distribution).
-#'
-#' Note that this function calculates overlap for the
-#'  lower-left triangle of the input matrices.
-#' @param A,G An adjacency matrix.
-#' @param p An integer. The number of overlaps to permute.
-#' @param fixseed A logical. If \code{TRUE}, the seed is fixed
-#' for reproducibility.
-permuteOR <- function(A, G, p = 500, fixseed=FALSE) {
-  Gstar <- G[lower.tri(G)]
-  res <- lapply(1:p, function(i) {
-    # Shuffle the adjacency matrix
-    if (fixseed) set.seed(i)
-    index <- sample(1:ncol(A))
-    A <- A[index, index]
-    Astar <- A[lower.tri(A)]
-    getOR(Astar, Gstar)
-  })
-
-  do.call("rbind", res)
-}
-
-#' Tabulate Overlap
-#'
-#' This function tabulates the overlap between
-#'  two vectors or two adjacency matrices.
-#'  It is a faster version of \code{table} that
-#'  only supports binary input.
-#' @param A,G A vector or adjacency matrix.
-#' @return A table of overlap.
-binTab <- function(A, G) {
-  diff <- A != G
-  only1 <- A[diff]
-  b <- as.numeric(sum(only1))
-  c <- as.numeric(length(only1) - b)
-
-  same <- !diff
-  double1 <- A[same]
-  a <- as.numeric(sum(double1))
-  d <- as.numeric(length(double1) - a)
-
-  matrix(c(d, b, c, a), 2, 2)
-}
-
-#' Calculate Odds Ratio
-#'
-#' This function calculates the overlap between
-#'  two vectors or two adjacency matrices.
-#'  It returns the OR as well as other metrics.
-#' @return A \code{data.frame} of results.
-getOR <- function(A, G) {
-  tab <- binTab(A, G)
-  or <- (tab[1, 1] * tab[2, 2]) / (tab[1, 2] * tab[2, 1])
-  data.frame(
-    "Neither" = tab[1, 1],
-    "G.only" = tab[1, 2],
-    "A.only" = tab[2, 1],
-    "Both" = tab[2, 2],
-    "Odds" = or,
-    "LogOR" = log(or)
-  )
-}
-
-#' Calculate Odds Ratio
-#'
-#' This function calculates an odds ratio for the
-#'  edge overlap between two non-random graphs.
-#'
-#' Note that this function calculates overlap for the
-#'  lower-left triangle of the input matrices.
-calculateOR <- function(A, G) {
-  Astar <- A[lower.tri(A)]
-  Gstar <- G[lower.tri(G)]
-  getOR(Astar, Gstar)
-}
-
-#' Calculate Odds Ratio FDR
-#'
-#' This function calculates the false discovery rate (FDR)
-#'  for over- and under-enrichment by counting the number of
-#'  times the actual OR was greater than
-#'  (or less than) a permuted OR.
-#' @param actual A result from \code{\link{calculateOR}}.
-#' @param permuted A result from \code{\link{permuteOR}}.
-#' @return A \code{data.frame} of the FDRs for over-
-#'  and under- enrichment.
-getFDR <- function(actual, permuted) {
-  actual$FDR.under <-
-    sum(permuted$Odds <= actual$Odds) / nrow(permuted)
-  actual$FDR.over <-
-    sum(permuted$Odds >= actual$Odds) / nrow(permuted)
-  actual
-}
-
-#' Permute FDR for Multiple Concepts
-#'
-#' This function calls \code{\link{permuteOR}} for each
-#'  concept (i.e., column) in the database \code{K}.
-#'
-#' For each concept, this function calculates the
-#'  false discovery rate (FDR) by counting the number of
-#'  times the actual OR was greater than
-#'  (or less than) a permuted OR.
-#' @param A An adjacency matrix.
-#' @param K A knowledge database where each row is a graph node
-#'  and each column is a concept.
-#' @param p An integer. The number of overlaps to permute.
-#' @param fixseed A logical. If \code{TRUE}, the seed is fixed
-#' for reproducibility.
-#' @export
-runGraflex <- function(A, K, p = 500, fixseed=FALSE) {
-  if (nrow(A) != nrow(K))
-    stop("'A' and 'K' must have identical rows.")
-
-  numTicks <- 0
-  res <- lapply(1:ncol(K), function(k) {
-    numTicks <<- progress(k, ncol(K), numTicks)
-    Gk <- K[, k] %*% t(K[, k])
-    actual <- calculateOR(A, Gk)
-    permuted <- permuteOR(A, Gk, p = p, fixseed=fixseed)
-    actual <- getFDR(actual, permuted)
-    actual$Permutes <- p
-    actual$Concept <- colnames(K)[k]
-    actual
-  })
-
-  do.call("rbind", res)
-}
diff --git a/R/3-shared-getAdjacency.R b/R/3-shared-getAdjacency.R
new file mode 100644
index 0000000..edd01a6
--- /dev/null
+++ b/R/3-shared-getAdjacency.R
@@ -0,0 +1,135 @@
+#' Get Adjacency Matrix as indicated by permutation tests
+#'
+#' This function gets the significant pairs, according to the permutation tests. 
+#' Then it fills the adjacency matrix with 1 if pair is significant, otherwise 0.
+#' @param object A \code{propd} or \code{propr} object.
+#' @param fdr A float value for the false discovery rate. Default is 0.05.
+#' @param window_size An integer. Default is 1. When it is greater than 1, the FDR
+#' values would be smoothed out by a moving average of the given window size.
+#' @param tails 1 for one-sided on the right, and 2 for two-sided. When NULL, use default
+#' test for the given metric. This is only relevant for \code{propr} objects.
+#' @return An adjacency matrix.
+#'
+#' @export
+getAdjacencyFDR <- function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+  
+  if (inherits(object, "propr")){
+    adj <- getAdjacencyFDR.propr(object, fdr=fdr, window_size=window_size, tails=tails)
+
+  }else if(inherits(object, "propd")){
+    adj <- getAdjacencyFDR.propd(object, fdr=fdr, window_size=window_size)
+
+  }else{
+    stop("Please provide a 'propr' or 'propd' object.")
+  }
+
+  return(adj)
+}
+
+#' @rdname getAdjacencyFDR
+#' @section Methods:
+#' \code{getAdjacencyFDR.propd:}
+#' This function creates an adjacency matrix with the significant
+#' pairs from a \code{propd} object.
+#' @export
+getAdjacencyFDR.propd <- 
+  function(object, fdr = 0.05, window_size = 1) {
+
+    # get cutoff
+    cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size)
+
+    # get theta matrix
+    mat <- getMatrix(object)
+
+    # create empty matrix
+    adj <- matrix(0, nrow = nrow(mat), ncol = ncol(mat))
+    rownames(adj) <- rownames(mat)
+    colnames(adj) <- colnames(mat)
+    diag(adj) <- 1
+
+    # fill in significant pairs
+    if (cutoff) adj[mat <= cutoff] <- 1
+
+    return(adj)
+  }
+
+#' @rdname getAdjacencyFDR
+#' @section Methods:
+#' \code{getAdjacencyFDR.propr:}
+#' This function creates an adjacency matrix with the significant
+#' pairs from a \code{propd} object.
+#' @export
+getAdjacencyFDR.propr <-
+  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+
+    if (is.null(tails)) {
+      tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
+    }
+    if (tails != 1 & tails != 2) {
+        stop("Please provide a valid value for tails: 1 or 2.")
+    }
+    if (tails == 1 & object@has_meaningful_negative_values) {
+      warning("Significant pairs are chosen based on one-sided FDR test.")
+    }
+    if (tails == 2 & !object@direct) {
+      stop("Two-sided FDR is not available for this metric.")
+    }
+
+    # correct fdr when two-sided
+    if (tails == 2) fdr <- fdr / 2
+
+    # create empty matrix
+    adj <- matrix(0, nrow = ncol(object@matrix), ncol = ncol(object@matrix))
+    rownames(adj) <- rownames(object@matrix)
+    colnames(adj) <- colnames(object@matrix)
+    diag(adj) <- 1
+
+    # for positive tail
+    cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size)
+    if (cutoff) {
+      if (object@direct) {
+        adj[object@matrix >= 0 & object@matrix >= cutoff] <- 1
+      } else {
+        adj[object@matrix <= cutoff] <- 1
+      }
+    }
+
+    # for negative tail
+    if (tails == 2){
+      cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size, positive=F)
+      if (cutoff) adj[object@matrix < 0 & object@matrix <= cutoff] <- 1
+    }
+
+    return(adj)
+  }
+
+#' Get Adjacency Matrix as indicated by F-statistics
+#'
+#' This function gets the significant pairs, according to the F-statistics. 
+#' Then it fills the adjacency matrix with 1 if pair is significant, otherwise 0.
+#' Note that it can only be applied to theta_d, as updateF only works for theta_d.
+#' 
+#' @param object A \code{propd} or \code{propr} object.
+#' @param pval A float value for the p-value. Default is 0.05.
+#' @param fdr_adjusted A boolean. If TRUE, use the the FDR- adjusted p-values. 
+#' Otherwise, get significant pairs based on the theoretical F-statistic cutoff.
+#' @return An adjacency matrix.
+#'
+#' @export
+getAdjacencyFstat <- 
+  function(object, pval = 0.05, fdr_adjusted = TRUE) {
+  
+  # get matrix
+  mat <- getMatrix(object)
+
+  # create empty adjacency matrix
+  adj <- matrix(0, nrow = nrow(mat), ncol = ncol(mat))
+  rownames(adj) <- rownames(mat)
+  colnames(adj) <- colnames(mat)
+
+  # fill in significant pairs
+  cutoff <- getCutoffFstat(object, pval, fdr_adjusted = fdr_adjusted)
+  if (cutoff) adj[mat <= cutoff] <- 1
+
+  return(adj)  
+}
diff --git a/R/3-shared-getCutoff.R b/R/3-shared-getCutoff.R
new file mode 100644
index 0000000..c8cd691
--- /dev/null
+++ b/R/3-shared-getCutoff.R
@@ -0,0 +1,192 @@
+#' Get a meaningful cutoff based on the FDR values from permutation tests.
+#' 
+#' @param object A \code{propd} or \code{propr} object.
+#' @param fdr A float value for the false discovery rate.
+#' Default is 0.05.
+#' @param window_size An integer. Default is 1. When it is greater than 1,
+#' the function will return a meaningful cutoff based on the moving
+#' average of the FDR values. This is useful when the FDR values are
+#' noisy and the user wants to smooth them out.
+#' @param positive A boolean. In the case of \code{propr} object, 
+#' toggles whether to return the cutoff for positive values (>=0). 
+#' Otherwise, return cutoff for negative values (<0). Default is TRUE.
+#' This argument is ignored for \code{propd} objects.
+#' @return A cutoff value.
+#' @export
+getCutoffFDR <- function(object, fdr = 0.05, window_size = 1, positive = TRUE) {
+  if (!"fdr" %in% slotNames(object)) {
+    stop("Please run updateCutoffs() before calling this function.")
+  }
+  if (nrow(object@fdr) == 0) {
+    stop("No FDR values found. Please run updateCutoffs() before calling this function.")
+  }
+  if (fdr < 0 | fdr > 1) {
+    stop("Provide a FDR cutoff from [0, 1].")
+  }
+
+  if (inherits(object, "propr")) {
+    return(getCutoffFDR.propr(object, fdr=fdr, window_size=window_size, positive=positive))
+
+  } else if (inherits(object, "propd")) {
+    return(getCutoffFDR.propd(object, fdr=fdr, window_size=window_size))
+    
+  } else{
+    stop("Provided 'object' not recognized.")
+  }
+}
+
+#' @rdname getCutoffFDR
+#' @section Methods:
+#' \code{getCutoffFDR.propr:}
+#' Same as \code{getCutoffFDR}, but for \code{propr} objects.
+#' @export
+getCutoffFDR.propr <- function(object, fdr = 0.05, window_size = 1, positive = TRUE) {
+  if (!positive & !object@direct) {
+    stop("Negative cutoffs are not allowed for this metric.")
+  }
+  if (!positive & !any(object@fdr$cutoff < 0, na.rm=T)){
+    stop("There are not negative values. Please check.")
+  }
+
+  # subset FDR data frame, to only include positive/negative tail
+  df <- object@fdr
+  if (positive) { df <- df[df$cutoff >= 0,] } else { df <- df[df$cutoff < 0,] }
+
+  # apply moving average to FDR values
+  if (window_size > 1) {
+    message("Applying moving average to FDR values.")
+    df$FDR <- getMovingAverage(df$FDR, window_size)
+  }
+
+  # get index of FDR values below the threshold
+  index <- (df$FDR <= fdr) & (is.finite(df$FDR))
+  if (!any(index)) {
+    warning("No significant cutoff found for the given FDR, when positive = ", positive)
+    return(FALSE)
+  }
+
+  # get cutoff
+  if (object@direct) {
+    cutoff <- min(abs(df$cutoff[index]))
+    if (!positive) cutoff <- -cutoff
+  } else{
+    cutoff <- max(df$cutoff[index])
+  }
+
+  return(cutoff)
+}
+
+#' @rdname getCutoffFDR
+#' @section Methods:
+#' \code{getCutoffFDR.propd:}
+#' Same as \code{getCutoffFDR}, but for \code{propd} objects.
+#' @export
+getCutoffFDR.propd <- function(object, fdr = 0.05, window_size = 1) {
+
+  # get df
+  df <- object@fdr
+  if (any(df$cutoff < 0)) stop("Negative cutoffs were found for theta values. Please check")
+
+  # apply moving average to FDR values
+  if (window_size > 1) {
+    message("Applying moving average to FDR values.")
+    df$FDR <- getMovingAverage(df$FDR, window_size)
+  }
+  
+  # get index of FDR values below the threshold
+  index <- (df$FDR <= fdr) & (is.finite(df$FDR))
+  if (!any(index)) {
+    warning("No significant cutoff found for the given FDR.")
+    return(FALSE)
+  }
+
+  # get cutoff
+  cutoff <- max(df$cutoff[index])
+  return(cutoff)
+}
+
+#' Calculate a theta Cutoff based on the F-statistic.
+#'
+#' This function uses the F distribution to calculate a cutoff of
+#'  theta for a p-value given by the \code{pval} argument.
+#'
+#' If the argument \code{fdr = TRUE}, this function returns the
+#'  empiric cutoff that corresponds to the FDR-adjusted p-value
+#'  stored in the \code{@@results$FDR} slot.
+#'
+#' @param object A \code{\link{propd}} object.
+#' @param pval A p-value at which to calculate a theta cutoff.
+#' @param fdr_adjusted A boolean. Toggles whether to calculate the theta
+#' cutoff for an FDR-adjusted p-value.
+#' @return A cutoff of theta from [0, 1].
+#' @export
+getCutoffFstat <- function(object, pval = 0.05, fdr_adjusted = FALSE) {
+  if (!"Fstat" %in% colnames(object@results)) {
+    stop("Please run updateF() on propd object before.")
+  }
+  if (pval < 0 | pval > 1) {
+    stop("Provide a p-value cutoff from [0, 1].")
+  }
+
+  if (fdr_adjusted) {
+    message("Alert: Returning an empiric cutoff based on the $FDR slot.")
+    index <- (object@results$FDR < pval) & (is.finite(object@results$FDR))
+    if (any(index)) {
+      cutoff <- max(object@results$theta[index])
+    } else{
+      warning("No significant cutoff found for the given p-value.")
+      cutoff <- FALSE
+    }
+
+  } else{
+    message("Alert: Returning an cutoff based on the F-statistic.")
+    # Compute based on theory
+    K <- length(unique(object@group))
+    N <- length(object@group) + object@dfz # population-level metric (i.e., N)
+    Q <- stats::qf(pval, K - 1, N - K, lower.tail = FALSE)
+    # # Fstat <- (N - 2) * (1 - object@theta$theta) / object@theta$theta
+    # # Q = Fstat
+    # # Q = (N-2) * (1-theta) / theta
+    # # Q / (N-2) = (1/theta) - 1
+    # # 1/theta = Q / (N-2) + 1 = Q(N-2)/(N-2)
+    # # theta = (N-2)/(Q+(N-2))
+    cutoff <- (N - 2) / (Q + (N - 2))
+  }
+
+  return(cutoff)
+}
+
+#' Caclulate the moving average of a vector.
+#' @param values A numeric vector.
+#' @param window_size An integer. The size of the window to calculate the
+#' moving average. Default is 1.
+getMovingAverage <- function(values, window_size = 1) {
+
+  if (any(is.na(values))) {
+    message("Moving averages are calculated for a vector containing NAs.")
+  }
+
+  # Initialize the result vector
+  n <- length(values)
+  result <- numeric(n)
+  
+  for (i in 1:n) {
+    # Determine the window indices
+    if (window_size %% 2 == 0) {
+      start_idx <- max(1, i - (window_size / 2 - 1))
+      end_idx <- min(n, i + window_size / 2)
+    }else{
+      start_idx <- max(1, i - floor(window_size / 2))
+      end_idx <- min(n, i + floor(window_size / 2))
+    }
+    
+    # Calculate the average for the current window
+    if (is.finite(values[i])){
+      result[i] <- mean(values[start_idx:end_idx], na.rm=TRUE)  # NA values are removed, to avoid propagation of NAs
+    }else{
+      result[i] <- values[i]  # this keeps the NA values corresponding to that position
+    }
+  }
+  
+  return(result)
+}
\ No newline at end of file
diff --git a/R/3-shared-getMatrix.R b/R/3-shared-getMatrix.R
index bb30a5b..8739d0e 100644
--- a/R/3-shared-getMatrix.R
+++ b/R/3-shared-getMatrix.R
@@ -8,8 +8,20 @@
 #' @return A matrix.
 #'
 #' @export
-getMatrix <-
-  function(object) {
-    df <- object@matrix
-    return(df)
+getMatrix <- function(object) {
+
+  if(class(object) == "propr"){
+    mat <- object@matrix
+
+  }else if(class(object) == "propd"){
+    mat <- vector2mat(object@results$theta, object@results$Partner, object@results$Pair)
+    rownames(mat) <- colnames(object@counts)
+    colnames(mat) <- colnames(object@counts)
+    diag(mat) <- 0
+
+  }else{
+    stop("Provided 'object' not recognized.")
+  }
+
+  return(mat)
   }
diff --git a/R/3-shared-getRatios.R b/R/3-shared-getRatios.R
index 07b121d..36ac592 100644
--- a/R/3-shared-getRatios.R
+++ b/R/3-shared-getRatios.R
@@ -56,11 +56,9 @@ getRatios <- function(object, switch = TRUE) {
 ratios <- function(matrix, alpha = NA) {
   lab <- labRcpp(ncol(matrix))
 
-  # Replace count zeros with 1 if appropriate
+  # Replace count zeros if appropriate
   if (any(as.matrix(matrix) == 0) & is.na(alpha)) {
-    message("Alert: Replacing 0s with next smallest value.")
-    zeros <- matrix == 0
-    matrix[zeros] <- min(matrix[!zeros])
+    matrix <- simple_zero_replacement(matrix)
   }
 
   # Get (log-)ratios [based on alpha]
diff --git a/R/3-shared-getResults.R b/R/3-shared-getResults.R
index 3e3c5b1..24cf451 100644
--- a/R/3-shared-getResults.R
+++ b/R/3-shared-getResults.R
@@ -10,9 +10,151 @@
 #' @export
 getResults <-
   function(object) {
-    df <- object@results
+    results <- object@results
     names <- colnames(object@counts)
-    df$Partner <- names[df$Partner]
-    df$Pair <- names[df$Pair]
-    return(df)
+    results$Partner <- names[results$Partner]
+    results$Pair <- names[results$Pair]
+    return(results)
   }
+
+#' Get Significant Results from Object based on the permutation tests.
+#'
+#' This function provides a unified wrapper to retrieve results
+#'  from a \code{propr} or \code{propd} object keeping only the
+#'  statistically significant pairs.
+#'
+#' @inheritParams getAdjacencyFDR
+#' @return A \code{data.frame} of results.
+#'
+#' @export
+getSignificantResultsFDR <-
+  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+
+    if (inherits(object, "propr")) {
+      results <- getSignificantResultsFDR.propr(object, fdr=fdr, window_size=window_size, tails=tails)
+
+    } else if(inherits(object, "propd")) {
+      results <- getSignificantResultsFDR.propd(object, fdr=fdr, window_size=window_size)
+
+    } else {
+      stop("Please provide a 'propr' or 'propd' object.")
+    }
+
+    return(results)
+}
+
+#' @rdname getSignificantResultsFDR
+#' @section Methods:
+#' \code{getSignificantResultsFDR.propd:}
+#' This function retrieves results from a \code{propd} object keeping 
+#' only the statistically significant pairs.
+#' @export
+getSignificantResultsFDR.propd <- 
+  function(object, fdr = 0.05, window_size = 1) {
+
+    results <- getResults(object)
+    cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size)
+    if (cutoff) {
+      return(results[which(results$theta <= cutoff), ])
+    } else {
+      return(results[0,])
+    }
+}
+
+#' @rdname getSignificantResultsFDR
+#' @section Methods:
+#' \code{getSignificantResultsFDR.propr:}
+#' This function retrieves results from a \code{propr} object keeping 
+#' only the statistically significant pairs.
+#' @export
+getSignificantResultsFDR.propr <- 
+  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+
+    if (is.null(tails)) {
+      tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
+    }
+    if (tails != 1 & tails != 2) {
+        stop("Please provide a valid value for tails: 1 or 2.")
+    }
+    if (tails == 1 & object@has_meaningful_negative_values) {
+      warning("Significant pairs are chosen based on one-sided FDR test.")
+    }
+    if (tails == 2 & !object@direct) {
+      stop("Two-sided FDR is not available for this metric.")
+    }
+
+    # function to subset the results data frame based on the cutoff
+    subsetBeyondCutoff <- function(data, cutoff) {
+      if (!cutoff) return(data[0,])  # return empty data frame when no cutoff found
+      if (object@direct) {
+        return(data[which(abs(data$propr) >= abs(cutoff)), ])
+      } else {
+        return(data[which(data$propr <= cutoff), ])
+      }
+    }
+
+    # correct fdr when two-sided
+    if (tails == 2) fdr <- fdr / 2
+
+    # define results data frame
+    df <- getResults(object)
+    results <- data.frame(matrix(ncol = ncol(df), nrow = 0))
+    colnames(results) <- colnames(df)
+
+    # get the significant positive values
+    cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size)
+    part <- df[which(df$propr >= 0),]
+    results <- rbind(results, subsetBeyondCutoff(part, cutoff))
+
+    # get the significant negative values
+    if (tails == 2) {
+      cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size, positive=F)
+      part <- df[which(df$propr < 0),]
+      results <- rbind(results, subsetBeyondCutoff(part, cutoff))
+    }
+
+    return(results)
+}
+
+#' Get Significant Results based on the F-stats.
+#'
+#' This function provides a unified wrapper to retrieve results
+#'  from a \code{propd} object keeping only the statistically 
+#'  significant pairs. Note that it can only be applied to theta_d,
+#'  as updateF only works for theta_d.
+#'
+#' @inheritParams getAdjacencyFstat
+#' @return A \code{data.frame} of results.
+#'
+#' @export
+getSignificantResultsFstat <- 
+  function(object, pval = 0.05, fdr_adjusted = TRUE) {
+
+    if (!"Fstat" %in% colnames(object@results)) {
+      stop("Please run updateF() on propd object before.")
+    }
+    if (pval < 0 | pval > 1) {
+      stop("Provide a p-value cutoff from [0, 1].")
+    }
+
+    # get results data frame
+    results <- getResults(object)
+
+    # get significant theta based on the FDR adjusted empirical p-values
+    if (fdr_adjusted) {
+      message("Alert: Returning the significant pairs based on the FDR adjusted p-values.")
+      results <- results[which(results$FDR <= pval), ]
+
+    # get siniginicant theta based on the F-statistic cutoff
+    } else {
+      message("Alert: Returning the significant pairs based on the F-statistic cutoff.")
+      cutoff <- getCutoffFstat(object, pval = pval, fdr_adjusted = FALSE)
+      if (cutoff) {
+        results <- results[which(results$theta <= cutoff), ]
+      } else {
+        results <- results[0,]
+      }
+    }
+
+    return(results)
+}
diff --git a/R/3-shared-graflex.R b/R/3-shared-graflex.R
new file mode 100644
index 0000000..01d7dc7
--- /dev/null
+++ b/R/3-shared-graflex.R
@@ -0,0 +1,51 @@
+#' Calculate odds ratio and FDR 
+#'
+#' This function calls \code{\link{graflex}} for each
+#' concept (i.e., column) in the database \code{K}.
+#'
+#' For each concept, this function calculates the odds ratio
+#' and determines the false discovery rate (FDR) by counting 
+#' the number of the actual OR was greater or less than a 
+#' permuted OR.
+#' 
+#' @param A An adjacency matrix.
+#' @param K A knowledge database where each row is a graph node
+#'  and each column is a concept.
+#' @param p An integer. The number of permutation.
+#' 
+#' @export
+runGraflex <- function(A, K, p=100, ncores=1) {
+  if (nrow(A) != nrow(K))
+    stop("'A' and 'K' must have identical rows.")
+  if (nrow(A) != ncol(A))
+    stop("'A' must be a square matrix.")
+
+  if (ncores == 1){
+
+    # for each knowledge network, calculate odds ratio and FDR
+    res <- lapply(1:ncol(K), function(k) {
+      Gk <- K[, k] %*% t(K[, k])      # converts the k column into an adjacency matrix (genes x genes)
+      graflex(A, Gk, p=p)  # this calls the graflex function implemented in Rcpp C++
+    })
+
+  }else{
+    packageCheck("parallel")
+    cl <- parallel::makeCluster(ncores)
+    res <- parallel::parLapply(cl, 1:ncol(K), function(k) {
+      Gk <- K[, k] %*% t(K[, k])
+      graflex(A, Gk, p=p)
+    })
+    parallel::stopCluster(cl)
+  }
+
+  # parse resulting data frame
+  res <- do.call("rbind", res)
+  res <- cbind(res, rep(p, ncol(K)))
+  res <- cbind(res, colnames(K))
+  res <- as.data.frame(res)
+  colnames(res) <- c("Neither", "G.only", "A.only", "Both", "Odds", "LogOR", "FDR.under", "FDR.over", "Permutes", "Concept")
+  # change the values to numeric, except for the concept column
+  res[,1:9] <- lapply(res[,1:9], as.numeric)
+
+  return(res)
+}
diff --git a/R/3-shared-updateCutoffs.R b/R/3-shared-updateCutoffs.R
index 23ee0a4..32f94a7 100644
--- a/R/3-shared-updateCutoffs.R
+++ b/R/3-shared-updateCutoffs.R
@@ -6,51 +6,111 @@
 #'  \code{updateCutoffs.propd}.
 #'
 #' @param object A \code{propr} or \code{propd} object.
-#' @param cutoff For \code{updateCutoffs}, a numeric vector.
-#'  this argument provides the FDR cutoffs to test.
+#' @param number_of_cutoffs An integer. The number of cutoffs to test.
+#' Given this number, the cutoffs will be determined based on the quantile of the data.
+#' In this way, the cutoffs will be evenly spaced across the data.
+#' @param custom_cutoffs A numeric vector. When provided, this vector is used
+#' as the FDR cutoffs to test, and number_of_cutoffs is ignored.
+#' @param tails 1 for one-sided on the right, and 2 for two-sided. When two-sided,
+#' the FDR is calculated for both positive and negative values. When NULL, use default
+#' option for the given metric. This is only relevant for \code{propr} objects.
 #' @param ncores An integer. The number of parallel cores to use.
 #' @return A \code{propr} or \code{propd} object with the FDR slot updated.
 #' @export
 updateCutoffs <-
   function(object,
-           cutoff = seq(.05, .95, .3),
+           number_of_cutoffs = 100,
+           custom_cutoffs = NULL,
+           tails = NULL,
            ncores = 1) {
-    if (inherits(object, "propr")) {
-      if (ncores == 1) {
-        message("Alert: Try parallelizing updateCutoffs with ncores > 1.")
-      }
 
-      updateCutoffs.propr(object, cutoff, ncores)
+    if (!is.null(tails)){
+      if (tails != 1 & tails != 2) stop("Tails only accept values 1 or 2, if not NULL.")
+    }
 
-    } else if (inherits(object, "propd")) {
-      if (ncores > 1) {
-        message("Alert: Parallel updateCutoffs not yet supported.")
-      }
+    get_cutoffs <- function(values, number_of_cutoffs, custom_cutoffs, tails=1) {
+      if (!is.null(custom_cutoffs)) return(custom_cutoffs)
+      if (tails == 1) values <- values[values >= 0]
+      return(as.numeric(quantile(values, probs = seq(0, 1, length.out = number_of_cutoffs))))
+    }
 
-      updateCutoffs.propd(object, cutoff)
+    if (inherits(object, "propr")) {
+      if (is.null(tails)) tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
+      if (tails == 1 & object@has_meaningful_negative_values) warning("One-sided FDR test is performed.")
+      if (tails == 2 & !object@direct) stop("Two-sided FDR is not available for this metric.")
+      values <- object@results$propr
+      cutoffs <- get_cutoffs(values, number_of_cutoffs, custom_cutoffs, tails)
+      updateCutoffs.propr(object, cutoffs, ncores)
+
+    } else if (inherits(object, "propd")) {
+      values <- object@results$theta
+      cutoffs <- get_cutoffs(values, number_of_cutoffs, custom_cutoffs, tails=1)
+      updateCutoffs.propd(object, cutoffs, ncores)
 
     } else{
       stop("Provided 'object' not recognized.")
     }
+
   }
 
 #' @rdname updateCutoffs
 #' @section Methods:
 #' \code{updateCutoffs.propr:}
-#'  Use the \code{propr} object to permute proportionality
+#'  Use the \code{propr} object to permute correlation-like metrics
+#'  (ie. rho, phi, phs, cor, pcor, pcor.shrink, pcor.bshrink),
 #'  across a number of cutoffs. Since the permutations get saved
 #'  when the object is created, calling \code{updateCutoffs}
 #'  will use the same random seed each time.
 #' @export
 updateCutoffs.propr <-
-  function(object, cutoff, ncores) {
+  function(object, cutoffs, ncores) {
 
-    metric_is_up <- function(metric) {
-      metrics <- c("rho", "cor", "pcor", "pcor.shrink", "pcor.bshrink")
-      return(metric %in% metrics)
+    if (identical(object@permutes, list(NULL))) {
+      stop("Permutation testing is disabled.")
     }
+    if (object@metric == "phi") {
+      warning("We recommend using the symmetric phi 'phs' for FDR permutation.")
+    }
+
+    # Set up FDR cutoff table
+    FDR <- as.data.frame(matrix(0, nrow = length(cutoffs), ncol = 4))
+    colnames(FDR) <- c("cutoff", "randcounts", "truecounts", "FDR")
+    FDR$cutoff <- cutoffs
+
+    # count the permuted values greater or less than each cutoff
+    if (ncores > 1) {
+      FDR$randcounts <- getFdrRandcounts.propr.parallel(object, cutoffs, ncores)
+    } else{
+      FDR$randcounts <- getFdrRandcounts.propr.run(object, cutoffs)
+    }
+
+    # count actual values greater or less than each cutoff
+    FDR$truecounts <- sapply(FDR$cutoff, function(cutoff) {
+      countValuesBeyondThreshold(object@results$propr, cutoff, object@direct)
+    })
+
+    # calculate FDR
+    FDR$FDR <- FDR$randcounts / FDR$truecounts
+
+    # Initialize @fdr
+    object@fdr <- FDR
+
+    return(object)
+  }
 
+#' This function counts the permuted values greater or less than each cutoff,
+#' using parallel processing, for a propr object
+getFdrRandcounts.propr.parallel <- 
+  function(object, cutoffs, ncores) {
+
+    # Set up the cluster
+    packageCheck("parallel")
+    cl <- parallel::makeCluster(ncores)
+    # parallel::clusterEvalQ(cl, requireNamespace(propr, quietly = TRUE))
+
+    # define function to parallelize
     getFdrRandcounts <- function(ct.k) {
+      # calculate permuted propr
       pr.k <- suppressMessages(propr::propr(
         ct.k,
         object@metric,
@@ -58,115 +118,60 @@ updateCutoffs.propr <-
         alpha = object@alpha,
         p = 0
       ))
-
       # Vector of propr scores for each pair of taxa.
       pkt <- pr.k@results$propr
-
-      # Find number of permuted theta less than cutoff
-      sapply(FDR$cutoff, function(cut) {
-        if (metric_is_up(object@metric)) {
-          count_greater_than(pkt, cut)
-        } else{
-          count_less_than(pkt, cut)
-        }
-      })
+      # Find number of permuted theta more or less than cutoff
+      sapply(cutoffs, function(cutoff) countValuesBeyondThreshold(pkt, cutoff, object@direct))
     }
 
-    if (object@metric == "rho") {
-      message("Alert: Estimating FDR for largely positive proportional pairs only.")
-    }
+    # Each element of this list will be a vector whose elements
+    # are the count of theta values less than the cutoff.
+    randcounts <- parallel::parLapply(cl = cl,
+                                      X = object@permutes,
+                                      fun = getFdrRandcounts)
 
-    if (object@metric == "phi") {
-      warning("We recommend using the symmetric phi 'phs' for FDR permutation.")
-    }
+    # get the average randcounts across all permutations
+    randcounts <- apply(as.data.frame(randcounts), 1, sum)
+    randcounts <- randcounts / length(object@permutes)
 
-    if (identical(object@permutes, list(NULL)))
-      stop("Permutation testing is disabled.")
-
-    # Let NA cutoff skip function
-    if (identical(cutoff, NA))
-      return(object)
+    # Explicitly stop the cluster.
+    parallel::stopCluster(cl)
 
-    # Set up FDR cutoff table
-    FDR <- as.data.frame(matrix(0, nrow = length(cutoff), ncol = 4))
-    colnames(FDR) <- c("cutoff", "randcounts", "truecounts", "FDR")
-    FDR$cutoff <- cutoff
-    p <- length(object@permutes)
-
-    if (ncores > 1) {
-      packageCheck("parallel")
-
-      # Set up the cluster and require propr
-      cl <- parallel::makeCluster(ncores)
-      # parallel::clusterEvalQ(cl, requireNamespace(propr, quietly = TRUE))
-
-      # Each element of this list will be a vector whose elements
-      # are the count of theta values less than the cutoff.
-      randcounts <- parallel::parLapply(cl = cl,
-                                        X = object@permutes,
-                                        fun = getFdrRandcounts)
+    return(randcounts)
+  }
 
-      # Sum across cutoff values
-      FDR$randcounts <- apply(as.data.frame(randcounts), 1, sum)
+#' This function counts the permuted values greater or less than each cutoff,
+#' using a single core, for a propr object.
+getFdrRandcounts.propr.run <-
+  function(object, cutoffs) {
 
-      # Explicitly stop the cluster.
-      parallel::stopCluster(cl)
+    # create empty randcounts
+    randcounts <- rep(0, length(cutoffs))
 
-    } else{
-      # Calculate propr for each permutation -- NOTE: `select` and `subset` disable permutation testing
-      for (k in 1:p) {
-        numTicks <- progress(k, p, numTicks)
-
-        # Calculate propr exactly based on @metric, @ivar, and @alpha
-        ct.k <- object@permutes[[k]]
-        pr.k <- suppressMessages(propr(
-          ct.k,
-          object@metric,
-          ivar = object@ivar,
-          alpha = object@alpha,
-          p = 0
-        ))
-        pkt <- pr.k@results$propr
-
-        # Find number of permuted theta less than cutoff
-        for (cut in 1:nrow(FDR)) {
-          # randcounts as cumsum
-
-          # Count positives as rho > cutoff, cor > cutoff, phi < cutoff, phs < cutoff
-          if (metric_is_up(object@metric)) {
-            FDR[cut, "randcounts"] <-
-              FDR[cut, "randcounts"] + count_greater_than(pkt, FDR[cut, "cutoff"])
-          } else{
-            # phi & phs
-            FDR[cut, "randcounts"] <-
-              FDR[cut, "randcounts"] + count_less_than(pkt, FDR[cut, "cutoff"])
-          }
-        }
-      }
-    }
+    # Calculate propr for each permutation -- NOTE: `select` and `subset` disable permutation testing
+    p <- length(object@permutes)
+    for (k in 1:p) {
+      numTicks <- progress(k, p, numTicks)
 
-    # Calculate FDR based on real and permuted tallys
-    FDR$randcounts <- FDR$randcounts / p # randcounts as mean
+      # Calculate propr exactly based on @metric, @ivar, and @alpha
+      ct.k <- object@permutes[[k]]
+      pr.k <- suppressMessages(propr(
+        ct.k,
+        object@metric,
+        ivar = object@ivar,
+        alpha = object@alpha,
+        p = 0
+      ))
+      pkt <- pr.k@results$propr
 
-    for (cut in 1:nrow(FDR)) {
-      # Count positives as rho > cutoff, cor > cutoff, phi < cutoff, phs < cutoff
-      if (metric_is_up(object@metric)) {
-        FDR[cut, "truecounts"] <-
-          sum(object@results$propr > FDR[cut, "cutoff"])
-      } else{
-        # phi & phs
-        FDR[cut, "truecounts"] <-
-          sum(object@results$propr < FDR[cut, "cutoff"])
+      # calculate the cumulative (across permutations) number of permuted values more or less than cutoff
+      for (cut in 1:length(cutoffs)){
+        randcounts[cut] <- randcounts[cut] + countValuesBeyondThreshold(pkt, cutoffs[cut], object@direct)
       }
-
-      FDR[cut, "FDR"] <-
-        FDR[cut, "randcounts"] / FDR[cut, "truecounts"]
     }
 
-    # Initialize @fdr
-    object@fdr <- FDR
-
-    return(object)
+    randcounts <- randcounts / p  # averaged across permutations
+    return(randcounts)
   }
 
 #' @rdname updateCutoffs
@@ -178,78 +183,186 @@ updateCutoffs.propr <-
 #'  will use the same random seed each time.
 #' @export
 updateCutoffs.propd <-
-  function(object, cutoff = seq(.05, .95, .3)) {
+  function(object, cutoffs, ncores) {
     if (identical(object@permutes, data.frame()))
       stop("Permutation testing is disabled.")
 
-    # Let NA cutoff skip function
-    if (identical(cutoff, NA))
-      return(object)
-
     # Set up FDR cutoff table
-    FDR <- as.data.frame(matrix(0, nrow = length(cutoff), ncol = 4))
+    FDR <- as.data.frame(matrix(0, nrow = length(cutoffs), ncol = 4))
     colnames(FDR) <- c("cutoff", "randcounts", "truecounts", "FDR")
-    FDR$cutoff <- cutoff
-    p <- ncol(object@permutes)
-    lrv <- object@results$lrv
+    FDR$cutoff <- cutoffs
 
-    # Use calculateTheta to permute active theta
-    for (k in 1:p) {
-      numTicks <- progress(k, p, numTicks)
+    # Count the permuted values greater or less than each cutoff
+    if (ncores > 1) {
+      FDR$randcounts <- getFdrRandcounts.propd.parallel(object, cutoffs, ncores)
+    } else{
+      FDR$randcounts <- getFdrRandcounts.propd.run(object, cutoffs)
+    }
 
-      # Tally k-th thetas that fall below each cutoff
-      shuffle <- object@permutes[, k]
+    # count actual values greater or less than each cutoff
+    FDR$truecounts <- sapply(1:nrow(FDR), function(cut) {
+        countValuesBeyondThreshold(object@results$theta, FDR[cut, "cutoff"], direct=FALSE)
+    })
 
-      if (object@active == "theta_mod") {
-        # Calculate theta_mod with updateF (using i-th permuted object)
-        if (is.na(object@Fivar))
-          stop("Please re-run 'updateF' with 'moderation = TRUE'.")
-        propdi <- suppressMessages(
-          propd(
-            object@counts[shuffle,],
-            group = object@group,
-            alpha = object@alpha,
-            p = 0,
-            weighted = object@weighted
-          )
+    # Calculate FDR
+    FDR$FDR <- FDR$randcounts / FDR$truecounts
+
+    # Initialize @fdr
+    object@fdr <- FDR
+
+    return(object)
+  }
+
+#' This function counts the permuted values greater or less than each cutoff,
+#' using parallel processing, for a propd object.
+getFdrRandcounts.propd.parallel <- 
+  function(object, cutoffs, ncores) {
+
+    # Set up the cluster
+    packageCheck("parallel")
+    cl <- parallel::makeCluster(ncores)
+    # parallel::clusterEvalQ(cl, requireNamespace(propr, quietly = TRUE))
+
+    # define functions to parallelize
+    getFdrRandcountsMod <- function(k) {
+      if (is.na(object@Fivar)) stop("Please re-run 'updateF' with 'moderation = TRUE'.")
+      shuffle <- object@permutes[, k]
+      propdi <- suppressMessages(
+        propd(
+          object@counts[shuffle,],
+          group = object@group,
+          alpha = object@alpha,
+          p = 0,
+          weighted = object@weighted
         )
-        propdi <-
-          suppressMessages(updateF(propdi, moderated = TRUE, ivar = object@Fivar))
-        pkt <- propdi@results$theta_mod
+      )
+      propdi <- suppressMessages(updateF(propdi, moderated = TRUE, ivar = Fivar))
+      pkt <- propdi@results$theta_mod
+      sapply(cutoffs, function(cutoff) countValuesBeyondThreshold(pkt, cutoff, direct=FALSE))
+    }
+    getFdrRandcounts <- function(k) {
+      shuffle <- object@permutes[, k]
+      pkt <- suppressMessages(
+        calculate_theta(
+          object@counts[shuffle,],
+          object@group,
+          object@alpha,
+          object@results$lrv,
+          only = object@active,
+          weighted = object@weighted
+        )
+      )
+      sapply(cutoffs, function(cutoff) countValuesBeyondThreshold(pkt, cutoff, direct=FALSE))
+    }
+
+    # Each element of this list will be a vector whose elements
+    # are the count of theta values less than the cutoff.
+    func = ifelse(object@active == "theta_mod", getFdrRandcountsMod, getFdrRandcounts)
+    randcounts <- parallel::parLapply(cl = cl,
+                                      X = 1:ncol(object@permutes),
+                                      fun = func)
+
+    # get the average randcounts across all permutations
+    randcounts <- apply(as.data.frame(randcounts), 1, sum)
+    randcounts <- randcounts / length(object@permutes)
+
+    # Explicitly stop the cluster.
+    parallel::stopCluster(cl)
+
+    return(randcounts)
+  }
+
+#' This function counts the permuted values greater or less than each cutoff,
+#' using a single core, for a propd object.
+getFdrRandcounts.propd.run <- 
+  function(object, cutoffs) {
+
+    # create empty randcounts
+    randcounts <- rep(0, length(cutoffs))
+
+    # use calculateTheta to permute active theta
+    p <- ncol(object@permutes)
+    for (k in 1:p) {
+      numTicks <- progress(k, p, numTicks)
 
+      # calculate permuted theta values
+      if (object@active == "theta_mod") {
+        pkt <- suppressMessages(getPermutedThetaMod(object, k))
       } else{
-        # Calculate all other thetas directly (using calculateTheta)
-        pkt <- suppressMessages(
-          calculate_theta(
-            object@counts[shuffle,],
-            object@group,
-            object@alpha,
-            lrv,
-            only = object@active,
-            weighted = object@weighted
-          )
-        )
+        pkt <- suppressMessages(getPermutedTheta(object, k))
       }
 
-      # Find number of permuted theta less than cutoff
-      for (cut in 1:nrow(FDR)) {
-        # randcounts as cumsum
-        FDR[cut, "randcounts"] <-
-          FDR[cut, "randcounts"] + sum(pkt < FDR[cut, "cutoff"])
+      # calculate the cumulative (across permutations) number of permuted values more or less than cutoff
+      for (cut in 1:length(cutoffs)){
+        randcounts[cut] <- randcounts[cut] + countValuesBeyondThreshold(pkt, cutoffs[cut], direct=FALSE)
       }
     }
 
-    # Calculate FDR based on real and permuted tallys
-    FDR$randcounts <- FDR$randcounts / p # randcounts as mean
-    for (cut in 1:nrow(FDR)) {
-      FDR[cut, "truecounts"] <-
-        sum(object@results$theta < FDR[cut, "cutoff"])
-      FDR[cut, "FDR"] <-
-        FDR[cut, "randcounts"] / FDR[cut, "truecounts"]
-    }
+    randcounts <- randcounts / p  # averaged across permutations
+    return(randcounts)
+}
 
-    # Initialize @fdr
-    object@fdr <- FDR
+#' Get the theta mod values for a given permutation
+getPermutedThetaMod <- 
+  function(object, k) {
 
-    return(object)
+    if (is.na(object@Fivar)) stop("Please re-run 'updateF' with 'moderation = TRUE'.")
+
+    # Tally k-th thetas that fall below each cutoff
+    shuffle <- object@permutes[, k]
+
+    # Calculate theta_mod with updateF (using k-th permuted object)
+    propdi <- suppressMessages(
+      propd(
+        object@counts[shuffle,],
+        group = object@group,
+        alpha = object@alpha,
+        p = 0,
+        weighted = object@weighted
+      )
+    )
+    propdi <- suppressMessages(updateF(propdi, moderated = TRUE, ivar = Fivar))
+
+    return(propdi@results$theta_mod)
+}
+
+#' Get the theta values for a given permutation
+getPermutedTheta <- 
+  function(object, k) {
+
+    # Tally k-th thetas that fall below each cutoff
+    shuffle <- object@permutes[, k]
+
+    # Calculate all other thetas directly (using calculateTheta)
+    pkt <- suppressMessages(
+      calculate_theta(
+        object@counts[shuffle,],
+        object@group,
+        object@alpha,
+        object@results$lrv,
+        only = object@active,
+        weighted = object@weighted
+      )
+    )
+
+    return(pkt)
+  }
+
+#' Count Values Greater or Less Than a Threshold
+#'
+#' This function counts the number of values greater or less than a threshold.
+#' The direction depends on if a direct or inverse relationship is asked,
+#' as well as the sign of the threshold.
+#'
+#' @param values A numeric vector.
+#' @param cutoff A numeric value.
+#' @direct A logical value. If \code{TRUE}, direct relationship is considered.
+#' @return The number of values greater or less than the threshold.
+countValuesBeyondThreshold <- function(values, cutoff, direct){
+  if (direct){
+    func <- ifelse(cutoff >= 0, count_greater_than, count_less_than)
+  } else {
+    func <- count_less_than
   }
+  return(func(values, cutoff))
+}
diff --git a/R/3-shared-updatePermutes.R b/R/3-shared-updatePermutes.R
index 225eb4a..c8d919d 100644
--- a/R/3-shared-updatePermutes.R
+++ b/R/3-shared-updatePermutes.R
@@ -6,16 +6,16 @@
 #'  \code{updatePermutes.propd}.
 #'
 #' @param object A \code{propr} or \code{propd} object.
-#' @param p The number of permutations to perform.
+#' @param p The number of permutations to perform. Default is 100.
 #' @return A \code{propr} or \code{propd} object with the permutes slot updated.
 #' @export
 updatePermutes <-
-  function(object, p, fixseed = FALSE) {
+  function(object, p=100) {
     if (inherits(object, "propr")) {
-      updatePermutes.propr(object, p, fixseed)
+      updatePermutes.propr(object, p)
 
     } else if (inherits(object, "propd")) {
-      updatePermutes.propd(object, p, fixseed)
+      updatePermutes.propd(object, p)
 
     } else{
       stop("Provided 'object' not recognized.")
@@ -23,13 +23,12 @@ updatePermutes <-
   }
 
 updatePermutes.propr <-
-  function(object, p = 100, fixseed = FALSE) {
+  function(object, p) {
     # Shuffle row-wise so that per-sample CLR does not change
     message("Alert: Fixing permutations to active random seed.")
     ct <- object@counts
     permutes <- vector("list", p)
     for (ins in 1:p){
-      if (fixseed) set.seed(ins)
       permutes[[ins]] <- t(apply(ct, 1, sample))
     }
     object@permutes <- permutes
@@ -37,12 +36,11 @@ updatePermutes.propr <-
   }
 
 updatePermutes.propd <-
-  function(object, p = 100, fixseed = FALSE) {
+  function(object, p) {
     message("Alert: Fixing permutations to active random seed.")
     ct <- object@counts
     permutes <- as.data.frame(matrix(0, nrow = nrow(ct), ncol = p))
     for (col in 1:p){
-      if (fixseed) set.seed(col)
       permutes[, col] <- sample(1:nrow(ct))
     }
     object@permutes <- permutes
diff --git a/R/3-shared-zero.R b/R/3-shared-zero.R
new file mode 100644
index 0000000..a1fea64
--- /dev/null
+++ b/R/3-shared-zero.R
@@ -0,0 +1,23 @@
+#' Simple Zero Replacement in a Count Matrix
+#'
+#' This function replaces zeros with the next smallest non-zero value in the
+#'  input count matrix. If the matrix contains no zeros, it produces an
+#'  informational message indicating that no replacements were made.
+#'
+#' @param ct A data matrix containing numerical values.
+#' @return A matrix with zero values replaced by the next smallest non-zero value.
+#' If no zeros are found, the function returns the original matrix.
+#' @examples
+#' # Sample input count data with zeros
+#' data <- matrix(c(0, 2, 3, 4, 5, 0), nrow = 2, byrow = TRUE)
+#' @export
+simple_zero_replacement <- function(ct) {
+  if (any(ct == 0)) {
+    message("Alert: replacing zeros with minimun value.")
+    zeros <- ct == 0
+    ct[zeros] <- min(ct[!zeros])
+  } else{
+    message("Alert: No 0s found that need replacement.")
+  }
+  return(ct)
+}
diff --git a/R/5-selectReference.R b/R/5-selectReference.R
index 7edc46c..5b34da7 100644
--- a/R/5-selectReference.R
+++ b/R/5-selectReference.R
@@ -22,6 +22,8 @@
 #'
 #' @export
 selectReference <- function(counts, ivar, alpha) {
+  # replace zeros
+  counts <- simple_zero_replacement(counts)
   # Transform data into log space
   lr <- logratio(counts, ivar, alpha)
 
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 28a2cfd..60f02dc 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -73,10 +73,18 @@ half2mat <- function(X) {
     .Call(`_propr_half2mat`, X)
 }
 
+vector2mat <- function(X, i, j) {
+    .Call(`_propr_vector2mat`, X, i, j)
+}
+
 ratiosRcpp <- function(X) {
     .Call(`_propr_ratiosRcpp`, X)
 }
 
+results2matRcpp <- function(results, n, diagonal = 0.0) {
+    .Call(`_propr_results2matRcpp`, results, n, diagonal)
+}
+
 count_less_than <- function(x, cutoff) {
     .Call(`_propr_count_less_than`, x, cutoff)
 }
@@ -85,10 +93,50 @@ count_greater_than <- function(x, cutoff) {
     .Call(`_propr_count_greater_than`, x, cutoff)
 }
 
+count_less_equal_than <- function(x, cutoff) {
+    .Call(`_propr_count_less_equal_than`, x, cutoff)
+}
+
+count_greater_equal_than <- function(x, cutoff) {
+    .Call(`_propr_count_greater_equal_than`, x, cutoff)
+}
+
 ctzRcpp <- function(X) {
     .Call(`_propr_ctzRcpp`, X)
 }
 
+get_lower_triangle <- function(mat) {
+    .Call(`_propr_get_lower_triangle`, mat)
+}
+
+shuffle_and_get_lower_triangle <- function(mat) {
+    .Call(`_propr_shuffle_and_get_lower_triangle`, mat)
+}
+
+binTab <- function(A, G) {
+    .Call(`_propr_binTab`, A, G)
+}
+
+getOR <- function(A, G) {
+    .Call(`_propr_getOR`, A, G)
+}
+
+permuteOR <- function(A, Gstar, p = 100L) {
+    .Call(`_propr_permuteOR`, A, Gstar, p)
+}
+
+getFDR_over <- function(actual, permuted) {
+    .Call(`_propr_getFDR_over`, actual, permuted)
+}
+
+getFDR_under <- function(actual, permuted) {
+    .Call(`_propr_getFDR_under`, actual, permuted)
+}
+
+graflex <- function(A, G, p = 100L) {
+    .Call(`_propr_graflex`, A, G, p)
+}
+
 lr2vlr <- function(lr) {
     .Call(`_propr_lr2vlr`, lr)
 }
diff --git a/man/basis_shrinkage.Rd b/man/basis_shrinkage.Rd
deleted file mode 100644
index 5a45b22..0000000
--- a/man/basis_shrinkage.Rd
+++ /dev/null
@@ -1,35 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/1a-propr-backend.R
-\name{basis_shrinkage}
-\alias{basis_shrinkage}
-\title{Covariance Shrinkage and Partial Correlation Calculation}
-\usage{
-basis_shrinkage(M, outtype = c("clr", "alr"))
-}
-\arguments{
-\item{M}{A data matrix representing the counts. Each row should represent
-observations, and each column should represent different variables.}
-
-\item{outtype}{A character vector specifying the output type. It can take
-either "clr" (centered log-ratio) or "alr" (additive log-ratio).}
-}
-\value{
-A matrix representing the partial correlation matrix.
-}
-\description{
-This function performs covariance shrinkage and calculates the partial
- correlation matrix based on the input count data matrix. The function can
- output the results in two formats: centered log-ratio (clr) or
- additive log-ratio (alr).
-}
-\examples{
-# Sample input count data
-data <- iris[,1:4]
-
-# Calculate partial correlation matrix using clr transformation
-result_clr <- basis_shrinkage(data, outtype = "clr")
-
-# Calculate partial correlation matrix using alr transformation
-result_alr <- basis_shrinkage(data, outtype = "alr")
-
-}
diff --git a/man/binTab.Rd b/man/binTab.Rd
deleted file mode 100644
index 34396db..0000000
--- a/man/binTab.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
-\name{binTab}
-\alias{binTab}
-\title{Tabulate Overlap}
-\usage{
-binTab(A, G)
-}
-\arguments{
-\item{A, G}{A vector or adjacency matrix.}
-}
-\value{
-A table of overlap.
-}
-\description{
-This function tabulates the overlap between
- two vectors or two adjacency matrices.
- It is a faster version of \code{table} that
- only supports binary input.
-}
diff --git a/man/calculateOR.Rd b/man/calculateOR.Rd
deleted file mode 100644
index 46556f7..0000000
--- a/man/calculateOR.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
-\name{calculateOR}
-\alias{calculateOR}
-\title{Calculate Odds Ratio}
-\usage{
-calculateOR(A, G)
-}
-\description{
-This function calculates an odds ratio for the
- edge overlap between two non-random graphs.
-}
-\details{
-Note that this function calculates overlap for the
- lower-left triangle of the input matrices.
-}
diff --git a/man/countValuesBeyondThreshold.Rd b/man/countValuesBeyondThreshold.Rd
new file mode 100644
index 0000000..8cd9af0
--- /dev/null
+++ b/man/countValuesBeyondThreshold.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{countValuesBeyondThreshold}
+\alias{countValuesBeyondThreshold}
+\title{Count Values Greater or Less Than a Threshold}
+\usage{
+countValuesBeyondThreshold(values, cutoff, direct)
+}
+\arguments{
+\item{values}{A numeric vector.}
+
+\item{cutoff}{A numeric value.}
+}
+\value{
+The number of values greater or less than the threshold.
+}
+\description{
+This function counts the number of values greater or less than a threshold.
+The direction depends on if a direct or inverse relationship is asked,
+as well as the sign of the threshold.
+}
diff --git a/man/getAdjacencyFDR.Rd b/man/getAdjacencyFDR.Rd
new file mode 100644
index 0000000..8ecbf01
--- /dev/null
+++ b/man/getAdjacencyFDR.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getAdjacency.R
+\name{getAdjacencyFDR}
+\alias{getAdjacencyFDR}
+\alias{getAdjacencyFDR.propd}
+\alias{getAdjacencyFDR.propr}
+\title{Get Adjacency Matrix as indicated by permutation tests}
+\usage{
+getAdjacencyFDR(object, fdr = 0.05, window_size = 1, tails = NULL)
+
+getAdjacencyFDR.propd(object, fdr = 0.05, window_size = 1)
+
+getAdjacencyFDR.propr(object, fdr = 0.05, window_size = 1, tails = NULL)
+}
+\arguments{
+\item{object}{A \code{propd} or \code{propr} object.}
+
+\item{fdr}{A float value for the false discovery rate. Default is 0.05.}
+
+\item{window_size}{An integer. Default is 1. When it is greater than 1, the FDR
+values would be smoothed out by a moving average of the given window size.}
+
+\item{tails}{1 for one-sided on the right, and 2 for two-sided. When NULL, use default
+test for the given metric. This is only relevant for \code{propr} objects.}
+}
+\value{
+An adjacency matrix.
+}
+\description{
+This function gets the significant pairs, according to the permutation tests. 
+Then it fills the adjacency matrix with 1 if pair is significant, otherwise 0.
+}
+\section{Methods}{
+
+\code{getAdjacencyFDR.propd:}
+This function creates an adjacency matrix with the significant
+pairs from a \code{propd} object.
+
+
+\code{getAdjacencyFDR.propr:}
+This function creates an adjacency matrix with the significant
+pairs from a \code{propd} object.
+}
+
diff --git a/man/getAdjacencyFstat.Rd b/man/getAdjacencyFstat.Rd
new file mode 100644
index 0000000..74abca1
--- /dev/null
+++ b/man/getAdjacencyFstat.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getAdjacency.R
+\name{getAdjacencyFstat}
+\alias{getAdjacencyFstat}
+\title{Get Adjacency Matrix as indicated by F-statistics}
+\usage{
+getAdjacencyFstat(object, pval = 0.05, fdr_adjusted = TRUE)
+}
+\arguments{
+\item{object}{A \code{propd} or \code{propr} object.}
+
+\item{pval}{A float value for the p-value. Default is 0.05.}
+
+\item{fdr_adjusted}{A boolean. If TRUE, use the the FDR- adjusted p-values. 
+Otherwise, get significant pairs based on the theoretical F-statistic cutoff.}
+}
+\value{
+An adjacency matrix.
+}
+\description{
+This function gets the significant pairs, according to the F-statistics. 
+Then it fills the adjacency matrix with 1 if pair is significant, otherwise 0.
+Note that it can only be applied to theta_d, as updateF only works for theta_d.
+}
diff --git a/man/getCutoffFDR.Rd b/man/getCutoffFDR.Rd
new file mode 100644
index 0000000..9964164
--- /dev/null
+++ b/man/getCutoffFDR.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getCutoff.R
+\name{getCutoffFDR}
+\alias{getCutoffFDR}
+\alias{getCutoffFDR.propr}
+\alias{getCutoffFDR.propd}
+\title{Get a meaningful cutoff based on the FDR values from permutation tests.}
+\usage{
+getCutoffFDR(object, fdr = 0.05, window_size = 1, positive = TRUE)
+
+getCutoffFDR.propr(object, fdr = 0.05, window_size = 1, positive = TRUE)
+
+getCutoffFDR.propd(object, fdr = 0.05, window_size = 1)
+}
+\arguments{
+\item{object}{A \code{propd} or \code{propr} object.}
+
+\item{fdr}{A float value for the false discovery rate.
+Default is 0.05.}
+
+\item{window_size}{An integer. Default is 1. When it is greater than 1,
+the function will return a meaningful cutoff based on the moving
+average of the FDR values. This is useful when the FDR values are
+noisy and the user wants to smooth them out.}
+
+\item{positive}{A boolean. In the case of \code{propr} object, 
+toggles whether to return the cutoff for positive values (>=0). 
+Otherwise, return cutoff for negative values (<0). Default is TRUE.
+This argument is ignored for \code{propd} objects.}
+}
+\value{
+A cutoff value.
+}
+\description{
+Get a meaningful cutoff based on the FDR values from permutation tests.
+}
+\section{Methods}{
+
+\code{getCutoffFDR.propr:}
+Same as \code{getCutoffFDR}, but for \code{propr} objects.
+
+
+\code{getCutoffFDR.propd:}
+Same as \code{getCutoffFDR}, but for \code{propd} objects.
+}
+
diff --git a/man/runCutoff.Rd b/man/getCutoffFstat.Rd
similarity index 66%
rename from man/runCutoff.Rd
rename to man/getCutoffFstat.Rd
index f1bf980..e66c6f4 100644
--- a/man/runCutoff.Rd
+++ b/man/getCutoffFstat.Rd
@@ -1,17 +1,17 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2c-propd-experimental.R
-\name{runCutoff}
-\alias{runCutoff}
-\title{Calculate a theta Cutoff}
+% Please edit documentation in R/3-shared-getCutoff.R
+\name{getCutoffFstat}
+\alias{getCutoffFstat}
+\title{Calculate a theta Cutoff based on the F-statistic.}
 \usage{
-runCutoff(object, pval = 0.05, fdr = FALSE)
+getCutoffFstat(object, pval = 0.05, fdr_adjusted = FALSE)
 }
 \arguments{
 \item{object}{A \code{\link{propd}} object.}
 
 \item{pval}{A p-value at which to calculate a theta cutoff.}
 
-\item{fdr}{A boolean. Toggles whether to calculate the theta
+\item{fdr_adjusted}{A boolean. Toggles whether to calculate the theta
 cutoff for an FDR-adjusted p-value.}
 }
 \value{
diff --git a/man/getFDR.Rd b/man/getFDR.Rd
deleted file mode 100644
index 7841325..0000000
--- a/man/getFDR.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
-\name{getFDR}
-\alias{getFDR}
-\title{Calculate Odds Ratio FDR}
-\usage{
-getFDR(actual, permuted)
-}
-\arguments{
-\item{actual}{A result from \code{\link{calculateOR}}.}
-
-\item{permuted}{A result from \code{\link{permuteOR}}.}
-}
-\value{
-A \code{data.frame} of the FDRs for over-
- and under- enrichment.
-}
-\description{
-This function calculates the false discovery rate (FDR)
- for over- and under-enrichment by counting the number of
- times the actual OR was greater than
- (or less than) a permuted OR.
-}
diff --git a/man/getFdrRandcounts.propd.parallel.Rd b/man/getFdrRandcounts.propd.parallel.Rd
new file mode 100644
index 0000000..15faf55
--- /dev/null
+++ b/man/getFdrRandcounts.propd.parallel.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getFdrRandcounts.propd.parallel}
+\alias{getFdrRandcounts.propd.parallel}
+\title{This function counts the permuted values greater or less than each cutoff,
+using parallel processing, for a propd object.}
+\usage{
+getFdrRandcounts.propd.parallel(object, cutoffs, ncores)
+}
+\description{
+This function counts the permuted values greater or less than each cutoff,
+using parallel processing, for a propd object.
+}
diff --git a/man/getFdrRandcounts.propd.run.Rd b/man/getFdrRandcounts.propd.run.Rd
new file mode 100644
index 0000000..5bed408
--- /dev/null
+++ b/man/getFdrRandcounts.propd.run.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getFdrRandcounts.propd.run}
+\alias{getFdrRandcounts.propd.run}
+\title{This function counts the permuted values greater or less than each cutoff,
+using a single core, for a propd object.}
+\usage{
+getFdrRandcounts.propd.run(object, cutoffs)
+}
+\description{
+This function counts the permuted values greater or less than each cutoff,
+using a single core, for a propd object.
+}
diff --git a/man/getFdrRandcounts.propr.parallel.Rd b/man/getFdrRandcounts.propr.parallel.Rd
new file mode 100644
index 0000000..e9db844
--- /dev/null
+++ b/man/getFdrRandcounts.propr.parallel.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getFdrRandcounts.propr.parallel}
+\alias{getFdrRandcounts.propr.parallel}
+\title{This function counts the permuted values greater or less than each cutoff,
+using parallel processing, for a propr object}
+\usage{
+getFdrRandcounts.propr.parallel(object, cutoffs, ncores)
+}
+\description{
+This function counts the permuted values greater or less than each cutoff,
+using parallel processing, for a propr object
+}
diff --git a/man/getFdrRandcounts.propr.run.Rd b/man/getFdrRandcounts.propr.run.Rd
new file mode 100644
index 0000000..e11d9c8
--- /dev/null
+++ b/man/getFdrRandcounts.propr.run.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getFdrRandcounts.propr.run}
+\alias{getFdrRandcounts.propr.run}
+\title{This function counts the permuted values greater or less than each cutoff,
+using a single core, for a propr object.}
+\usage{
+getFdrRandcounts.propr.run(object, cutoffs)
+}
+\description{
+This function counts the permuted values greater or less than each cutoff,
+using a single core, for a propr object.
+}
diff --git a/man/getMovingAverage.Rd b/man/getMovingAverage.Rd
new file mode 100644
index 0000000..a7e47bd
--- /dev/null
+++ b/man/getMovingAverage.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getCutoff.R
+\name{getMovingAverage}
+\alias{getMovingAverage}
+\title{Caclulate the moving average of a vector.}
+\usage{
+getMovingAverage(values, window_size = 1)
+}
+\arguments{
+\item{values}{A numeric vector.}
+
+\item{window_size}{An integer. The size of the window to calculate the
+moving average. Default is 1.}
+}
+\description{
+Caclulate the moving average of a vector.
+}
diff --git a/man/getOR.Rd b/man/getOR.Rd
deleted file mode 100644
index 8fdef45..0000000
--- a/man/getOR.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
-\name{getOR}
-\alias{getOR}
-\title{Calculate Odds Ratio}
-\usage{
-getOR(A, G)
-}
-\value{
-A \code{data.frame} of results.
-}
-\description{
-This function calculates the overlap between
- two vectors or two adjacency matrices.
- It returns the OR as well as other metrics.
-}
diff --git a/man/getPermutedTheta.Rd b/man/getPermutedTheta.Rd
new file mode 100644
index 0000000..86a07c9
--- /dev/null
+++ b/man/getPermutedTheta.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getPermutedTheta}
+\alias{getPermutedTheta}
+\title{Get the theta values for a given permutation}
+\usage{
+getPermutedTheta(object, k)
+}
+\description{
+Get the theta values for a given permutation
+}
diff --git a/man/getPermutedThetaMod.Rd b/man/getPermutedThetaMod.Rd
new file mode 100644
index 0000000..f266d84
--- /dev/null
+++ b/man/getPermutedThetaMod.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-updateCutoffs.R
+\name{getPermutedThetaMod}
+\alias{getPermutedThetaMod}
+\title{Get the theta mod values for a given permutation}
+\usage{
+getPermutedThetaMod(object, k)
+}
+\description{
+Get the theta mod values for a given permutation
+}
diff --git a/man/getSignificantResultsFDR.Rd b/man/getSignificantResultsFDR.Rd
new file mode 100644
index 0000000..27c7183
--- /dev/null
+++ b/man/getSignificantResultsFDR.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getResults.R
+\name{getSignificantResultsFDR}
+\alias{getSignificantResultsFDR}
+\alias{getSignificantResultsFDR.propd}
+\alias{getSignificantResultsFDR.propr}
+\title{Get Significant Results from Object based on the permutation tests.}
+\usage{
+getSignificantResultsFDR(object, fdr = 0.05, window_size = 1, tails = NULL)
+
+getSignificantResultsFDR.propd(object, fdr = 0.05, window_size = 1)
+
+getSignificantResultsFDR.propr(
+  object,
+  fdr = 0.05,
+  window_size = 1,
+  tails = NULL
+)
+}
+\arguments{
+\item{object}{A \code{propd} or \code{propr} object.}
+
+\item{fdr}{A float value for the false discovery rate. Default is 0.05.}
+
+\item{window_size}{An integer. Default is 1. When it is greater than 1, the FDR
+values would be smoothed out by a moving average of the given window size.}
+
+\item{tails}{1 for one-sided on the right, and 2 for two-sided. When NULL, use default
+test for the given metric. This is only relevant for \code{propr} objects.}
+}
+\value{
+A \code{data.frame} of results.
+}
+\description{
+This function provides a unified wrapper to retrieve results
+ from a \code{propr} or \code{propd} object keeping only the
+ statistically significant pairs.
+}
+\section{Methods}{
+
+\code{getSignificantResultsFDR.propd:}
+This function retrieves results from a \code{propd} object keeping 
+only the statistically significant pairs.
+
+
+\code{getSignificantResultsFDR.propr:}
+This function retrieves results from a \code{propr} object keeping 
+only the statistically significant pairs.
+}
+
diff --git a/man/getSignificantResultsFstat.Rd b/man/getSignificantResultsFstat.Rd
new file mode 100644
index 0000000..9362a6c
--- /dev/null
+++ b/man/getSignificantResultsFstat.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/3-shared-getResults.R
+\name{getSignificantResultsFstat}
+\alias{getSignificantResultsFstat}
+\title{Get Significant Results based on the F-stats.}
+\usage{
+getSignificantResultsFstat(object, pval = 0.05, fdr_adjusted = TRUE)
+}
+\arguments{
+\item{object}{A \code{propd} or \code{propr} object.}
+
+\item{pval}{A float value for the p-value. Default is 0.05.}
+
+\item{fdr_adjusted}{A boolean. If TRUE, use the the FDR- adjusted p-values. 
+Otherwise, get significant pairs based on the theoretical F-statistic cutoff.}
+}
+\value{
+A \code{data.frame} of results.
+}
+\description{
+This function provides a unified wrapper to retrieve results
+ from a \code{propd} object keeping only the statistically 
+ significant pairs. Note that it can only be applied to theta_d,
+ as updateF only works for theta_d.
+}
diff --git a/man/index_reference.Rd b/man/index_reference.Rd
index 3ea6547..eff5fe6 100644
--- a/man/index_reference.Rd
+++ b/man/index_reference.Rd
@@ -15,8 +15,9 @@ It can take the following values:
  - "clr" (default): Centered log-ratio transformation.
  - "alr": Additive log-ratio transformation ("pcor.bshrink" metric only).
  - "iqlr": Inter-quartile log-ratio transformation from ALDEx2.
- - The explicit name(s) of variable(s) to use as a reference.
-Use NA to skip log-ratio transformation.}
+ - The explicit name(s) or index(es) of variable(s) to use as a reference.
+ - Use NA to skip log-ratio transformation and any other pre-processing, like
+zero replacement. This is useful when the input data is already pre-processed.}
 }
 \value{
 A numeric vector representing the indices of features selected.
diff --git a/man/logratio.Rd b/man/logratio.Rd
index c93baeb..3605c22 100644
--- a/man/logratio.Rd
+++ b/man/logratio.Rd
@@ -15,8 +15,9 @@ It can take the following values:
  - "clr" (default): Centered log-ratio transformation.
  - "alr": Additive log-ratio transformation ("pcor.bshrink" metric only).
  - "iqlr": Inter-quartile log-ratio transformation from ALDEx2.
- - The explicit name(s) of variable(s) to use as a reference.
-Use NA to skip log-ratio transformation.}
+ - The explicit name(s) or index(es) of variable(s) to use as a reference.
+ - Use NA to skip log-ratio transformation and any other pre-processing, like
+zero replacement. This is useful when the input data is already pre-processed.}
 
 \item{alpha}{The alpha parameter used in the alpha log-ratio transformation.}
 }
diff --git a/man/logratio_without_alpha.Rd b/man/logratio_without_alpha.Rd
index ca3a8ce..4b39bef 100644
--- a/man/logratio_without_alpha.Rd
+++ b/man/logratio_without_alpha.Rd
@@ -25,7 +25,7 @@ This function applies a log-ratio transformation to a given data matrix
 # Sample input data
 data <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, byrow = TRUE)
 
-# Applying log-ratio transformation to rows using columns 2 and 3
+# Applying log-ratio transformation to rows using columns 2 and 3 as reference
 result <- logratio_without_alpha(data, use = c(2, 3))
 
 }
diff --git a/man/pcor.bshrink.Rd b/man/pcor.bshrink.Rd
new file mode 100644
index 0000000..a325f83
--- /dev/null
+++ b/man/pcor.bshrink.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/1a-propr-backend.R
+\name{pcor.bshrink}
+\alias{pcor.bshrink}
+\title{Basis Covariance Shrinkage and Partial Correlation Calculation}
+\usage{
+pcor.bshrink(ct, outtype = c("clr", "alr"))
+}
+\arguments{
+\item{ct}{A data matrix representing the counts. Each row should represent
+observations, and each column should represent different variables.}
+
+\item{outtype}{A character vector specifying the output type. It can take
+either "clr" (centered log-ratio) or "alr" (additive log-ratio). Since
+the reference does not affect the logratio partial correlation coefficients,
+the index of the reference is not needed for "alr". Also, "clr" is recommended
+because of the same reason, at the same time avoiding losing one dimension.}
+}
+\value{
+A matrix representing the shrunk partial correlation matrix.
+}
+\description{
+This function performs covariance shrinkage on the basis matrix
+ and calculates the partial correlation matrix. The function can
+ output the results in two formats: centered log-ratio (clr) or
+ additive log-ratio (alr).
+}
+\examples{
+# Sample input count data
+data <- iris[,1:4]
+
+# Calculate partial correlation matrix using clr transformation
+result_clr <- pcor.bshrink(data, outtype = "clr")
+
+# Calculate partial correlation matrix using alr transformation
+result_alr <- pcor.bshrink(data, outtype = "alr")
+
+}
diff --git a/man/permuteOR.Rd b/man/permuteOR.Rd
deleted file mode 100644
index 9162559..0000000
--- a/man/permuteOR.Rd
+++ /dev/null
@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
-\name{permuteOR}
-\alias{permuteOR}
-\title{Permute Odds Ratio}
-\usage{
-permuteOR(A, G, p = 500, fixseed = FALSE)
-}
-\arguments{
-\item{A, G}{An adjacency matrix.}
-
-\item{p}{An integer. The number of overlaps to permute.}
-
-\item{fixseed}{A logical. If \code{TRUE}, the seed is fixed
-for reproducibility.}
-}
-\description{
-This function permutes \code{p} odds ratios for the
- edge overlap between two non-random graphs.
- It does this by randomly shuffling the rows and
- columns of \code{A} (jointly, thus preserving
- the degree distribution).
-}
-\details{
-Note that this function calculates overlap for the
- lower-left triangle of the input matrices.
-}
diff --git a/man/propd.Rd b/man/propd.Rd
index 9ff08d7..bc1b1ec 100755
--- a/man/propd.Rd
+++ b/man/propd.Rd
@@ -14,7 +14,7 @@
 \usage{
 \S4method{show}{propd}(object)
 
-propd(counts, group, alpha = NA, p = 0, fixseed = FALSE, weighted = FALSE)
+propd(counts, group, alpha = NA, p = 0, weighted = FALSE)
 
 setDisjointed(propd)
 
@@ -38,10 +38,6 @@ assignment of each count to different groups.}
 \item{p}{The number of permutations to perform for calculating the false
 discovery rate (FDR). The default is 0.}
 
-\item{fixseed}{A logical value indicating whether to fix the random seed 
-for generating the pemurations. If `TRUE`, the function will fix the random
-The default is `FALSE`.}
-
 \item{weighted}{A logical value indicating whether weighted calculations
 should be performed. If \code{TRUE}, the function will use limma-based
 weights for the calculations.}
diff --git a/man/propr.Rd b/man/propr.Rd
index dd181bd..ff2a400 100755
--- a/man/propr.Rd
+++ b/man/propr.Rd
@@ -11,13 +11,13 @@
 
 propr(
   counts,
-  metric = c("rho", "phi", "phs", "cor", "vlr", "pcor", "pcor.shrink", "pcor.bshrink"),
+  metric = c("rho", "phi", "phs", "cor", "vlr", "ppcor", "pcor", "pcor.shrink",
+    "pcor.bshrink"),
   ivar = "clr",
   select = NA,
   symmetrize = FALSE,
   alpha = NA,
   p = 0,
-  fixseed = FALSE,
   ...
 )
 }
@@ -46,8 +46,9 @@ It can take the following values:
  - "clr" (default): Centered log-ratio transformation.
  - "alr": Additive log-ratio transformation ("pcor.bshrink" metric only).
  - "iqlr": Inter-quartile log-ratio transformation from ALDEx2.
- - The explicit name(s) of variable(s) to use as a reference.
-Use NA to skip log-ratio transformation.}
+ - The explicit name(s) or index(es) of variable(s) to use as a reference.
+ - Use NA to skip log-ratio transformation and any other pre-processing, like
+zero replacement. This is useful when the input data is already pre-processed.}
 
 \item{select}{A numeric vector representing the indices of features to be
 used for computing the Propr matrix. This argument is optional. If
@@ -62,10 +63,6 @@ will symmetrize the matrix; otherwise, it will return the original matrix.}
 \item{p}{The number of permutations to perform for calculating the false
 discovery rate (FDR). The default is 0.}
 
-\item{fixseed}{A logical value indicating whether to fix the random seed 
-for generating the pemurations. If `TRUE`, the function will fix the random
-The default is `FALSE`.}
-
 \item{...}{Additional arguments passed to \code{corpcor::pcor.shrink},
 if "pcor.shrink" metric is selected.}
 }
diff --git a/man/runGraflex.Rd b/man/runGraflex.Rd
index 3184dd1..fe18119 100644
--- a/man/runGraflex.Rd
+++ b/man/runGraflex.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/2d-propd-graflex.R
+% Please edit documentation in R/3-shared-graflex.R
 \name{runGraflex}
 \alias{runGraflex}
-\title{Permute FDR for Multiple Concepts}
+\title{Calculate odds ratio and FDR}
 \usage{
-runGraflex(A, K, p = 500, fixseed = FALSE)
+runGraflex(A, K, p = 100, ncores = 1)
 }
 \arguments{
 \item{A}{An adjacency matrix.}
@@ -12,18 +12,15 @@ runGraflex(A, K, p = 500, fixseed = FALSE)
 \item{K}{A knowledge database where each row is a graph node
 and each column is a concept.}
 
-\item{p}{An integer. The number of overlaps to permute.}
-
-\item{fixseed}{A logical. If \code{TRUE}, the seed is fixed
-for reproducibility.}
+\item{p}{An integer. The number of permutation.}
 }
 \description{
-This function calls \code{\link{permuteOR}} for each
- concept (i.e., column) in the database \code{K}.
+This function calls \code{\link{graflex}} for each
+concept (i.e., column) in the database \code{K}.
 }
 \details{
-For each concept, this function calculates the
- false discovery rate (FDR) by counting the number of
- times the actual OR was greater than
- (or less than) a permuted OR.
+For each concept, this function calculates the odds ratio
+and determines the false discovery rate (FDR) by counting 
+the number of the actual OR was greater or less than a 
+permuted OR.
 }
diff --git a/man/selectReference.Rd b/man/selectReference.Rd
index e2749f4..e7ee62f 100644
--- a/man/selectReference.Rd
+++ b/man/selectReference.Rd
@@ -15,8 +15,9 @@ It can take the following values:
  - "clr" (default): Centered log-ratio transformation.
  - "alr": Additive log-ratio transformation ("pcor.bshrink" metric only).
  - "iqlr": Inter-quartile log-ratio transformation from ALDEx2.
- - The explicit name(s) of variable(s) to use as a reference.
-Use NA to skip log-ratio transformation.}
+ - The explicit name(s) or index(es) of variable(s) to use as a reference.
+ - Use NA to skip log-ratio transformation and any other pre-processing, like
+zero replacement. This is useful when the input data is already pre-processed.}
 
 \item{alpha}{The alpha parameter used in the alpha log-ratio transformation.}
 }
diff --git a/man/simple_zero_replacement.Rd b/man/simple_zero_replacement.Rd
index 451b737..89f2015 100644
--- a/man/simple_zero_replacement.Rd
+++ b/man/simple_zero_replacement.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/1a-propr-backend.R
+% Please edit documentation in R/3-shared-zero.R
 \name{simple_zero_replacement}
 \alias{simple_zero_replacement}
 \title{Simple Zero Replacement in a Count Matrix}
@@ -7,11 +7,11 @@
 simple_zero_replacement(ct)
 }
 \arguments{
-\item{ct}{A data matrix for which the log-ratio transformation will be
-performed. It is assumed that the matrix contains numerical values only.}
+\item{ct}{A data matrix containing numerical values.}
 }
 \value{
 A matrix with zero values replaced by the next smallest non-zero value.
+If no zeros are found, the function returns the original matrix.
 }
 \description{
 This function replaces zeros with the next smallest non-zero value in the
@@ -21,5 +21,4 @@ This function replaces zeros with the next smallest non-zero value in the
 \examples{
 # Sample input count data with zeros
 data <- matrix(c(0, 2, 3, 4, 5, 0), nrow = 2, byrow = TRUE)
-
 }
diff --git a/man/updateCutoffs.Rd b/man/updateCutoffs.Rd
index 990bee3..13b119b 100644
--- a/man/updateCutoffs.Rd
+++ b/man/updateCutoffs.Rd
@@ -6,17 +6,31 @@
 \alias{updateCutoffs.propd}
 \title{Update FDR by Permutation}
 \usage{
-updateCutoffs(object, cutoff = seq(0.05, 0.95, 0.3), ncores = 1)
+updateCutoffs(
+  object,
+  number_of_cutoffs = 100,
+  custom_cutoffs = NULL,
+  tails = NULL,
+  ncores = 1
+)
 
-updateCutoffs.propr(object, cutoff, ncores)
+updateCutoffs.propr(object, cutoffs, ncores)
 
-updateCutoffs.propd(object, cutoff = seq(0.05, 0.95, 0.3))
+updateCutoffs.propd(object, cutoffs, ncores)
 }
 \arguments{
 \item{object}{A \code{propr} or \code{propd} object.}
 
-\item{cutoff}{For \code{updateCutoffs}, a numeric vector.
-this argument provides the FDR cutoffs to test.}
+\item{number_of_cutoffs}{An integer. The number of cutoffs to test.
+Given this number, the cutoffs will be determined based on the quantile of the data.
+In this way, the cutoffs will be evenly spaced across the data.}
+
+\item{custom_cutoffs}{A numeric vector. When provided, this vector is used
+as the FDR cutoffs to test, and number_of_cutoffs is ignored.}
+
+\item{tails}{1 for one-sided on the right, and 2 for two-sided. When two-sided,
+the FDR is calculated for both positive and negative values. When NULL, use default
+option for the given metric. This is only relevant for \code{propr} objects.}
 
 \item{ncores}{An integer. The number of parallel cores to use.}
 }
@@ -33,7 +47,8 @@ This function wraps \code{updateCutoffs.propr} and
 \section{Methods}{
 
 \code{updateCutoffs.propr:}
- Use the \code{propr} object to permute proportionality
+ Use the \code{propr} object to permute correlation-like metrics
+ (ie. rho, phi, phs, cor, pcor, pcor.shrink, pcor.bshrink),
  across a number of cutoffs. Since the permutations get saved
  when the object is created, calling \code{updateCutoffs}
  will use the same random seed each time.
diff --git a/man/updatePermutes.Rd b/man/updatePermutes.Rd
index 042a0ad..a2b0cbb 100644
--- a/man/updatePermutes.Rd
+++ b/man/updatePermutes.Rd
@@ -4,12 +4,12 @@
 \alias{updatePermutes}
 \title{Create permuted data}
 \usage{
-updatePermutes(object, p, fixseed = FALSE)
+updatePermutes(object, p = 100)
 }
 \arguments{
 \item{object}{A \code{propr} or \code{propd} object.}
 
-\item{p}{The number of permutations to perform.}
+\item{p}{The number of permutations to perform. Default is 100.}
 }
 \value{
 A \code{propr} or \code{propd} object with the permutes slot updated.
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 45e951f..1e4572b 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -221,6 +221,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// vector2mat
+NumericMatrix vector2mat(NumericVector X, IntegerVector i, IntegerVector j);
+RcppExport SEXP _propr_vector2mat(SEXP XSEXP, SEXP iSEXP, SEXP jSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type X(XSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type i(iSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type j(jSEXP);
+    rcpp_result_gen = Rcpp::wrap(vector2mat(X, i, j));
+    return rcpp_result_gen;
+END_RCPP
+}
 // ratiosRcpp
 NumericMatrix ratiosRcpp(NumericMatrix& X);
 RcppExport SEXP _propr_ratiosRcpp(SEXP XSEXP) {
@@ -232,6 +245,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// results2matRcpp
+NumericMatrix results2matRcpp(DataFrame& results, int n, double diagonal);
+RcppExport SEXP _propr_results2matRcpp(SEXP resultsSEXP, SEXP nSEXP, SEXP diagonalSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< DataFrame& >::type results(resultsSEXP);
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< double >::type diagonal(diagonalSEXP);
+    rcpp_result_gen = Rcpp::wrap(results2matRcpp(results, n, diagonal));
+    return rcpp_result_gen;
+END_RCPP
+}
 // count_less_than
 int count_less_than(NumericVector x, double cutoff);
 RcppExport SEXP _propr_count_less_than(SEXP xSEXP, SEXP cutoffSEXP) {
@@ -256,6 +282,30 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// count_less_equal_than
+int count_less_equal_than(NumericVector x, double cutoff);
+RcppExport SEXP _propr_count_less_equal_than(SEXP xSEXP, SEXP cutoffSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type cutoff(cutoffSEXP);
+    rcpp_result_gen = Rcpp::wrap(count_less_equal_than(x, cutoff));
+    return rcpp_result_gen;
+END_RCPP
+}
+// count_greater_equal_than
+int count_greater_equal_than(NumericVector x, double cutoff);
+RcppExport SEXP _propr_count_greater_equal_than(SEXP xSEXP, SEXP cutoffSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type cutoff(cutoffSEXP);
+    rcpp_result_gen = Rcpp::wrap(count_greater_equal_than(x, cutoff));
+    return rcpp_result_gen;
+END_RCPP
+}
 // ctzRcpp
 NumericVector ctzRcpp(NumericMatrix& X);
 RcppExport SEXP _propr_ctzRcpp(SEXP XSEXP) {
@@ -267,6 +317,102 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// get_lower_triangle
+IntegerVector get_lower_triangle(IntegerMatrix& mat);
+RcppExport SEXP _propr_get_lower_triangle(SEXP matSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP);
+    rcpp_result_gen = Rcpp::wrap(get_lower_triangle(mat));
+    return rcpp_result_gen;
+END_RCPP
+}
+// shuffle_and_get_lower_triangle
+IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat);
+RcppExport SEXP _propr_shuffle_and_get_lower_triangle(SEXP matSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP);
+    rcpp_result_gen = Rcpp::wrap(shuffle_and_get_lower_triangle(mat));
+    return rcpp_result_gen;
+END_RCPP
+}
+// binTab
+NumericMatrix binTab(IntegerVector& A, IntegerVector& G);
+RcppExport SEXP _propr_binTab(SEXP ASEXP, SEXP GSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP);
+    rcpp_result_gen = Rcpp::wrap(binTab(A, G));
+    return rcpp_result_gen;
+END_RCPP
+}
+// getOR
+NumericMatrix getOR(IntegerVector& A, IntegerVector& G);
+RcppExport SEXP _propr_getOR(SEXP ASEXP, SEXP GSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP);
+    rcpp_result_gen = Rcpp::wrap(getOR(A, G));
+    return rcpp_result_gen;
+END_RCPP
+}
+// permuteOR
+NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p);
+RcppExport SEXP _propr_permuteOR(SEXP ASEXP, SEXP GstarSEXP, SEXP pSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< IntegerVector& >::type Gstar(GstarSEXP);
+    Rcpp::traits::input_parameter< int >::type p(pSEXP);
+    rcpp_result_gen = Rcpp::wrap(permuteOR(A, Gstar, p));
+    return rcpp_result_gen;
+END_RCPP
+}
+// getFDR_over
+double getFDR_over(double actual, NumericVector permuted);
+RcppExport SEXP _propr_getFDR_over(SEXP actualSEXP, SEXP permutedSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< double >::type actual(actualSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP);
+    rcpp_result_gen = Rcpp::wrap(getFDR_over(actual, permuted));
+    return rcpp_result_gen;
+END_RCPP
+}
+// getFDR_under
+double getFDR_under(double actual, NumericVector permuted);
+RcppExport SEXP _propr_getFDR_under(SEXP actualSEXP, SEXP permutedSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< double >::type actual(actualSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP);
+    rcpp_result_gen = Rcpp::wrap(getFDR_under(actual, permuted));
+    return rcpp_result_gen;
+END_RCPP
+}
+// graflex
+NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p);
+RcppExport SEXP _propr_graflex(SEXP ASEXP, SEXP GSEXP, SEXP pSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< IntegerMatrix& >::type G(GSEXP);
+    Rcpp::traits::input_parameter< int >::type p(pSEXP);
+    rcpp_result_gen = Rcpp::wrap(graflex(A, G, p));
+    return rcpp_result_gen;
+END_RCPP
+}
 // lr2vlr
 NumericMatrix lr2vlr(NumericMatrix lr);
 RcppExport SEXP _propr_lr2vlr(SEXP lrSEXP) {
@@ -385,10 +531,22 @@ static const R_CallMethodDef CallEntries[] = {
     {"_propr_urtRcpp", (DL_FUNC) &_propr_urtRcpp, 1},
     {"_propr_labRcpp", (DL_FUNC) &_propr_labRcpp, 1},
     {"_propr_half2mat", (DL_FUNC) &_propr_half2mat, 1},
+    {"_propr_vector2mat", (DL_FUNC) &_propr_vector2mat, 3},
     {"_propr_ratiosRcpp", (DL_FUNC) &_propr_ratiosRcpp, 1},
+    {"_propr_results2matRcpp", (DL_FUNC) &_propr_results2matRcpp, 3},
     {"_propr_count_less_than", (DL_FUNC) &_propr_count_less_than, 2},
     {"_propr_count_greater_than", (DL_FUNC) &_propr_count_greater_than, 2},
+    {"_propr_count_less_equal_than", (DL_FUNC) &_propr_count_less_equal_than, 2},
+    {"_propr_count_greater_equal_than", (DL_FUNC) &_propr_count_greater_equal_than, 2},
     {"_propr_ctzRcpp", (DL_FUNC) &_propr_ctzRcpp, 1},
+    {"_propr_get_lower_triangle", (DL_FUNC) &_propr_get_lower_triangle, 1},
+    {"_propr_shuffle_and_get_lower_triangle", (DL_FUNC) &_propr_shuffle_and_get_lower_triangle, 1},
+    {"_propr_binTab", (DL_FUNC) &_propr_binTab, 2},
+    {"_propr_getOR", (DL_FUNC) &_propr_getOR, 2},
+    {"_propr_permuteOR", (DL_FUNC) &_propr_permuteOR, 3},
+    {"_propr_getFDR_over", (DL_FUNC) &_propr_getFDR_over, 2},
+    {"_propr_getFDR_under", (DL_FUNC) &_propr_getFDR_under, 2},
+    {"_propr_graflex", (DL_FUNC) &_propr_graflex, 3},
     {"_propr_lr2vlr", (DL_FUNC) &_propr_lr2vlr, 1},
     {"_propr_lr2phi", (DL_FUNC) &_propr_lr2phi, 1},
     {"_propr_lr2rho", (DL_FUNC) &_propr_lr2rho, 1},
diff --git a/src/backend.cpp b/src/backend.cpp
index a55a9e1..69f723a 100755
--- a/src/backend.cpp
+++ b/src/backend.cpp
@@ -427,6 +427,27 @@ NumericMatrix half2mat(NumericVector X){
   return A;
 }
 
+// Function to recast a vector into matrix, given the positions i and j
+// [[Rcpp::export]]
+NumericMatrix vector2mat(NumericVector X, IntegerVector i, IntegerVector j){
+  
+  int nfeats = sqrt(2 * X.length() + .25) + .5;
+  NumericMatrix A(nfeats, nfeats);
+
+  int ni = i.length();
+  int nj = j.length();
+  if (ni != nj){
+    stop("i and j must be the same length.");
+  }
+
+  for (int counter = 0; counter < ni; counter++){
+    A(i[counter]-1, j[counter]-1) = X[counter];   // NOTE that R is 1-indexed whereas C is 0-indexed
+    A(j[counter]-1, i[counter]-1) = X[counter];
+  }
+
+  return A;  
+}
+
 // Function to recast matrix as feature ratios
 // [[Rcpp::export]]
 NumericMatrix ratiosRcpp(NumericMatrix & X){
@@ -446,3 +467,25 @@ NumericMatrix ratiosRcpp(NumericMatrix & X){
 
   return result;
 }
+
+// Function to recast results data frame as gene-gene matrix
+// [[Rcpp::export]]
+NumericMatrix results2matRcpp(DataFrame& results, int n, double diagonal = 0.0){
+  NumericMatrix mat(n, n);
+
+  // fill in the values for each pair
+  int npairs = results.nrows();
+  for (int i = 0; i < npairs; i++){
+    int row = int(results(i, 0)) - 1;
+    int col = int(results(i, 1)) - 1;
+    mat(row, col) = results(i, 2);
+    mat(col, row) = results(i, 2);
+  }
+
+  // fill in the diagonal values
+  for (int i = 0; i < n; i++){
+    mat(i, i) = diagonal;
+  }
+  
+  return mat;
+}
\ No newline at end of file
diff --git a/src/comparison.cpp b/src/comparison.cpp
index 4fa8aed..d154cac 100644
--- a/src/comparison.cpp
+++ b/src/comparison.cpp
@@ -36,3 +36,33 @@ int count_greater_than(NumericVector x, double cutoff) {
 
   return count;
 }
+
+// Count the number of elements in vector `x` that are less or equal than the
+// `cutoff`.
+//
+// [[Rcpp::export]]
+int count_less_equal_than(NumericVector x, double cutoff) {
+  int count = 0;
+  int len = x.size();
+
+  for (int i = 0; i < len; ++i) {
+    count += x[i] <= cutoff;
+  }
+
+  return count;
+}
+
+// Count the number of elements in vector `x` that are greater or equal than the
+// `cutoff`.
+//
+// [[Rcpp::export]]
+int count_greater_equal_than(NumericVector x, double cutoff) {
+  int count = 0;
+  int len = x.size();
+
+  for (int i = 0; i < len; ++i) {
+    count += x[i] >= cutoff;
+  }
+
+  return count;
+}
\ No newline at end of file
diff --git a/src/graflex.cpp b/src/graflex.cpp
new file mode 100644
index 0000000..255e6e3
--- /dev/null
+++ b/src/graflex.cpp
@@ -0,0 +1,137 @@
+#include <Rcpp.h>
+#include <numeric>
+using namespace Rcpp;
+
+
+// Function to extract the triangle of a square and symmetric IntegerMatrix
+// [[Rcpp::export]]
+IntegerVector get_lower_triangle(IntegerMatrix& mat) {
+  int nrow = mat.nrow();
+  int ncol = mat.ncol();
+  int n = ncol * (ncol - 1) / 2;
+
+  IntegerVector triangle(n);
+
+  int k = 0;
+  for (int j = 0; j < ncol; ++j) {
+      for (int i = j+1; i < nrow; ++i) {
+      triangle[k++] = mat(i, j);
+    }
+  }
+
+  return triangle;
+}
+
+// Function to shuffle a square and symmetric IntegerMatrix, and get the lower triangle
+// [[Rcpp::export]]
+IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat) {
+  int nrow = mat.nrow();
+  int ncol = mat.ncol();
+  int n = ncol * (ncol - 1) / 2;
+
+  IntegerVector shuffled_triangle(n);
+  IntegerVector index = sample(ncol, ncol, false) - 1;
+
+  int k = 0;
+  for (int i = 0; i < ncol; ++i) {
+    int index_i = index[i];
+    for (int j = 0; j < i; ++j) {
+      int index_j = index[j];
+      shuffled_triangle[k++] = mat(index_i, index_j);
+    }
+  }
+
+  return shuffled_triangle;
+}
+
+// Function to calculate the contingency table
+// [[Rcpp::export]]
+NumericMatrix binTab(IntegerVector& A, IntegerVector& G) {
+  NumericMatrix tab(2, 2);
+  int n = A.size();
+
+  tab(0, 0) = sum((1 - A) * (1 - G));  // not in A and not in G
+  tab(0, 1) = sum((1 - A) * G);        // not in A but in G
+  tab(1, 0) = sum(A * (1 - G));        // in A but not G
+  tab(1, 1) = sum(A * G);              // in A and G
+
+  return tab;
+}
+
+// Function to calculate the odds ratio, and parse all information into a matrix
+// [[Rcpp::export]]
+NumericMatrix getOR(IntegerVector& A, IntegerVector& G) {
+
+  NumericMatrix tab = binTab(A, G);
+  double odds_ratio = (tab(0, 0) * tab(1, 1)) / (tab(0, 1) * tab(1, 0));
+
+  NumericMatrix or_table(1, 8);
+  or_table(0, 0) = tab(0, 0);
+  or_table(0, 1) = tab(0, 1);
+  or_table(0, 2) = tab(1, 0);
+  or_table(0, 3) = tab(1, 1);
+  or_table(0, 4) = odds_ratio;
+  or_table(0, 5) = log(odds_ratio);
+
+  return or_table;
+}
+
+// Function to calculate the odds ratio and other relevant info for each permutation
+// [[Rcpp::export]]
+NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p = 100) {
+
+  NumericMatrix or_table(p, 8);
+
+  // calculate odds ratio for each permutation
+  for (int i = 0; i < p; ++i) {
+    IntegerVector Astar = shuffle_and_get_lower_triangle(A);
+    or_table(i, _) = getOR(Astar, Gstar);
+  }
+
+  return(or_table);
+}
+
+// [[Rcpp::export]]
+double getFDR_over(double actual, NumericVector permuted) {
+  double n = permuted.size();
+  double fdr = 0.0;
+  for (int i = 0; i < n; ++i) {
+    if (permuted[i] >= actual) {
+      fdr += 1.0;
+    }
+  }
+  fdr /= n;
+  return fdr;
+}
+
+// [[Rcpp::export]]
+double getFDR_under(double actual, NumericVector permuted) {
+  double n = permuted.size();
+  double fdr = 0.0;
+  for (int i = 0; i < n; ++i) {
+    if (permuted[i] <= actual) {
+      fdr += 1.0;
+    }
+  }
+  fdr /= n;
+  return fdr;
+}
+
+// Function to calculate the odds ratio and FDR, given the adjacency matrix A and the knowledge graph G
+// [[Rcpp::export]]
+NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p = 100) {
+
+  // get the actual odds ratio
+  IntegerVector Astar = get_lower_triangle(A);
+  IntegerVector Gstar = get_lower_triangle(G);
+  NumericMatrix actual = getOR(Astar, Gstar);
+
+  // get distribution of odds ratios on permuted data
+  NumericMatrix permuted = permuteOR(A, Gstar, p);
+
+  // calculate the FDR
+  actual(0, 6) = getFDR_under(actual(0, 4), permuted(_, 4));
+  actual(0, 7) = getFDR_over(actual(0, 4), permuted(_, 4));
+
+  return actual;
+}
diff --git a/tests/testthat/test-GET-getMatrix.R b/tests/testthat/test-GET-getMatrix.R
new file mode 100644
index 0000000..97f9578
--- /dev/null
+++ b/tests/testthat/test-GET-getMatrix.R
@@ -0,0 +1,49 @@
+library(testthat)
+library(propr)
+
+test_that('get propd matrix work',{
+
+  # compute differential proportionality
+  x <- iris[,1:4]  # data matrix with 4 variables
+  y <- iris[,5]    # group vector
+  pd <- propd(x, as.character(y))
+  mat <- getMatrix(pd)
+
+  # check matrix is square and have 4 columns
+  expect_equal(ncol(mat), nrow(mat))
+  expect_equal(ncol(mat), 4)
+
+  # check diagonal is zero
+  expect_equal(mat[1,1], 0)
+  expect_equal(mat[2,2], 0)
+
+  # check matrix values
+  results <- pd@results
+  expect_equal(
+    mat[1,2], 
+    results[which(results$Partner==2 & results$Pair==1),3]
+  )
+  expect_equal(
+    mat[2,3], 
+    results[which(results$Partner==3 & results$Pair==2),3]
+  )
+})
+
+test_that('get propr matrix work', {
+
+  # create matrix and compute proportionality
+  N <- 100
+  a <- seq(from = 5, to = 15, length.out = N)
+  b <- a * rnorm(N, mean = 1, sd = 0.1)
+  c <- rnorm(N, mean = 10)
+  d <- rnorm(N, mean = 10)
+  e <- rep(10, N)
+  X <- data.frame(a, b, c, d, e)
+  pr <- propr(X, metric = "rho")
+
+  # check
+  expect_equal(
+    pr@matrix,
+    getMatrix(pr)
+  )
+})
\ No newline at end of file
diff --git a/tests/testthat/test-GRAFLEX.R b/tests/testthat/test-GRAFLEX.R
new file mode 100644
index 0000000..e06fd55
--- /dev/null
+++ b/tests/testthat/test-GRAFLEX.R
@@ -0,0 +1,217 @@
+library(testthat)
+library(propr)
+library(Rcpp)
+
+# ===================== #
+# old code for graflex  #
+# ===================== #
+
+permuteOR_old <- function(A, G, p = 500) {
+  Gstar <- G[lower.tri(G)]
+  res <- lapply(1:p, function(i) {
+    # Shuffle the adjacency matrix
+    # index <- sample(1:ncol(A))
+    # A <- A[index, index]
+    # Astar <- A[lower.tri(A)]
+    Astar <- propr:::shuffle_and_get_lower_triangle(A)
+    getOR_old(Astar, Gstar)
+  })
+
+  do.call("rbind", res)
+}
+
+binTab_old <- function(A, G) {
+  diff <- A != G
+  only1 <- A[diff]
+  b <- as.numeric(sum(only1))
+  c <- as.numeric(length(only1) - b)
+
+  same <- !diff
+  double1 <- A[same]
+  a <- as.numeric(sum(double1))
+  d <- as.numeric(length(double1) - a)
+
+  matrix(c(d, b, c, a), 2, 2)
+}
+
+getOR_old <- function(A, G) {
+  tab <- binTab_old(A, G)
+  or <- (tab[1, 1] * tab[2, 2]) / (tab[1, 2] * tab[2, 1])
+  data.frame(
+    "Neither" = tab[1, 1],
+    "G.only" = tab[1, 2],
+    "A.only" = tab[2, 1],
+    "Both" = tab[2, 2],
+    "Odds" = or,
+    "LogOR" = log(or)
+  )
+}
+
+calculateOR_old <- function(A, G) {
+  Astar <- A[lower.tri(A)]
+  Gstar <- G[lower.tri(G)]
+  getOR_old(Astar, Gstar)
+}
+
+getFDR_old <- function(actual, permuted) {
+  actual$FDR.under <-
+    sum(permuted$Odds <= actual$Odds) / nrow(permuted)
+  actual$FDR.over <-
+    sum(permuted$Odds >= actual$Odds) / nrow(permuted)
+  actual
+}
+
+runGraflex_old <- function(A, K, p = 500) {
+  if (nrow(A) != nrow(K))
+    stop("'A' and 'K' must have identical rows.")
+
+  numTicks <- 0
+  res <- lapply(1:ncol(K), function(k) {
+    Gk <- K[, k] %*% t(K[, k])
+    actual <- calculateOR_old(A, Gk)
+    permuted <- permuteOR_old(A, Gk, p = p)
+    actual <- getFDR_old(actual, permuted)
+    actual$Permutes <- p
+    actual$Concept <- colnames(K)[k]
+    actual
+  })
+
+  do.call("rbind", res)
+}
+
+# ===================== #
+# run tests             #
+# ===================== #
+
+
+test_that("check if odds ratio is computed correctly", {
+
+  # Create two vectors of 1 and 0
+  Astar <- c(1, 0, 1, 1, 0, 1, 0, 0, 1, 0)
+  Gstar <- c(1, 0, 0, 1, 1, 0, 0, 1, 0, 1)
+  
+  # Expected values
+  a <- 2  # not in A not in G
+  b <- 3  # not in A but in G
+  c <- 3  # in A but not in G
+  d <- 2  # in A in G
+  expected_odds_ratio <- (a * d) / (b * c)
+
+  # compute odds ratio table
+  or_table <- propr:::getOR(Astar, Gstar)
+
+  # check
+  expect_equal(a, or_table[1,1])
+  expect_equal(b, or_table[1,2])
+  expect_equal(c, or_table[1,3])
+  expect_equal(d, or_table[1,4])
+  expect_equal(expected_odds_ratio, or_table[1,5])
+})
+
+test_that("check if runGraflex produce the expected results", {
+
+  # Create matrices of 0 and 1
+  A <- matrix(c(1, 1, 0, 1, 0, 
+                1, 1, 1, 0, 1, 
+                0, 1, 1, 0, 1, 
+                1, 0, 0, 1, 1, 
+                0, 1, 1, 1, 1), 
+              nrow = 5, byrow = TRUE)
+  K <- matrix(c(1, 0, 
+                0, 1, 
+                1, 0, 
+                0, 1, 
+                1, 1), 
+              nrow = 5, byrow = TRUE)
+  colnames(K) <- c("C1", "C2")
+  
+  # Expected values for C1
+  a1 <- 2  # not in A not in G
+  b1 <- 2  # not in A but in G
+  c1 <- 5  # in A but not in G
+  d1 <- 1  # in A in G
+  expected_odds_ratio1 <- (a1 * d1) / (b1 * c1)
+
+  # Expected values for C2
+  a2 <- 3  # not in A not in G
+  b2 <- 1  # not in A but in G
+  c2 <- 4  # in A but not in G
+  d2 <- 2  # in A in G
+  expected_odds_ratio2 <- (a2 * d2) / (b2 * c2)
+
+  # compute graflex
+  res <- propr:::runGraflex(A, K)
+
+  # check
+  expect_equal(a1, as.numeric(res[1,1]))
+  expect_equal(b1, as.numeric(res[1,2]))
+  expect_equal(c1, as.numeric(res[1,3]))
+  expect_equal(d1, as.numeric(res[1,4]))
+  expect_equal(expected_odds_ratio1, as.numeric(res[1,5]))
+  expect_equal(a2, as.numeric(res[2,1]))
+  expect_equal(b2, as.numeric(res[2,2]))
+  expect_equal(c2, as.numeric(res[2,3]))
+  expect_equal(d2, as.numeric(res[2,4]))
+  expect_equal(expected_odds_ratio2, as.numeric(res[2,5]))
+})
+
+test_that("check reproducibility seed works", {
+  
+    # Create a matrix of 0 and 1
+    A <- matrix(c(1, 1, 0, 1, 0, 
+                  1, 1, 1, 0, 1, 
+                  0, 1, 1, 0, 1, 
+                  1, 0, 0, 1, 1, 
+                  0, 1, 1, 1, 1), 
+                nrow = 5, byrow = TRUE)
+    K <- matrix(c(1, 0, 
+                  0, 1, 
+                  1, 0, 
+                  0, 1, 
+                  1, 1), 
+                nrow = 5, byrow = TRUE)
+    colnames(K) <- c("C1", "C2")
+    
+    # compute graflex
+    set.seed(0)
+    res1 <- propr:::runGraflex(A, K, p=100)
+    set.seed(0)
+    res2 <- propr:::runGraflex(A, K, p=100)
+    set.seed(123)
+    res3 <- propr:::runGraflex(A, K, p=100)
+  
+    # check
+    expect_equal(res1, res2)
+    expect_false(isTRUE(all.equal(res1, res3)))
+})
+
+test_that("check if same results are obtained compared to old code", {
+    
+    # Create matrices of 0 and 1
+    A <- matrix(c(1, 1, 0, 1, 0, 
+                  1, 1, 1, 0, 1, 
+                  0, 1, 1, 0, 1, 
+                  1, 0, 0, 1, 1, 
+                  0, 1, 1, 1, 1), 
+                nrow = 5, byrow = TRUE)
+    K <- matrix(c(1, 0, 
+                  0, 1, 
+                  1, 0, 
+                  0, 1, 
+                  1, 1), 
+                nrow = 5, byrow = TRUE)
+    colnames(K) <- c("C1", "C2")
+
+    # compute graflex
+    set.seed(0)
+    res <- propr:::runGraflex(A, K, p=100)
+    res$LogOR <- round(res$LogOR, 5)
+    
+    # compute graflex with old code
+    set.seed(0)
+    res_old <- runGraflex_old(A, K, p=100)
+    res_old$LogOR <- round(res_old$LogOR, 5)
+    
+    # check
+    expect_equal(res, res_old)
+})
diff --git a/tests/testthat/test-PROPR-backend.R b/tests/testthat/test-PROPR-backend.R
index 664fd88..2864c05 100644
--- a/tests/testthat/test-PROPR-backend.R
+++ b/tests/testthat/test-PROPR-backend.R
@@ -1,27 +1,53 @@
 library(testthat)
 library(propr)
 
-test_that("simple_zero_replacement replaces zeros correctly", {
-  # Example data with zeros
-  data_with_zeros <- matrix(c(1, 2, 3, 0, 5, 6, 0, 8, 9), nrow = 3, byrow = TRUE)
+# define data
+data <- matrix(c(1:12), ncol=4)
+colnames(data) <- c("A", "B", "C", "D")
 
-  # Apply the simple_zero_replacement function to the example data
-  replaced_data <- simple_zero_replacement(data_with_zeros)
 
-  # Define the expected output after zero replacement
-  expected_output <- matrix(c(1, 2, 3, 1, 5, 6, 1, 8, 9), nrow = 3, byrow = TRUE)
+test_that("index_reference works properly", {
+
+  expect_identical(
+    as.integer(index_reference(data, 2)),
+    as.integer(2)
+  )
+  expect_identical(
+    as.integer(index_reference(data, c("A", "C"))),
+    as.integer(c(1, 3))
+  )
+  expect_identical(
+    as.integer(index_reference(data, "clr")),
+    as.integer(c(1, 2, 3, 4))
+  )
 
-  # Compare the output with the expected output
-  expect_identical(replaced_data, expected_output)
 })
 
 test_that("logratio_without_alpha performs log-ratio transformation correctly", {
-  # Example data
-  example_data <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 3, byrow = TRUE)
 
-  # Apply the logratio_without_alpha function to the example data
-  transformed_data <- logratio_without_alpha(example_data, use = 2)
+  # when ivar is 2
+  transformed_data <- logratio_without_alpha(data, use = 2)
+  expect_identical(
+    as.numeric(round(transformed_data[1,1], 6)), 
+    round(log(1/4), 6)
+  )
+  expect_identical(
+    as.numeric(round(transformed_data[2,1], 6)),
+    round(log(2/5), 6)
+  )
+
+  # when ivar is 1,3
+  transformed_data <- logratio_without_alpha(data, use = c(1,3))
+  expected <- t(apply(data, 1, function(x) log(x) - mean(log(x[c(1,3)]))))
+  expect_true(
+    all(round(transformed_data, 6) == round(expected, 6))
+  )
+
+  # when ivar is clr
+  transformed_data <- logratio_without_alpha(data, use = c(1,2,3,4))
+  expected <- t(apply(data, 1, function(x) log(x) - mean(log(x))))
+  expect_true(
+    all(round(transformed_data, 6) == round(expected, 6))
+  )
 
-  # Compare the output with the expected output
-  expect_identical(transformed_data[1,1], log(1/2))
 })
diff --git a/tests/testthat/test-PROPR-ivar.R b/tests/testthat/test-PROPR-ivar.R
new file mode 100644
index 0000000..4ccf2d1
--- /dev/null
+++ b/tests/testthat/test-PROPR-ivar.R
@@ -0,0 +1,153 @@
+library(testthat)
+library(propr)
+
+data(mtcars)
+
+test_that("the data is properly handled when ivar is NA", {
+
+    # compute propr object
+    pr <- propr(mtcars, ivar=NA)
+
+    # check the counts remain the same as the original input
+    expect_equal(
+        pr@counts, 
+        mtcars
+    )
+
+    # check the logratios are the same
+    expect_equal(
+        pr@logratio, 
+        mtcars
+    )
+})
+
+test_that("the data is properly handled when ivar is clr", {
+
+    # compute propr object
+    pr <- propr(mtcars, ivar="clr")
+
+    # compute expected data
+    ct <- simple_zero_replacement(mtcars)
+    clr <- logratio_without_alpha(ct, c(1:ncol(ct)))
+
+    # check the counts remain the same as the zero-replaced input
+    expect_equal(
+        pr@counts, 
+        ct
+    )
+
+    # check the logratios are the same
+    expect_equal(
+        pr@logratio, 
+        clr
+    )
+})
+
+test_that("the data is properly handled when ivar is 3", {
+
+    # compute propr object
+    pr <- propr(mtcars, ivar=3)
+
+    # compute expected data
+    ct <- simple_zero_replacement(mtcars)
+    alr <- logratio_without_alpha(ct, 3)
+
+    # check the counts remain the same as the original input
+    expect_equal(
+        pr@counts, 
+        ct
+    )
+
+    # check the logratios are the same
+    expect_equal(
+        pr@logratio, 
+        alr
+    )
+})
+
+test_that("the data is properly handled when ivar is 1,6", {
+
+    # compute propr object
+    pr <- propr(mtcars, ivar=c(1,6))
+
+    # compute expected data
+    ct <- simple_zero_replacement(mtcars)
+    alr <- logratio_without_alpha(ct, c(1,6))
+
+    # check the counts remain the same as the original input
+    expect_equal(
+        pr@counts, 
+        ct
+    )
+
+    # check the logratios are the same
+    expect_equal(
+        pr@logratio, 
+        alr
+    )
+})
+
+test_that("pearson correlation is correct when ivar is NA", {
+
+    # compute propr object
+    cor_propr <- propr(mtcars, metric = "cor", ivar=NA)@matrix
+    
+    # get correlation using cor
+    cor_cor <- cor(mtcars, method = "pearson")
+
+    expect_equal(
+        round(cor_propr, 6), 
+        round(cor_cor, 6)
+    )
+})
+
+test_that("pearson correlation is correct when ivar is clr", {
+
+    # get correlation using propr
+    pr <- propr(mtcars, metric = "cor", ivar="clr")
+
+    # get correlation using cor
+    ct <- simple_zero_replacement(mtcars)
+    clr <- logratio_without_alpha(ct, c(1:ncol(ct)))
+    ccor <- cor(clr, method = "pearson")
+
+    # check the correlations are the same
+    expect_equal(
+        round(pr@matrix, 6), 
+        round(ccor, 6)
+    )
+})
+
+test_that("pearson correlation is correct when ivar is 1", {
+
+    # get correlation using propr
+    pr <- suppressWarnings(propr(mtcars, metric = "cor", ivar=1))
+
+    # get correlation using cor
+    ct <- simple_zero_replacement(mtcars)
+    alr <- logratio_without_alpha(ct, 1)
+    ccor <- suppressWarnings(cor(alr, method = "pearson"))
+
+    # check the correlations are the same
+    expect_equal(
+        round(pr@matrix, 6), 
+        round(ccor, 6)
+    )
+})
+
+test_that("pearson correlation is correct when ivar is 1,3", {
+
+    # get correlation using propr
+    pr <- suppressWarnings(propr(mtcars, metric = "cor", ivar=c(1,3)))
+
+    # get correlation using cor
+    ct <- simple_zero_replacement(mtcars)
+    alr <- logratio_without_alpha(ct, c(1,3))
+    ccor <- suppressWarnings(cor(alr, method = "pearson"))
+
+    # check the correlations are the same
+    expect_equal(
+        round(pr@matrix, 6), 
+        round(ccor, 6)
+    )
+})
diff --git a/tests/testthat/test-PROPR-pcor.R b/tests/testthat/test-PROPR-pcor.R
new file mode 100644
index 0000000..2d0437c
--- /dev/null
+++ b/tests/testthat/test-PROPR-pcor.R
@@ -0,0 +1,110 @@
+library(testthat)
+library(propr)
+
+# define data
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+for (i in 1:4){
+  X[sample(1:N, 10),i] <- 0
+}
+
+test_that("pcor is correct when ivar is given", {
+
+  d1 <- list("clr", 5, c(1,3))
+  d2 <- list(c(1:5), 5, c(1,3))
+
+  for (i in 1:length(d1)){
+
+    # compute pcor manually
+    ct <- simple_zero_replacement(X)
+    lr <- logratio_without_alpha(ct, d2[[i]])
+    cov <- cov(lr)
+    mat <- suppressWarnings(corpcor::cor2pcor(cov))
+
+    # compute pcor
+    pr <- suppressWarnings(propr(X, metric = "pcor", ivar=d1[[i]]))
+
+    # expect counts to have zeros replaced
+    expect_true(
+      all(pr@counts == ct)
+    )
+
+    # expect that logratio are equal
+    expect_true(
+      all(pr@logratio == lr)
+    )
+
+    # expect computed coefficients are equal
+    expect_true(
+      all(round(pr@matrix, 8) == round(mat, 8), na.rm=T)
+    )
+
+    # check shrinkage is not applied
+    expect_equal(
+      pr@lambda,
+      NULL
+    )
+  }
+})
+
+
+test_that("pcor is correct when ivar is NA using previously transformed data", {
+
+  d1 <- list("clr", 5, c(1,3))
+  d2 <- list(c(1:5), 5, c(1,3))
+
+  for (i in 1:length(d1)){
+
+    # compute pcor manually
+    ct <- simple_zero_replacement(X)
+    lr <- logratio_without_alpha(ct, d2[[i]])
+    cov <- cov(lr)
+    mat <- suppressWarnings(corpcor::cor2pcor(cov))
+
+    # compute pcor
+    pr <- suppressWarnings(propr(lr, metric = "pcor", ivar=NA))
+
+    # expect counts to contain the previously transformed data
+    expect_true(
+      all(pr@counts == lr)
+    )
+
+    # expect that the logratio slot also contains the exact input data
+    expect_true(
+      all(pr@logratio == lr)
+    )
+
+    # expect computed coefficients are equal
+    expect_true(
+      all(round(pr@matrix, 8) == round(mat, 8), na.rm=T)
+    )
+
+    # check shrinkage is not applied
+    expect_equal(
+      pr@lambda,
+      NULL
+    )
+  }
+})
+
+test_that("pcor with alr and clr are the same", {
+
+    # compute pcor with alr
+    pr <- suppressWarnings(propr(X, metric = "pcor", ivar=5))
+    pcor_alr <- getMatrix(pr)[1:4, 1:4]
+
+    # compute pcor with clr
+    pr <- propr(X, metric = "pcor", ivar="clr")
+    pcor_clr <- getMatrix(pr)[1:4, 1:4]
+
+    expect_equal(
+      round(pcor_alr, 8),
+      round(pcor_clr, 8)
+    )
+
+})
diff --git a/tests/testthat/test-PROPR-pcorbshrink.R b/tests/testthat/test-PROPR-pcorbshrink.R
index f320a7e..fc29c1b 100644
--- a/tests/testthat/test-PROPR-pcorbshrink.R
+++ b/tests/testthat/test-PROPR-pcorbshrink.R
@@ -1,6 +1,8 @@
 library(testthat)
 library(propr)
+library(corpcor)
 
+# define data
 N <- 100
 a <- seq(from = 5, to = 15, length.out = N)
 b <- a * rnorm(N, mean = 1, sd = 0.1)
@@ -8,20 +10,73 @@ c <- rnorm(N, mean = 10)
 d <- rnorm(N, mean = 10)
 e <- rep(10, N)
 X <- data.frame(a, b, c, d, e)
+for (i in 1:4){
+  X[sample(1:N, 10),i] <- 0
+}
 
-# compute pcor.bshrink with alr
-pr <- propr(X, metric = "pcor.bshrink", ivar="alr")
-pcor_alr <- getMatrix(pr)[1:4, 1:4]
+test_that("pcor.bshrink is correct when ivar is clr or alr",{
 
-# compute pcor.bshrink with clr
-pr <- propr(X, metric = "pcor.bshrink", ivar="clr")
-pcor_clr <- getMatrix(pr)[1:4, 1:4]
+  for (ivar in c("clr", "alr")){
+
+    # compute pcor manually
+    ct <- simple_zero_replacement(X)
+    out <- propr:::pcor.bshrink(ct, ivar)
+    mat <- out$matrix
+    lambda <- out$lambda
+
+    # compute pcor
+    pr <- propr(X, metric = "pcor.bshrink", ivar=ivar)
+
+    # expect counts to have zeros replaced
+    expect_true(
+      all(pr@counts == ct)
+    )
+
+    # NOTE that the data is not logratio transformed while computing pcor.bshrink
+    # it is internally handled by covariance conversion.
+    # so pr@logratio should be equal to pr@counts
+    expect_true(
+      all(pr@logratio == ct)
+    )
+  
+    # expect computed coefficients are equal
+    expect_true(
+      all(round(pr@matrix, 8) == round(mat, 8))
+    )
+
+    # expect same lambda
+    expect_equal(
+      pr@lambda,
+      lambda
+    )
+
+    # check dimensions are correct
+    expect_equal(ncol(pr@matrix), 5)
+    expect_equal(nrow(pr@matrix), 5)
+  }
+  
+})
+
+test_that("test that pcor.bshrink gives error when ivar is NA", {
+    expect_error(
+        propr(X, metric = "pcor.bshrink", ivar=NA)
+    )
+})
 
-# test that the results are the same
 test_that("pcor.bshrink with alr and clr are the same", {
 
-  expect_equal(
-    pcor_alr,
-    pcor_clr
-  )
+    # compute pcor with alr
+    pr_alr <- suppressWarnings(propr(X, metric = "pcor.bshrink", ivar='alr'))
+    pcor_alr <- getMatrix(pr_alr)[1:4, 1:4]
+
+    # compute pcor with clr
+    pr_clr <- propr(X, metric = "pcor.bshrink", ivar="clr")
+    pcor_clr <- getMatrix(pr_clr)[1:4, 1:4]
+
+    # expect that the coefficients are the same
+    expect_equal(
+        round(pcor_alr, 8),
+        round(pcor_clr, 8)
+    )
+
 })
diff --git a/tests/testthat/test-PROPR-pcorshrink.R b/tests/testthat/test-PROPR-pcorshrink.R
new file mode 100644
index 0000000..4f087e3
--- /dev/null
+++ b/tests/testthat/test-PROPR-pcorshrink.R
@@ -0,0 +1,141 @@
+library(testthat)
+library(propr)
+library(corpcor)
+
+# define data
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+for (i in 1:4){
+  X[sample(1:N, 10),i] <- 0
+}
+
+test_that("pcor.shrink is the same when ivar is 5, with or without the reference column", {
+
+  # compute lr
+  ct <- simple_zero_replacement(X)
+  lr <- logratio_without_alpha(ct, 5)
+
+  # compute pcor.shrink
+  pr1 <- suppressWarnings(propr(lr, metric = "pcor.shrink", ivar=NA))
+
+  # compute pcor.shrink without 5th column
+  pr2 <- suppressWarnings(propr(lr[1:4], metric = "pcor.shrink", ivar=NA))
+
+  # expect computed coefficients are equal
+  expect_true(
+    all(round(pr1@matrix, 8)[1:4,1:4] == round(pr2@matrix, 8)[1:4,1:4])
+  )
+
+  # check shrinkage is the same
+  expect_equal(
+    pr1@lambda,
+    pr2@lambda
+  )
+
+})
+
+test_that("pcor.shrink is correct when ivar is given", {
+
+  d1 <- list("clr", 5, c(1,3))
+  d2 <- list(c(1:5), 5, c(1,3))
+
+  for (i in 1:length(d1)){
+
+    # compute pcor manually
+    ct <- simple_zero_replacement(X)
+    lr <- logratio_without_alpha(ct, d2[[i]])
+    cov <- suppressWarnings(cov.shrink(lr))
+    mat <- cor2pcor(cov)
+    mat <- matrix(mat, ncol=ncol(lr), nrow=ncol(lr))
+    class(mat) <- "matrix"
+    lambda <- attr(cov, "lambda")
+
+    # compute pcor
+    pr <- suppressWarnings(propr(X, metric = "pcor.shrink", ivar=d1[[i]]))
+
+    # expect counts to have zeros replaced
+    expect_true(
+      all(pr@counts == ct)
+    )
+
+    # expect that logratio are equal
+    expect_true(
+      all(pr@logratio == lr)
+    )
+
+    # expect computed coefficients are equal
+    expect_true(
+      all(round(pr@matrix, 8) == round(mat, 8), na.rm=T)
+    )
+
+    # check shrinkage is not applied
+    expect_equal(
+      pr@lambda,
+      lambda
+    )
+  }
+})
+
+
+test_that("pcor.shrink is correct when ivar is NA using previously transformed data", {
+
+  d1 <- list("clr", 5, c(1,3))
+  d2 <- list(c(1:5), 5, c(1,3))
+
+  for (i in 1:length(d1)){
+
+    # compute pcor manually
+    ct <- simple_zero_replacement(X)
+    lr <- logratio_without_alpha(ct, d2[[i]])
+    cov <- suppressWarnings(cov.shrink(lr))
+    mat <- cor2pcor(cov)
+    mat <- matrix(mat, ncol=ncol(lr), nrow=ncol(lr))
+    class(mat) <- "matrix"
+    lambda <- attr(cov, "lambda")
+
+    # compute pcor
+    pr <- suppressWarnings(propr(lr, metric = "pcor.shrink", ivar=NA))
+
+    # expect counts to contain the previously transformed data
+    expect_true(
+      all(pr@counts == lr)
+    )
+
+    # expect that the logratio slot also contains the exact input data
+    expect_true(
+      all(pr@logratio == lr)
+    )
+
+    # expect computed coefficients are equal
+    expect_true(
+      all(round(pr@matrix, 8) == round(mat, 8), na.rm=T)
+    )
+
+    # check shrinkage is not applied
+    expect_equal(
+      pr@lambda,
+      lambda
+    )
+  }
+})
+
+test_that("pcor.shrink with alr and clr are not the same", {
+
+    # compute pcor.shrink with alr
+    pr <- suppressWarnings(propr(X, metric = "pcor.shrink", ivar=5))
+    pcor_alr <- getMatrix(pr)[1:4, 1:4]
+
+    # compute pcor.shrink with clr
+    pr <- propr(X, metric = "pcor.shrink", ivar="clr")
+    pcor_clr <- getMatrix(pr)[1:4, 1:4]
+
+    expect_false(
+      all(round(pcor_alr, 8) == round(pcor_clr, 8))
+    )
+
+})
diff --git a/tests/testthat/test-PROPR-phi-rho-phs-cor.R b/tests/testthat/test-PROPR-phi-rho-phs-cor.R
index 96c9490..c18065c 100644
--- a/tests/testthat/test-PROPR-phi-rho-phs-cor.R
+++ b/tests/testthat/test-PROPR-phi-rho-phs-cor.R
@@ -53,4 +53,5 @@ test_that("calculating rho from phi matches rho", {
     cor_,
     cor
   )
+
 })
diff --git a/tests/testthat/test-PROPR-updatepermutes.R b/tests/testthat/test-PROPR-updatepermutes.R
deleted file mode 100644
index 88b956f..0000000
--- a/tests/testthat/test-PROPR-updatepermutes.R
+++ /dev/null
@@ -1,52 +0,0 @@
-library(testthat)
-library(propr)
-
-N <- 100
-a <- seq(from = 5, to = 15, length.out = N)
-b <- a * rnorm(N, mean = 1, sd = 0.1)
-c <- rnorm(N, mean = 10)
-d <- rnorm(N, mean = 10)
-e <- rep(10, N)
-X <- data.frame(a, b, c, d, e)
-
-# compute pcor.bshrink without fixed seed
-pcorbshrink1 <- propr(X, metric = "pcor.bshrink", p=10)
-pcorbshrink1 <- updateCutoffs(pcorbshrink1)
-
-pcorbshrink2 <- propr(X, metric = "pcor.bshrink", p=10)
-pcorbshrink2 <- updateCutoffs(pcorbshrink2)
-
-# compute pcor.bshrink with fixed seed
-pcorbshrink1_ <- propr(X, metric = "pcor.bshrink", p=10, fixseed=TRUE)
-pcorbshrink1_ <- updateCutoffs(pcorbshrink1_)
-
-pcorbshrink2_ <- propr(X, metric = "pcor.bshrink", p=10, fixseed=TRUE)
-pcorbshrink2_ <- updateCutoffs(pcorbshrink2_)
-
-# test that the results are as expected
-test_that("test that fdr will stay the same only if fixseed=TRUE", {
-
-  expect_false(
-    isTRUE(all.equal(
-        pcorbshrink1@permutes,
-        pcorbshrink2@permutes
-    ))
-  )
-
-  expect_false(
-    isTRUE(all.equal(
-        pcorbshrink1@fdr,
-        pcorbshrink2@fdr
-    ))
-  )
-
-  expect_equal(
-    pcorbshrink1_@permutes,
-    pcorbshrink2_@permutes
-  )
-
-  expect_equal(
-    pcorbshrink1_@fdr,
-    pcorbshrink2_@fdr
-  )
-})
\ No newline at end of file
diff --git a/tests/testthat/test-RCPP-count.R b/tests/testthat/test-RCPP-count.R
new file mode 100644
index 0000000..9a18fa5
--- /dev/null
+++ b/tests/testthat/test-RCPP-count.R
@@ -0,0 +1,18 @@
+library(testthat)
+library(propr)
+
+test_that("count_greater_than and count_less_than work properly", {
+
+    # define values
+    values <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
+
+    expect_equal(propr:::count_greater_than(values, 5), 5)
+    expect_equal(propr:::count_greater_than(values, 8), 2)
+    expect_equal(propr:::count_less_than(values, 5),4)
+    expect_equal(propr:::count_less_than(values, 8), 7)
+
+    expect_equal(propr:::count_greater_equal_than(values, 5), 6)
+    expect_equal(propr:::count_greater_equal_than(values, 8), 3)
+    expect_equal(propr:::count_less_equal_than(values, 5), 5)
+    expect_equal(propr:::count_less_equal_than(values, 8), 8)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-SHARED-getAdjacency-propd.R b/tests/testthat/test-SHARED-getAdjacency-propd.R
new file mode 100644
index 0000000..9b0ce51
--- /dev/null
+++ b/tests/testthat/test-SHARED-getAdjacency-propd.R
@@ -0,0 +1,90 @@
+library(testthat)
+library(propr)
+library(MASS)
+
+# define data
+data(crabs)
+x <- crabs[,4:8]  # data matrix with 5 variables
+y <- crabs[,1]    # group vector
+
+test_that("getAdjacencyFDR works properly for theta", {
+
+    # get propd object
+    pr <- propd(x, as.character(y), p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(x), ncol = ncol(x))
+    adj_expected[propr:::getMatrix(pr) <= getCutoffFDR(pr)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(x)
+    colnames(adj_expected) <- colnames(x)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR and getSignificantResultsFDR return coherent results", {
+
+    # get propd object
+    pr <- propd(x, as.character(y), p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get significant results
+    results <- getSignificantResultsFDR(pr)
+
+    # check that the values are correct
+    for (i in 1:nrow(results)){
+        expect_equal(adj[results[i,1], results[i,2]], 1)
+    }
+})
+
+test_that("getAdjacencyFstat works properly", {
+
+    for (fdr_adjusted in c(TRUE, FALSE)){
+
+        # get propd object
+        pr <- propd(x, as.character(y), p=10)
+        pr <- updateF(pr)
+
+        # get adjacency matrix
+        adj <- getAdjacencyFstat(pr, fdr_adjusted=fdr_adjusted)
+
+        # get expected adjacency matrix
+        adj_expected <- matrix(0, nrow = ncol(x), ncol = ncol(x))
+        adj_expected[propr:::getMatrix(pr) <= getCutoffFstat(pr, fdr_adjusted=fdr_adjusted)] <- 1
+        adj_expected[diag(adj_expected)] <- 1
+        rownames(adj_expected) <- colnames(x)
+        colnames(adj_expected) <- colnames(x)
+
+        # check that the values are correct
+        expect_equal(adj, adj_expected)
+    }
+})
+
+test_that("getAdjacencyFstat and getSignificantResultsFstat return coherent results", {
+
+    for (fdr in c(TRUE, FALSE)){
+
+        # get propd object
+        pr <- propd(x, as.character(y), p=10)
+        pr <- updateF(pr)
+
+        # get adjacency matrix
+        adj <- getAdjacencyFstat(pr, fdr=fdr)
+
+        # get significant results
+        results <- getSignificantResultsFstat(pr, fdr=fdr)
+
+        # check that the values are correct
+        for (i in 1:nrow(results)){
+            expect_equal(adj[results[i,1], results[i,2]], 1)
+        }
+    }
+})
diff --git a/tests/testthat/test-SHARED-getAdjacency-propr.R b/tests/testthat/test-SHARED-getAdjacency-propr.R
new file mode 100644
index 0000000..4bd2d7d
--- /dev/null
+++ b/tests/testthat/test-SHARED-getAdjacency-propr.R
@@ -0,0 +1,155 @@
+library(testthat)
+library(propr)
+
+# define data matrix
+set.seed(123)
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+
+
+test_that("getAdjacencyFDR returns the expected values for pcor.bshrink - clr", {
+
+    # get propr object
+    pr <- propr(X, metric = "pcor.bshrink", ivar='clr', p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix <= getCutoffFDR(pr, fdr=0.025, positive = FALSE) | pr@matrix >= getCutoffFDR(pr, fdr=0.025, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR returns the expected values for pcor.bshrink - alr", {
+
+    # get propr object
+    pr <- propr(X, metric = "pcor.bshrink", ivar='alr', p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix <= getCutoffFDR(pr, fdr=0.025, positive = FALSE) | pr@matrix >= getCutoffFDR(pr, fdr=0.025, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR returns the expected values for rho - clr", {
+
+    # get propr object
+    pr <- propr(X, metric = "rho", ivar='clr', p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix >= getCutoffFDR(pr, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR returns the expected values for rho - 5", {
+
+    # get propr object
+    pr <- propr(X, metric = "rho", ivar=5, p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix >= getCutoffFDR(pr, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR returns the expected values for phs - clr", {
+
+    # get propr object
+    pr <- propr(X, metric = "phs", ivar='clr', p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix <= getCutoffFDR(pr, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR returns the expected values for phs - 5", {
+
+    # get propr object
+    pr <- propr(X, metric = "phs", ivar=5, p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+    # get adjacency matrix
+    adj <- getAdjacencyFDR(pr)
+
+    # get expected adjacency matrix
+    adj_expected <- matrix(0, nrow = ncol(X), ncol = ncol(X))
+    adj_expected[pr@matrix <= getCutoffFDR(pr, positive = TRUE)] <- 1
+    adj_expected[diag(adj_expected)] <- 1
+    rownames(adj_expected) <- colnames(X)
+    colnames(adj_expected) <- colnames(X)
+
+    # check that the values are correct
+    expect_equal(adj, adj_expected)
+})
+
+test_that("getAdjacencyFDR and getSignificantResultsFDR return coherent results", {
+    
+    for (metric in c('rho', 'phi', 'phs', 'pcor', 'pcor.shrink', 'pcor.bshrink')){
+        print(metric)
+
+        # get propr object
+        pr <- propr(X, metric=metric, p=10)
+        pr <- updateCutoffs(pr, number_of_cutoffs=100)
+
+        # get adjacency matrix
+        adj <- getAdjacencyFDR(pr)
+
+        # get significant results
+        results <- getSignificantResultsFDR(pr)
+
+        # check that the values are correct
+        for (i in 1:nrow(results)){
+            expect_equal(adj[results[i,1], results[i,2]], 1)
+        }
+    }
+})
diff --git a/tests/testthat/test-SHARED-getCutoff-propd.R b/tests/testthat/test-SHARED-getCutoff-propd.R
new file mode 100644
index 0000000..d8ee367
--- /dev/null
+++ b/tests/testthat/test-SHARED-getCutoff-propd.R
@@ -0,0 +1,56 @@
+library(testthat)
+library(propr)
+library(MASS)
+
+# define data
+data(crabs)
+x <- crabs[,4:8]  # data matrix with 5 variables
+y <- crabs[,1]    # group vector
+
+test_that("getCutoffFDR gets the correct cutoff", {
+
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr <- propd(x, as.character(y), p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+    # get cutoff
+    cutoff <- getCutoffFDR(pr, fdr = 0.05, window_size = 1)
+
+    # check that the cutoff is correct
+    expect_equal(round(cutoff, 4), 0.9738)  # for the moment the expected value is manually calculated
+})
+
+test_that("getCutoffFstat gets the correct cutoff when FDR correction is considered", {
+
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr <- propd(x, as.character(y), p=10)
+    pr <- updateF(pr)
+
+    # get cutoff 
+    cutoff_expected <- max(pr@results$theta[pr@results$FDR <= 0.05])
+    cutoff_actual <- getCutoffFstat(pr, pval = 0.05, fdr = TRUE)
+
+    # check that the cutoff is correct
+    expect_equal(round(cutoff_actual, 6), round(cutoff_expected, 6))
+})
+
+test_that("getCutoffFstat gets the correct cutoff when FDR correction is not considered", {
+
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr <- propd(x, as.character(y), p=10)
+    pr <- updateF(pr)
+
+    # get cutoff 
+    pval <- 0.05
+    K <- length(unique(pr@group))
+    N <- length(pr@group) + pr@dfz # population-level metric (i.e., N)
+    Q <- stats::qf(pval, K - 1, N - K, lower.tail = FALSE)
+    cutoff_expected <- (N - 2) / (Q + (N - 2))
+    cutoff_actual <- getCutoffFstat(pr, fdr = FALSE)
+
+    # check that the cutoff is correct
+    expect_equal(round(cutoff_actual, 6), round(cutoff_expected, 6))
+})
\ No newline at end of file
diff --git a/tests/testthat/test-SHARED-getCutoff-propr.R b/tests/testthat/test-SHARED-getCutoff-propr.R
new file mode 100644
index 0000000..03a348c
--- /dev/null
+++ b/tests/testthat/test-SHARED-getCutoff-propr.R
@@ -0,0 +1,45 @@
+library(testthat)
+library(propr)
+
+# define data matrix
+set.seed(123)
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+
+test_that("test that getCutoff gets the correct cutoff", {
+
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr <- propr(X, metric = "pcor.bshrink", p=10)
+    pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+    # check cutoff is correct for positive values
+    cutoff <- getCutoffFDR(pr, fdr = 0.05, window_size = 1, positive = TRUE)
+    expect_equal(round(cutoff, 4), 0.6573)  # for the moment the expected value is manually calculated
+
+    # check cutoff is correct for negative values
+    cutoff <- getCutoffFDR(pr, fdr = 0.05, window_size = 1, positive = FALSE)
+    expect_equal(round(cutoff, 4), -0.0313)  # for the moment the expected value is manually calculated
+})
+
+test_that("test that moving average works as expected", {
+
+    x <- c(NA, 1, 1, 2, NA, 2, 3, 3, 3, 4, 4, 5)
+
+    # define expected moving average
+    y1 <- x 
+    y2 <- c(NA, 1, 1.5, 2, NA, 2.5, 3, 3, 3.5, 4, 4.5, 5)
+    y3 <- c(NA, 1, 1.33, 1.5, NA, 2.5, 2.67, 3, 3.33, 3.67, 4.33, 4.5)
+    y4 <- c(NA, 1.33, 1.33, 1.67, NA, 2.67, 2.75, 3.25, 3.5, 4, 4.33, 4.5)
+
+    # check
+    expect_equal(round(propr:::getMovingAverage(x, 1),2), round(y1,2))
+    expect_equal(round(propr:::getMovingAverage(x, 2),2), round(y2,2))
+    expect_equal(round(propr:::getMovingAverage(x, 3),2), round(y3,2))
+    expect_equal(round(propr:::getMovingAverage(x, 4),2), round(y4,2))
+})
diff --git a/tests/testthat/test-SHARED-getResults-propd.R b/tests/testthat/test-SHARED-getResults-propd.R
new file mode 100644
index 0000000..695f66f
--- /dev/null
+++ b/tests/testthat/test-SHARED-getResults-propd.R
@@ -0,0 +1,60 @@
+library(testthat)
+library(propr)
+library(MASS)
+
+# define data
+data(crabs)
+x <- crabs[,4:8]  # data matrix with 5 variables
+y <- crabs[,1]    # group vector
+
+test_that("getResults works as expected", {
+
+  # get propd object
+  pr <- propd(x, as.character(y), p=10)
+  
+  # get results
+  results <- getResults(pr)
+  
+  # check that the values are correct
+  expect_equal(pr@results[,-c(1:2)], results[,-c(1:2)])
+
+  # check that the variable names are corretly replaced
+  expect_equal(results[,1], colnames(x)[pr@results[,1]])
+  expect_equal(results[,2], colnames(x)[pr@results[,2]])
+
+  # check that the order of pairs are as expected
+  ord <- list(c(2,1), c(3,1), c(3,2), c(4,1), c(4,2), c(4,3), c(5,1), c(5,2), c(5,3), c(5,4))
+  for (i in 1:10) {
+    expect_equal(results[i,1], colnames(x)[ord[[i]][1]])
+    expect_equal(results[i,2], colnames(x)[ord[[i]][2]])
+  }
+})
+
+test_that("getSignificantResultsFDR works as expected", {
+
+  # get propd object
+  pr <- propd(x, as.character(y), p=10)
+  pr <- updateCutoffs(pr, number_of_cutoffs=10)
+  
+  # get expected results
+  cutoff <- getCutoffFDR(pr, fdr = 0.05, window_size = 1)
+  expected <- pr@results[which(pr@results$theta <= cutoff),]
+  
+  # get significant results
+  results <- getSignificantResultsFDR(pr, fdr = 0.05)
+  
+  # check that the values are correct
+  expect_equal(results$theta, expected$theta)
+})
+
+test_that("getSignificantResultsFstat works as expected", {
+
+  # get propd object
+  pr <- propd(x, as.character(y), p=10)
+  pr <- updateF(pr)
+
+  # expect that the Fstat values are smaller or equal than the cutoff
+  expect_true(all(getSignificantResultsFstat(pr, fdr=F)$theta <= getCutoffFstat(pr, fdr=F)))
+  expect_true(all(getSignificantResultsFstat(pr, fdr=T)$fdr <= 0.05))
+
+})
diff --git a/tests/testthat/test-SHARED-getResults-propr.R b/tests/testthat/test-SHARED-getResults-propr.R
new file mode 100644
index 0000000..f4b5992
--- /dev/null
+++ b/tests/testthat/test-SHARED-getResults-propr.R
@@ -0,0 +1,55 @@
+library(testthat)
+library(propr)
+
+# define data matrix
+set.seed(123)
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+
+test_that("test that getResults works as expected", {
+  
+  # get propr object
+  pr <- propr(X, metric = "pcor.bshrink", p=10)
+  
+  # get results
+  results <- getResults(pr)
+  
+  # check that the values are correct
+  expect_equal(pr@results[,c(3:7)], results[,c(3:7)])
+
+  # check that the variable names are corretly replaced
+  expect_equal(results[,1], colnames(X)[pr@results[,1]])
+  expect_equal(results[,2], colnames(X)[pr@results[,2]])
+
+  # check that the order of pairs are as expected
+  ord <- list(c(2,1), c(3,1), c(3,2), c(4,1), c(4,2), c(4,3), c(5,1), c(5,2), c(5,3), c(5,4))
+  for (i in 1:10) {
+    expect_equal(results[i,1], colnames(X)[ord[[i]][1]])
+    expect_equal(results[i,2], colnames(X)[ord[[i]][2]])
+  }
+})
+
+test_that("test that getSignificantResultsFDR works as expected", {
+  
+  # get propr object
+  pr <- propr(X, metric = "pcor.bshrink", p=10)
+  pr <- updateCutoffs(pr, number_of_cutoffs=10)
+  
+  # get expected results
+  cutoff_positive <- getCutoffFDR(pr, fdr = 0.05, window_size = 1, positive = TRUE)
+  cutoff_negative <- getCutoffFDR(pr, fdr = 0.05, window_size = 1, positive = FALSE)
+  expected_positive <- pr@results$propr[which(pr@results$propr >= cutoff_positive)]
+  expected_negative <- pr@results$propr[which(pr@results$propr <= cutoff_negative)]
+
+  # get significant results
+  results <- getSignificantResultsFDR(pr, fdr = 0.05)
+
+  # check that the values are correct
+  expect_equal(results$propr[results$propr>=0], expected_positive)
+  expect_equal(results$propr[results$propr<0], expected_negative)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-SHARED-updateCutoffs-propd.R b/tests/testthat/test-SHARED-updateCutoffs-propd.R
new file mode 100644
index 0000000..cc338d1
--- /dev/null
+++ b/tests/testthat/test-SHARED-updateCutoffs-propd.R
@@ -0,0 +1,103 @@
+library(testthat)
+library(propr)
+library(MASS)
+
+# define data
+data(crabs)
+x <- crabs[,4:8]  # data matrix with 5 variables
+y <- crabs[,1]    # group vector
+
+
+test_that("updateCutoffs.propd properly set up cutoffs", {
+
+  # get propd object and update cutoffs
+  pd <- propd(x, as.character(y), p=10)
+  pd <- updateCutoffs(pd, number_of_cutoffs=10)
+
+  # get cutoffs
+  cutoffs <- as.numeric( quantile(pd@results$theta, probs = seq(0, 1, length.out = 10)) )
+
+  # check that cutoffs are properly defined
+  expect_equal(pd@fdr$cutoff, cutoffs)
+})
+
+test_that("updateCutoffs.propd properly calculates truecounts", {
+
+  # get propd object and update cutoffs
+  pd <- propd(x, as.character(y), p=10)
+  pd <- updateCutoffs(pd, number_of_cutoffs=10)
+
+  # get truecounts
+  truecounts <- sapply(
+    pd@fdr$cutoff, 
+    function(cut) sum(pd@results$theta < cut)
+  )
+  truecounts_manual <- c(0:9)
+
+  # check that truecounts are properly defined
+  expect_equal(pd@fdr$truecounts, truecounts)
+  expect_equal(pd@fdr$truecounts, truecounts_manual)
+})
+
+test_that("updateCutoffs.propd properly calculates randcounts", {
+
+  # get propd object and update cutoffs
+  set.seed(0)
+  pd <- propd(x, as.character(y), p=10)
+  pd <- updateCutoffs(pd, number_of_cutoffs=10)
+
+  # get permuted values
+  randcounts <- rep(0, 10)
+  for (k in 1:10){
+    shuffle <- pd@permutes[,k]
+    ct.k <- pd@counts[shuffle,]
+    pd.k <- suppressMessages(propd(
+      ct.k,
+      group = pd@group,
+      alpha = pd@alpha,
+      p = 0,
+      weighted = pd@weighted
+    ))
+    pkt <- pd.k@results$theta
+    for (cut in 1:length(pd@fdr$cutoff)){
+        randcounts[cut] <- randcounts[cut] + sum(pkt < pd@fdr$cutoff[cut])
+    }
+  }
+  randcounts <- randcounts / 10
+  randcounts_manual <- c(0, 0, 0, 0, 0, 0, 0, 0.1, 0.1, 5.4)
+
+  # check that the permutation tests work properly
+  expect_equal(pd@fdr$randcounts, randcounts)
+  expect_equal(round(pd@fdr$randcounts, 1), randcounts_manual)
+})
+
+test_that("updateCutoffs.propd is reproducible when seed is set", {
+  
+    # get propd object and update cutoffs
+    set.seed(0)
+    pr1 <- propd(x, as.character(y), p=10)
+    pr1 <- updateCutoffs(pr1, number_of_cutoffs=10)
+    set.seed(0)
+    pr2 <- propd(x, as.character(y), p=10)
+    pr2 <- updateCutoffs(pr2, number_of_cutoffs=10)
+    pr3 <- propd(x, as.character(y), p=10)
+    pr3 <- updateCutoffs(pr3, number_of_cutoffs=10)
+  
+    # check that fdr are the same only when seed is the same
+    expect_equal(pr1@fdr, pr2@fdr)
+    expect_false(isTRUE(all.equal(pr1@fdr, pr3@fdr)))
+})
+
+test_that("updateCutoffs.propd works when ncores > 1", {
+  
+    # get propd object and update cutoffs
+    set.seed(0)
+    pr1 <- propd(x, as.character(y), p=10)
+    pr1 <- updateCutoffs(pr1, number_of_cutoffs=10, ncores=1)
+    set.seed(0)
+    pr2 <- propd(x, as.character(y), p=10)
+    pr2 <- updateCutoffs(pr2, number_of_cutoffs=10, ncores=2)
+
+    # check that fdr are the same
+    expect_equal(pr1@fdr, pr2@fdr)
+})
diff --git a/tests/testthat/test-SHARED-updateCutoffs-propr.R b/tests/testthat/test-SHARED-updateCutoffs-propr.R
new file mode 100644
index 0000000..f6155b1
--- /dev/null
+++ b/tests/testthat/test-SHARED-updateCutoffs-propr.R
@@ -0,0 +1,114 @@
+library(testthat)
+library(propr)
+
+# define data matrix
+set.seed(123)
+N <- 100
+a <- seq(from = 5, to = 15, length.out = N)
+b <- a * rnorm(N, mean = 1, sd = 0.1)
+c <- rnorm(N, mean = 10)
+d <- rnorm(N, mean = 10)
+e <- rep(10, N)
+X <- data.frame(a, b, c, d, e)
+
+test_that("updateCutoffs.propr properly set up cutoffs", {
+
+  # get propr object and update cutoffs
+  pr <- propr(X, metric = "pcor.bshrink", p=10)
+  pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+  # get cutoffs
+  cutoffs <- as.numeric( quantile(pr@matrix[lower.tri(pr@matrix)], probs = seq(0, 1, length.out = 10)) )
+
+  # check that cutoffs are properly defined
+  expect_equal(pr@fdr$cutoff, cutoffs)
+})
+
+test_that("updateCutoffs.propr properly calculates truecounts", {
+
+  # get propr object and update cutoffs
+  pr <- propr(X, metric = "pcor.bshrink", p=10)
+  pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+  # get truecounts
+  truecounts1 <- sapply(
+    pr@fdr$cutoff[pr@fdr$cutoff < 0], 
+    function(cut) sum(pr@results$propr < cut)
+  )
+  truecounts2 <- sapply(
+    pr@fdr$cutoff[pr@fdr$cutoff >= 0], 
+    function(cut) sum(pr@results$propr > cut)
+  )
+  truecounts <- c(truecounts1, truecounts2)
+  truecounts_manual <- c(0,1,2, 6,5,4,3,3,1,0)
+
+  # check that truecounts are properly defined
+  expect_equal(pr@fdr$truecounts, truecounts)
+  expect_equal(pr@fdr$truecounts, truecounts_manual)
+})
+
+test_that("updateCutoffs.propr properly calculates randcounts", {
+
+  # get propr object and update cutoffs
+  set.seed(0)
+  pr <- propr(X, metric = "pcor.bshrink", p=10)
+  pr <- updateCutoffs(pr, number_of_cutoffs=10)
+
+  # get permuted values
+  randcounts <- rep(0, 10)
+  for (k in 1:10){
+    ct.k <- pr@permutes[[k]]
+    pr.k <- suppressMessages(propr(
+      ct.k,
+      pr@metric,
+      ivar = pr@ivar,
+      alpha = pr@alpha,
+      p = 0
+    ))
+    pkt <- pr.k@results$propr
+    for (cut in 1:length(pr@fdr$cutoff)){
+      cutoff <- pr@fdr$cutoff[cut]
+      if (cutoff >= 0) {
+        randcounts[cut] <- randcounts[cut] + sum(pkt > cutoff)
+      } else {
+        randcounts[cut] <- randcounts[cut] + sum(pkt < cutoff)
+      }
+    }
+  }
+  randcounts <- randcounts / 10
+  randcounts_manual <- c(0,0,0, 9.7,9.7,9.2,8.2,0,0,0)
+
+  # check that the permutation tests work properly
+  expect_equal(pr@fdr$randcounts, randcounts)
+})
+
+test_that("updateCutoffs.propr is reproducible when seed is set", {
+  
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr1 <- propr(X, metric = "pcor.bshrink", p=10)
+    pr1 <- updateCutoffs(pr1, number_of_cutoffs=10)
+    set.seed(0)
+    pr2 <- propr(X, metric = "pcor.bshrink", p=10)
+    pr2 <- updateCutoffs(pr2, number_of_cutoffs=10)
+    pr3 <- propr(X, metric = "pcor.bshrink", p=10)
+    pr3 <- updateCutoffs(pr3, number_of_cutoffs=10)
+  
+    # check that fdr are the same only when seed is the same
+    expect_equal(pr1@fdr, pr2@fdr)
+    expect_false(isTRUE(all.equal(pr1@fdr, pr3@fdr)))
+})
+
+test_that("updateCutoffs.propr works when ncores > 1", {
+
+    # get propr object and update cutoffs
+    set.seed(0)
+    pr1 <- propr(X, metric = "pcor.bshrink", p=10)
+    pr1 <- updateCutoffs(pr1, number_of_cutoffs=10, ncores=1)
+    set.seed(0)
+    pr2 <- propr(X, metric = "pcor.bshrink", p=10)
+    pr2 <- updateCutoffs(pr2, number_of_cutoffs=10, ncores=2)
+
+    # check that fdr are the same
+    expect_equal(pr1@fdr, pr2@fdr)
+})
diff --git a/tests/testthat/test-SHARED-updatepermutes.R b/tests/testthat/test-SHARED-updatepermutes.R
new file mode 100644
index 0000000..092c5b4
--- /dev/null
+++ b/tests/testthat/test-SHARED-updatepermutes.R
@@ -0,0 +1,64 @@
+library(testthat)
+library(propr)
+
+
+test_that("propr - test that permute stay the same only when seed is given",{
+
+  # define data matrix
+  N <- 100
+  a <- seq(from = 5, to = 15, length.out = N)
+  b <- a * rnorm(N, mean = 1, sd = 0.1)
+  c <- rnorm(N, mean = 10)
+  d <- rnorm(N, mean = 10)
+  e <- rep(10, N)
+  X <- data.frame(a, b, c, d, e)
+
+  # test when seed is given
+  set.seed(0)
+  pr1 <- propr(X, metric = "pcor.bshrink", p=10)
+  set.seed(0)
+  pr2 <- propr(X, metric = "pcor.bshrink", p=10)
+  expect_equal(
+    pr1@permutes,
+    pr2@permutes
+  )
+
+  # test when seed is not given
+  pr1 <- propr(X, metric = "pcor.bshrink", p=10)
+  pr2 <- propr(X, metric = "pcor.bshrink", p=10)
+  expect_false(
+    isTRUE(all.equal(
+        pr1@permutes,
+        pr2@permutes
+    ))
+  )
+})
+
+test_that("propd - test that permute stay the same only when seed is given",{
+
+  # define data
+  x <- iris[,1:4]  # data matrix with 4 variables
+  y <- iris[,5]    # group vector
+  v <- vegan::rda(log(x[,1]/x[,2]) ~ y)
+
+  # test when seed is given
+  set.seed(0)
+  pd1 <- propd(x, as.character(y), p=10)
+  set.seed(0)
+  pd2 <- propd(x, as.character(y), p=10)
+  expect_equal(
+    pd1@permutes,
+    pd2@permutes
+  )
+
+  # test when seed is not given
+  pd1 <- propd(x, as.character(y), p=10)
+  pd2 <- propd(x, as.character(y), p=10)
+  expect_false(
+    isTRUE(all.equal(
+        pd1@permutes,
+        pd2@permutes
+    ))
+  )
+})
+
diff --git a/tests/testthat/test-SHARED-zero.R b/tests/testthat/test-SHARED-zero.R
new file mode 100644
index 0000000..4303b54
--- /dev/null
+++ b/tests/testthat/test-SHARED-zero.R
@@ -0,0 +1,18 @@
+library(testthat)
+library(propr)
+
+test_that("simple_zero_replacement replaces zeros correctly", {
+
+  # Example data with zeros
+  data_with_zeros <- matrix(c(1, 2, 3, 0, 5, 6, 0, 8, 9), nrow = 3, byrow = TRUE)
+
+  # Apply the simple_zero_replacement function to the example data
+  replaced_data <- simple_zero_replacement(data_with_zeros)
+
+  # Define the expected output after zero replacement
+  expected_output <- matrix(c(1, 2, 3, 1, 5, 6, 1, 8, 9), nrow = 3, byrow = TRUE)
+
+  # Compare the output with the expected output
+  expect_identical(replaced_data, expected_output)
+
+})
\ No newline at end of file