From 66e3e7e4d5b682d2f7a1e88b8e88ab0a93c7d489 Mon Sep 17 00:00:00 2001 From: pribiller Date: Mon, 17 Aug 2026 14:29:37 +0900 Subject: [PATCH 1/4] Add new estimate for the expected number of inversions Add new estimate for the expected number of inversions based on Theorem 3 from Berestycki and Durrett (2006). --- R/inversionEstimate_BD.R | 342 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 342 insertions(+) create mode 100644 R/inversionEstimate_BD.R diff --git a/R/inversionEstimate_BD.R b/R/inversionEstimate_BD.R new file mode 100644 index 00000000..3a07525e --- /dev/null +++ b/R/inversionEstimate_BD.R @@ -0,0 +1,342 @@ + +#' Computes the expected number of cycles in a random graph +#' after k edges. +#' +#' Implements Theorem 3 from Berestycki and Durrett (2006), +#' also described in Equation 6 in Biller et al. (2015). +#' +#' @param k An integer: The number of edges in the random graph, +#' equivalent to the number of inversions. +#' +#' @param N An integer: The number of vertices in the random graph, +#' equivalent to the number of markers+1. +#' +#' @references Berestycki, Nathanaël, and Rick Durrett. "A phase transition in the random transposition random walk." Discrete Mathematics and Theoretical Computer Science. Discrete Mathematics and Theoretical Computer Science, 2003. +#' @references Hannenhalli, Sridhar, and Pavel A. Pevzner. "Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals." Journal of the ACM (JACM) 46.1 (1999): 1-27. +#' @references Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evolution by double cut-and-join." BMC bioinformatics 16.Suppl 14 (2015): S7. +#' +#' @author Priscila Biller +nb_cycles_k <- function(k, N){ + + # Invalid value (negative number of inversions or not enough markers). + if ((N < 2) | (k < 0)){ + NA + + # No inversions. Genomes are the same. + } else if(k == 0) { + N + + # At least one inversion. + } else { + + nb_cyc_k <- 0 + # This loop is supposed to go until infinity (sum until infinity), + # but as i grows, the contribution of the i-th term gets smaller and smaller. + # Here, the computation is capped at 10 for efficiency. + # It does not seem to affect much the accuracy. + for(i in 1:10) { + term_base <- 2.0*k/N + term_k <- term_base*exp(-term_base) + nb_cyc_k <- nb_cyc_k + ((1.0 / term_base) * ( i**(i-2) / factorial(i) )* term_k**i) + } + nb_cyc_k*N + } +} + +#' The setup step computes the expected number of cycles +#' after k edges for various values of k. +#' +#' The expected number of cycles can take any value between 1 +#' and n (the number of markers). For each of these values, +#' the method computes at least one k. +#' +#' Notice that this is the most time-consuming step +#' in the estimate. Since it depends only on the number +#' n of markers, it is recommended that it be computed once +#' for all permutations of n markers, saving time. +#' +#' @param n An integer: The number of markers. +#' A marker represents an aligned block, a gene, +#' or any other genomic region being rearranged. +#' +#' @author Priscila Biller +inversionEstimate_BD_setup <- function(n){ + + # The components of the breakpoint graph for random reversals + # on n-1 markers are related to the cycles of random + # transposition on n markers. + N <- n+1 + + cyc_all <- list(list(k=0,nb_cyc_k=N)) + k <- 1 + pos_k <- 2 + + while(cyc_all[[length(cyc_all)]]$nb_cyc_k > 1) { + + # Computes the expected number of cycles in a random graph after k edges. + cyc_all[[pos_k]] <- list(k=k,nb_cyc_k=nb_cycles_k(k, N)) + + # As the number of events increases, the expected number of cycles + # requires more events to be reduced by one. To avoid computing many values of k + # that yield approximately the same expected number of cycles, the method skips some ks: + # it jumps when the difference in the number of cycles between the most recently + # computed value and its predecessor is less than one. + k_diff <- cyc_all[[pos_k]]$k - cyc_all[[pos_k-1]]$k + # The number of cycles decreases as the number of events increases. + nb_cyc_diff <- as.integer(cyc_all[[pos_k-1]]$nb_cyc_k) - as.integer(cyc_all[[pos_k]]$nb_cyc_k) + + if((nb_cyc_diff < 2) || (k_diff < 2)){ + last_pos <- length(cyc_all) + k_diff <- cyc_all[[last_pos]]$k - cyc_all[[last_pos-1]]$k + nb_cyc_diff <- as.integer(cyc_all[[last_pos]]$nb_cyc_k) - as.integer(cyc_all[[last_pos-1]]$nb_cyc_k) + if(nb_cyc_diff == 0){ + k_diff <- k_diff*2 + } + k <- cyc_all[[last_pos]]$k + k_diff + pos_k <- last_pos + 1 + } else { + # Append a new value in the current position. + k <- cyc_all[[pos_k-1]]$k + as.integer(k_diff)/2 + cyc_all <- c(cyc_all[1:(pos_k-1)], NA, cyc_all[pos_k:length(cyc_all)]) + } + } + cyc_all +} + + +#' It finds the number of inversions whose expected number of cycles +#' is closest to the observed number of cycles. +#' +#' This function is an alternative to the function "inversionEstimate_BD_many", +#' which uses the function "inversionEstimate_BD_setup". +#' +#' This function should be preferred when there are few inversion estimates to be computed. +#' +#' @param n An integer: The number of markers. +#' A marker represents an aligned block, a gene, +#' or any other genomic region being rearranged. +#' +#' @param obs_nb_cycles An integer: The observed number of cycles +#' in the breakpoint graph of two genomes. +#' +#' @author Priscila Biller +inversionEstimate_BD_single <- function(n, obs_nb_cycles){ + + # The components of the breakpoint graph for random reversals + # on n-1 markers are related to the cycles of random + # transposition on n markers. + N <- n+1 + + # At the start of the evolutionary process, when the two genomes are identical, + # the number of cycles is equal to the number of markers, the maximum possible value. + # Knowing that an inversion can split one cycle into two or merge two cycles into one, + # if c cycles are observed, then at least N-c inversions must have occurred. + k_beg <- 0 + k_end <- N + + # These variables hold the previous k_beg and k_end (to avoid infinite loops): + # If this case occurs, there is a bug. + k_beg_prev <- k_beg-1 + k_end_prev <- k_end+1 + + # These variables hold the last valid k_beg and k_end: + # ------[k_beg ------ point ------ k_end]------ + k_beg_valid <- 0 + k_end_valid <- 100000*k_end + + nb_cyc_k_beg <- N + nb_cyc_k_end <- nb_cycles_k(k_end, N) + + # Identical permutations. No inversions occurred. + if(obs_nb_cycles == nb_cyc_k_beg){ + return(0) + } + + #cat(paste("START: [", k_beg, " (", nb_cyc_k_beg, "), ", k_end, " (", nb_cyc_k_end, ")] / ", obs_nb_cycles, "\n")) + while(!((as.integer(nb_cyc_k_beg) == obs_nb_cycles) & (as.integer(nb_cyc_k_end) == (obs_nb_cycles - 1))) & ((k_beg != k_beg_prev) | (k_end != k_end_prev))){ + + # Update previous values. + k_beg_prev <- k_beg + k_end_prev <- k_end + + # Case: k_end does not encompass point + # ------[k_beg ------ k_end]------ point ------ + if(as.integer(nb_cyc_k_end) >= obs_nb_cycles){ + + k_beg_valid <- k_end + + k_beg <- k_end + k_end <- k_end*2 + + nb_cyc_k_beg <- nb_cyc_k_end + nb_cyc_k_end <- nb_cycles_k(k_end, N) + + # Case: k_beg does not encompass point + # ------point ------ [k_beg ------ k_end]------ + } else if(as.integer(nb_cyc_k_beg) < obs_nb_cycles){ + + k_end_valid <- k_beg + + k_beg <- k_beg_valid + floor((k_beg-k_beg_valid)/2) + k_end <- k_end_valid + + nb_cyc_k_end <- nb_cyc_k_beg + nb_cyc_k_beg <- nb_cycles_k(k_beg, N) + + # Case: point is encompassed but interval is not tight on the left + # ------[k_beg ------ point ------ k_end]------ + } else if(as.integer(nb_cyc_k_beg) > obs_nb_cycles){ + + k_beg_valid <- max(k_beg, k_beg_valid) + + k_beg <- k_beg + floor((k_end-k_beg)/2) + nb_cyc_k_beg <- nb_cycles_k(k_beg, N) + + # Case: point is encompassed but interval is not tight on the right + # ------[k_beg point ------ k_end]------ + } else if(as.integer(nb_cyc_k_end) < (obs_nb_cycles-1)){ + k_end_valid <- min(k_end, k_end_valid) + + k_end <- k_beg + ceiling((k_end-k_beg)/2) + nb_cyc_k_end <- nb_cycles_k(k_end, N) + } + + #cat(paste(" [", k_beg, " (", nb_cyc_k_beg, "), ", k_end, " (", nb_cyc_k_end, ")] / ", obs_nb_cycles, "\n")) + } + + #cat(paste("RES: [", k_beg, " (", nb_cyc_k_beg, "), ", k_end, " (", nb_cyc_k_end, ")] / ", obs_nb_cycles, "\n")) + list(k_beg=k_beg, k_end=k_end, k_avg=as.integer((k_beg + k_end)/2), nb_cyc_beg=nb_cyc_k_beg, nb_cyc_end=nb_cyc_k_end) +} + +#' This function should be preferred when there are many inversion estimates to be computed. +#' +#' @param n An integer: The number of markers. +#' A marker represents an aligned block, a gene, +#' or any other genomic region being rearranged. +#' +#' @param obs_nb_cycles An integer: The observed number of cycles +#' in the breakpoint graph of two genomes. +#' +#' @author Priscila Biller +inversionEstimate_BD_many <- function(n, obs_nb_cycles, cyc_all=NA){ + # Computes for various values of k the relation: + # Number of inversions k <-> Expected number of cycles after k inversions + if(is.na(cyc_all)){ + cyc_all <- inversionEstimate_BD_setup(n) + } + # Estimates expected number of inversions. + all_nb_cycles <- sapply(cyc_all, '[[', "nb_cyc_k") + closest <- cyc_all[[which.min(abs(all_nb_cycles - obs_nb_cycles))]] + list(k_beg=closest$k, k_end=closest$k, k_avg=closest$k, nb_cyc_beg=closest$nb_cyc_k, nb_cyc_end=closest$nb_cyc_k) +} + +#' Inversion Estimate - Berestycki and Durrett (2006) +#' +#' Computes the expected number of inversions that explains the differences between the gene orders of two genomes. +#' +#' To compute the expected number of inversions, the method needs to first compute the *breakpoint graph*, +#' a classical data structure introduced by Hannenhalli and Pevzner and used in many genome +#' rearrangement problems. +#' +#' Given the number of cycles in the breakpoint graph, +#' it returns the estimated number of evolutionary events (inversions). +#' +#' This function implements Theorem 3 from Berestycki and Durrett (2006), +#' also described in Equation 6 in Biller et al. (2015). +#' +#' @references Berestycki, Nathanaël, and Rick Durrett. "A phase transition in the random transposition random walk." Discrete Mathematics and Theoretical Computer Science. Discrete Mathematics and Theoretical Computer Science, 2003. +#' @references Hannenhalli, Sridhar, and Pavel A. Pevzner. "Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals." Journal of the ACM (JACM) 46.1 (1999): 1-27. +#' @references Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evolution by double cut-and-join." BMC bioinformatics 16.Suppl 14 (2015): S7. +#' +#' @note If the input genomes are multichromosomal, it is recommended to call +#' the function \code{\link{matchPairs}} before this function, to create a mapping +#' between reference and query chromosomes. +#' +#' @param gb A [`GBreaks`] object. +#' @param mode Optional parameter that specifies how estimates are computed. Two options are available: +#' 1. `single` : recommended when computing a few values in one function call; +#' 2. `many` : use when computing many values. +#' @param chrom_sep An optional string that separates the chromosome names in each item of the returned list. +#' +#' @return A list where each item corresponds to a chromosome pair. +#' The name of each item is the matched reference and query chromosome names concatenated, separated by `chrom_sep`. +#' Each item contains the following information: +#' 1. Basic stats of the breakpoint graph: `N`, `nbBreakpoints`, `nbCycles`. See the documentation of `breakpointGraphProperties` for more details; +#' 2. Details of the _Berestycki and Durrett_ estimate stored in `expinv_BD`: +#' * `k_beg` and `k_end`: the lower and upper bounds for the expected number of inversions, respectively; +#' * `k_avg` : the average of `k_beg` and `k_end`; +#' * `nb_cyc_beg` and `nb_cyc_end`: the expected number of cycles in the breakpoint graph after `k_beg` and `k_end` inversions. +#' The observed number of cycles in the breakpoint graph of the two genomes should fall between those values. +#' +#' @examples +#' \dontrun{ +#' +#' # Create a chromosome mapping given a GBreaks object (useful when genomes are multichromosomal). +#' chrMapping <- matchPairs(exampleInversionBader2001) +#' # Compute the expected number of inversions using the method from Berestycki and Durrett (2006). +#' expNbInversions <- inversionEstimate_BD(chrMapping) +#' # Compute the minimum number of inversions using the method from Hannehalli and Pevzner (1999). +#' minNbInversions <- inversionDistance(chrMapping) +#' # Output: Example from Bader et al. (2001): Inversion distance = 7 , Expected nb. of inversions = 7 +#' cat(paste("Example from Bader et al. (2001): Inversion distance =", minNbInversions, ", Expected nb. of inversions =", expNbInversions[[1]]$expinv_BD$k_avg)) +#' +#' # Another example, this time without computing the chromosome mapping. +#' # The chromosome mapping is not needed if genomes are unichromosomal. +#' expNbInversions <- inversionEstimate_BD(exampleInversionBergeron2005b) +#' minNbInversions <- inversionDistance(exampleInversionBergeron2005b) +#' # Output: Example used in the book ``Mathematics of Evolution and Phylogeny`` (2005) (Figure 10.6): Inversion distance = 13 , Expected nb. of inversions = 15 +#' cat(paste("Example used in the book ``Mathematics of Evolution and Phylogeny`` (2005) (Figure 10.6): Inversion distance =", minNbInversions, ", Expected nb. of inversions =", expNbInversions[[1]]$expinv_BD$k_avg)) +#' } +#' +#' @seealso \code{\link{breakpointGraphProperties}} for computing key properties of the breakpoint graph needed for this estimate. +#' +#' @family Rearrangement distances +#' @family Similarity indexes +#' +#' @author Priscila Biller +#' +#' @export +inversionEstimate_BD <- function(gb, mode="single", chrom_sep="---"){ + + # Get matched chromosomes. + chrom_pairs <- lapply(split(gb, paste(seqnames(gb),seqnames(gb$query),sep=chrom_sep), drop=TRUE), function(x) {breakpointGraphProperties(x)}) + + # Sort chromosomes from the one with the smallest number + # of aligned blocks to the one with the highest number. + Ns <- sapply(chrom_pairs, `[[`, "N") + chrom_pairs <- chrom_pairs[order(Ns)] + + N <- -1 + cyc_all <- NA + + for(chrom_pair in names(chrom_pairs)){ + chr <- chrom_pairs[[chrom_pair]] + + # If the chromosome has only one aligned block, there is nothing to be rearranged. + if(chr$N == 1){ + chrom_pairs[[chrom_pair]]$expinv_BD <- NA + next + } + + # Mode: "many" + if(mode == "many"){ + + # If two chromosomes have the same number of markers, + # then the most computationally expensive part of the + # estimate can be reused, saving time. + # Otherwise, a new computation is required. + if(chr$N != N){ + cyc_all <- NA + N <- chr$N + } + + # Estimate expected number of inversions. + chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_many(chr$N, chr$nb_cycles, cyc_all=cyc_all) + + # Mode: "single" + } else { + chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_single(chr$N, chr$nb_cycles) + } + } + chrom_pairs +} From adf31d28760b287d258bdb199d2b71aa4fa6a83e Mon Sep 17 00:00:00 2001 From: pribiller Date: Mon, 17 Aug 2026 14:33:15 +0900 Subject: [PATCH 2/4] Add documentation for new inversion estimate Add documentation for all functions implemented in `inversionEstimate_BD.R`. Update the existing documentation for all functions in the 'Similarity indexes' and 'Rearrangement distances' families to include a reference to the new functions. --- DESCRIPTION | 1 + NAMESPACE | 1 + man/F81_distance.Rd | 1 + man/GOC.Rd | 1 + man/HKY85_distance.Rd | 1 + man/JC69_distance.Rd | 1 + man/K80_distance.Rd | 1 + man/K80_gap_distance.Rd | 1 + man/P_distance.Rd | 1 + man/T92_distance.Rd | 1 + man/TN93_distance.Rd | 1 + man/breakpointGraphProperties.Rd | 4 +- man/correlation_index.Rd | 1 + man/inversionDistance.Rd | 4 +- man/inversionEstimate_BD.Rd | 109 +++++++++++++++++++++++++++++ man/inversionEstimate_BD_many.Rd | 22 ++++++ man/inversionEstimate_BD_setup.Rd | 28 ++++++++ man/inversionEstimate_BD_single.Rd | 27 +++++++ man/karyotype_index.Rd | 1 + man/logDet_distance.Rd | 1 + man/nb_cycles_k.Rd | 30 ++++++++ man/slidingWindow.Rd | 1 + man/strand_randomisation_index.Rd | 1 + man/synteny_index.Rd | 1 + man/tau_index.Rd | 1 + 25 files changed, 240 insertions(+), 2 deletions(-) create mode 100644 man/inversionEstimate_BD.Rd create mode 100644 man/inversionEstimate_BD_many.Rd create mode 100644 man/inversionEstimate_BD_setup.Rd create mode 100644 man/inversionEstimate_BD_single.Rd create mode 100644 man/nb_cycles_k.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 70dca7c7..fbeeaf34 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -95,6 +95,7 @@ Collate: 'int_frag.R' 'inv2UCSCData.R' 'inversionDistance.R' + 'inversionEstimate_BD.R' 'karyotype_index.R' 'keepLongestPair.R' 'load_genomic_breaks.R' diff --git a/NAMESPACE b/NAMESPACE index e4d233f2..7d1c870f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -55,6 +55,7 @@ export(guessSeqLengths) export(int_frag) export(inv2UCSCData) export(inversionDistance) +export(inversionEstimate_BD) export(karyotype_index) export(keepLongestPair) export(leftInversionGaps) diff --git a/man/F81_distance.Rd b/man/F81_distance.Rd index 67ec962c..58873d76 100644 --- a/man/F81_distance.Rd +++ b/man/F81_distance.Rd @@ -64,6 +64,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/GOC.Rd b/man/GOC.Rd index 84f69339..cc07c5e3 100644 --- a/man/GOC.Rd +++ b/man/GOC.Rd @@ -69,6 +69,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/HKY85_distance.Rd b/man/HKY85_distance.Rd index aa2f31a3..b46716e3 100644 --- a/man/HKY85_distance.Rd +++ b/man/HKY85_distance.Rd @@ -87,6 +87,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/JC69_distance.Rd b/man/JC69_distance.Rd index 1761b945..f902a139 100644 --- a/man/JC69_distance.Rd +++ b/man/JC69_distance.Rd @@ -53,6 +53,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/K80_distance.Rd b/man/K80_distance.Rd index 462bec49..976c96d5 100644 --- a/man/K80_distance.Rd +++ b/man/K80_distance.Rd @@ -53,6 +53,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/K80_gap_distance.Rd b/man/K80_gap_distance.Rd index 0db0a8cd..68c1d831 100644 --- a/man/K80_gap_distance.Rd +++ b/man/K80_gap_distance.Rd @@ -85,6 +85,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/P_distance.Rd b/man/P_distance.Rd index d8fc424d..db94bf5e 100644 --- a/man/P_distance.Rd +++ b/man/P_distance.Rd @@ -81,6 +81,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/T92_distance.Rd b/man/T92_distance.Rd index 7623ab23..ff747921 100644 --- a/man/T92_distance.Rd +++ b/man/T92_distance.Rd @@ -58,6 +58,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/TN93_distance.Rd b/man/TN93_distance.Rd index 4e51302e..6e5586fe 100644 --- a/man/TN93_distance.Rd +++ b/man/TN93_distance.Rd @@ -74,6 +74,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/breakpointGraphProperties.Rd b/man/breakpointGraphProperties.Rd index b455e885..e7e2308a 100644 --- a/man/breakpointGraphProperties.Rd +++ b/man/breakpointGraphProperties.Rd @@ -36,7 +36,8 @@ Other Breakpoint graph functions: \code{\link[=superhurdles_count]{superhurdles_count()}} Other Rearrangement distances: -\code{\link[=inversionDistance]{inversionDistance()}} +\code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}} Other Similarity indexes: \code{\link[=F81_distance]{F81_distance()}}, @@ -50,6 +51,7 @@ Other Similarity indexes: \code{\link[=TN93_distance]{TN93_distance()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/correlation_index.Rd b/man/correlation_index.Rd index f831fa43..3e99b787 100644 --- a/man/correlation_index.Rd +++ b/man/correlation_index.Rd @@ -44,6 +44,7 @@ Other Similarity indexes: \code{\link[=TN93_distance]{TN93_distance()}}, \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/inversionDistance.Rd b/man/inversionDistance.Rd index bb4c7f7f..10abe484 100644 --- a/man/inversionDistance.Rd +++ b/man/inversionDistance.Rd @@ -46,7 +46,8 @@ Hannenhalli, Sridhar, and Pavel A. Pevzner. "Transforming cabbage into turnip: p \code{\link{permutationVector}} for generating the permutation vector. Other Rearrangement distances: -\code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}} +\code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}} Other Similarity indexes: \code{\link[=F81_distance]{F81_distance()}}, @@ -60,6 +61,7 @@ Other Similarity indexes: \code{\link[=TN93_distance]{TN93_distance()}}, \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/inversionEstimate_BD.Rd b/man/inversionEstimate_BD.Rd new file mode 100644 index 00000000..3920b62a --- /dev/null +++ b/man/inversionEstimate_BD.Rd @@ -0,0 +1,109 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{inversionEstimate_BD} +\alias{inversionEstimate_BD} +\title{Inversion Estimate - Berestycki and Durrett (2006)} +\usage{ +inversionEstimate_BD(gb, mode = "single", chrom_sep = "---") +} +\arguments{ +\item{gb}{A \code{\link{GBreaks}} object.} + +\item{mode}{Optional parameter that specifies how estimates are computed. Two options are available: +\enumerate{ +\item \code{single} : recommended when computing a few values in one function call; +\item \code{many} : use when computing many values. +}} + +\item{chrom_sep}{An optional string that separates the chromosome names in each item of the returned list.} +} +\value{ +A list where each item corresponds to a chromosome pair. +The name of each item is the matched reference and query chromosome names concatenated, separated by \code{chrom_sep}. +Each item contains the following information: +1. Basic stats of the breakpoint graph: \code{N}, \code{nbBreakpoints}, \code{nbCycles}. See the documentation of \code{breakpointGraphProperties} for more details; +2. Details of the \emph{Berestycki and Durrett} estimate stored in \code{expinv_BD}: +* \code{k_beg} and \code{k_end}: the lower and upper bounds for the expected number of inversions, respectively; +* \code{k_avg} : the average of \code{k_beg} and \code{k_end}; +* \code{nb_cyc_beg} and \code{nb_cyc_end}: the expected number of cycles in the breakpoint graph after \code{k_beg} and \code{k_end} inversions. +The observed number of cycles in the breakpoint graph of the two genomes should fall between those values. +} +\description{ +Computes the expected number of inversions that explains the differences between the gene orders of two genomes. +} +\details{ +To compute the expected number of inversions, the method needs to first compute the \emph{breakpoint graph}, +a classical data structure introduced by Hannenhalli and Pevzner and used in many genome +rearrangement problems. + +Given the number of cycles in the breakpoint graph, +it returns the estimated number of evolutionary events (inversions). + +This function implements Theorem 3 from Berestycki and Durrett (2006), +also described in Equation 6 in Biller et al. (2015). +} +\note{ +If the input genomes are multichromosomal, it is recommended to call +the function \code{\link{matchPairs}} before this function, to create a mapping +between reference and query chromosomes. +} +\examples{ +\dontrun{ + +# Create a chromosome mapping given a GBreaks object (useful when genomes are multichromosomal). +chrMapping <- matchPairs(exampleInversionBader2001) +# Compute the expected number of inversions using the method from Berestycki and Durrett (2006). +expNbInversions <- inversionEstimate_BD(chrMapping) +# Compute the minimum number of inversions using the method from Hannehalli and Pevzner (1999). +minNbInversions <- inversionDistance(chrMapping) +# Output: Example from Bader et al. (2001): Inversion distance = 7 , Expected nb. of inversions = 7 +cat(paste("Example from Bader et al. (2001): Inversion distance =", minNbInversions, ", Expected nb. of inversions =", expNbInversions[[1]]$expinv_BD$k_avg)) + +# Another example, this time without computing the chromosome mapping. +# The chromosome mapping is not needed if genomes are unichromosomal. +expNbInversions <- inversionEstimate_BD(exampleInversionBergeron2005b) +minNbInversions <- inversionDistance(exampleInversionBergeron2005b) +# Output: Example used in the book ``Mathematics of Evolution and Phylogeny`` (2005) (Figure 10.6): Inversion distance = 13 , Expected nb. of inversions = 15 +cat(paste("Example used in the book ``Mathematics of Evolution and Phylogeny`` (2005) (Figure 10.6): Inversion distance =", minNbInversions, ", Expected nb. of inversions =", expNbInversions[[1]]$expinv_BD$k_avg)) +} + +} +\references{ +Berestycki, Nathanaël, and Rick Durrett. "A phase transition in the random transposition random walk." Discrete Mathematics and Theoretical Computer Science. Discrete Mathematics and Theoretical Computer Science, 2003. + +Hannenhalli, Sridhar, and Pavel A. Pevzner. "Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals." Journal of the ACM (JACM) 46.1 (1999): 1-27. + +Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evolution by double cut-and-join." BMC bioinformatics 16.Suppl 14 (2015): S7. +} +\seealso{ +\code{\link{breakpointGraphProperties}} for computing key properties of the breakpoint graph needed for this estimate. + +Other Rearrangement distances: +\code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, +\code{\link[=inversionDistance]{inversionDistance()}} + +Other Similarity indexes: +\code{\link[=F81_distance]{F81_distance()}}, +\code{\link[=GOC]{GOC()}}, +\code{\link[=HKY85_distance]{HKY85_distance()}}, +\code{\link[=JC69_distance]{JC69_distance()}}, +\code{\link[=K80_distance]{K80_distance()}}, +\code{\link[=K80_gap_distance]{K80_gap_distance()}}, +\code{\link[=P_distance]{P_distance()}}, +\code{\link[=T92_distance]{T92_distance()}}, +\code{\link[=TN93_distance]{TN93_distance()}}, +\code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, +\code{\link[=correlation_index]{correlation_index()}}, +\code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=karyotype_index]{karyotype_index()}}, +\code{\link[=logDet_distance]{logDet_distance()}}, +\code{\link[=slidingWindow]{slidingWindow()}}, +\code{\link[=strand_randomisation_index]{strand_randomisation_index()}}, +\code{\link[=synteny_index]{synteny_index()}}, +\code{\link[=tau_index]{tau_index()}} +} +\author{ +Priscila Biller +} +\concept{Rearrangement distances} +\concept{Similarity indexes} diff --git a/man/inversionEstimate_BD_many.Rd b/man/inversionEstimate_BD_many.Rd new file mode 100644 index 00000000..19f7c785 --- /dev/null +++ b/man/inversionEstimate_BD_many.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{inversionEstimate_BD_many} +\alias{inversionEstimate_BD_many} +\title{This function should be preferred when there are many inversion estimates to be computed.} +\usage{ +inversionEstimate_BD_many(n, obs_nb_cycles, cyc_all = NA) +} +\arguments{ +\item{n}{An integer: The number of markers. +A marker represents an aligned block, a gene, +or any other genomic region being rearranged.} + +\item{obs_nb_cycles}{An integer: The observed number of cycles +in the breakpoint graph of two genomes.} +} +\description{ +This function should be preferred when there are many inversion estimates to be computed. +} +\author{ +Priscila Biller +} diff --git a/man/inversionEstimate_BD_setup.Rd b/man/inversionEstimate_BD_setup.Rd new file mode 100644 index 00000000..48a8f6a6 --- /dev/null +++ b/man/inversionEstimate_BD_setup.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{inversionEstimate_BD_setup} +\alias{inversionEstimate_BD_setup} +\title{The setup step computes the expected number of cycles +after k edges for various values of k.} +\usage{ +inversionEstimate_BD_setup(n) +} +\arguments{ +\item{n}{An integer: The number of markers. +A marker represents an aligned block, a gene, +or any other genomic region being rearranged.} +} +\description{ +The expected number of cycles can take any value between 1 +and n (the number of markers). For each of these values, +the method computes at least one k. +} +\details{ +Notice that this is the most time-consuming step +in the estimate. Since it depends only on the number +n of markers, it is recommended that it be computed once +for all permutations of n markers, saving time. +} +\author{ +Priscila Biller +} diff --git a/man/inversionEstimate_BD_single.Rd b/man/inversionEstimate_BD_single.Rd new file mode 100644 index 00000000..cb3b0328 --- /dev/null +++ b/man/inversionEstimate_BD_single.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{inversionEstimate_BD_single} +\alias{inversionEstimate_BD_single} +\title{It finds the number of inversions whose expected number of cycles +is closest to the observed number of cycles.} +\usage{ +inversionEstimate_BD_single(n, obs_nb_cycles) +} +\arguments{ +\item{n}{An integer: The number of markers. +A marker represents an aligned block, a gene, +or any other genomic region being rearranged.} + +\item{obs_nb_cycles}{An integer: The observed number of cycles +in the breakpoint graph of two genomes.} +} +\description{ +This function is an alternative to the function "inversionEstimate_BD_many", +which uses the function "inversionEstimate_BD_setup". +} +\details{ +This function should be preferred when there are few inversion estimates to be computed. +} +\author{ +Priscila Biller +} diff --git a/man/karyotype_index.Rd b/man/karyotype_index.Rd index 07077e3a..80394f0b 100644 --- a/man/karyotype_index.Rd +++ b/man/karyotype_index.Rd @@ -42,6 +42,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, \code{\link[=strand_randomisation_index]{strand_randomisation_index()}}, diff --git a/man/logDet_distance.Rd b/man/logDet_distance.Rd index cbc85897..1ca5f36b 100644 --- a/man/logDet_distance.Rd +++ b/man/logDet_distance.Rd @@ -93,6 +93,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=slidingWindow]{slidingWindow()}}, \code{\link[=strand_randomisation_index]{strand_randomisation_index()}}, diff --git a/man/nb_cycles_k.Rd b/man/nb_cycles_k.Rd new file mode 100644 index 00000000..323aeb80 --- /dev/null +++ b/man/nb_cycles_k.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{nb_cycles_k} +\alias{nb_cycles_k} +\title{Computes the expected number of cycles in a random graph +after k edges.} +\usage{ +nb_cycles_k(k, N) +} +\arguments{ +\item{k}{An integer: The number of edges in the random graph, +equivalent to the number of inversions.} + +\item{N}{An integer: The number of vertices in the random graph, +equivalent to the number of markers+1.} +} +\description{ +Implements Theorem 3 from Berestycki and Durrett (2006), +also described in Equation 6 in Biller et al. (2015). +} +\references{ +Berestycki, Nathanaël, and Rick Durrett. "A phase transition in the random transposition random walk." Discrete Mathematics and Theoretical Computer Science. Discrete Mathematics and Theoretical Computer Science, 2003. + +Hannenhalli, Sridhar, and Pavel A. Pevzner. "Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals." Journal of the ACM (JACM) 46.1 (1999): 1-27. + +Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evolution by double cut-and-join." BMC bioinformatics 16.Suppl 14 (2015): S7. +} +\author{ +Priscila Biller +} diff --git a/man/slidingWindow.Rd b/man/slidingWindow.Rd index 143774d1..923507f4 100644 --- a/man/slidingWindow.Rd +++ b/man/slidingWindow.Rd @@ -54,6 +54,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=strand_randomisation_index]{strand_randomisation_index()}}, diff --git a/man/strand_randomisation_index.Rd b/man/strand_randomisation_index.Rd index e81c6f8a..22a570d8 100644 --- a/man/strand_randomisation_index.Rd +++ b/man/strand_randomisation_index.Rd @@ -62,6 +62,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/synteny_index.Rd b/man/synteny_index.Rd index bbcd31a4..bd69db60 100644 --- a/man/synteny_index.Rd +++ b/man/synteny_index.Rd @@ -52,6 +52,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, diff --git a/man/tau_index.Rd b/man/tau_index.Rd index 6d23580e..d97a01e3 100644 --- a/man/tau_index.Rd +++ b/man/tau_index.Rd @@ -67,6 +67,7 @@ Other Similarity indexes: \code{\link[=breakpointGraphProperties]{breakpointGraphProperties()}}, \code{\link[=correlation_index]{correlation_index()}}, \code{\link[=inversionDistance]{inversionDistance()}}, +\code{\link[=inversionEstimate_BD]{inversionEstimate_BD()}}, \code{\link[=karyotype_index]{karyotype_index()}}, \code{\link[=logDet_distance]{logDet_distance()}}, \code{\link[=slidingWindow]{slidingWindow()}}, From dfc66a3272cf18df7ab6728aac4c575112ece83a Mon Sep 17 00:00:00 2001 From: pribiller Date: Mon, 17 Aug 2026 15:07:43 +0900 Subject: [PATCH 3/4] Fix documentation Fix documentation for some of the new functions: "inversionEstimate_BD_many", "inversionEstimate_BD_setup", "inversionEstimate_BD_single", and "nb_cycles_k". --- NAMESPACE | 4 ++++ R/inversionEstimate_BD.R | 21 ++++++++++++++++++++- man/inversionEstimate_BD_many.Rd | 3 ++- man/inversionEstimate_BD_setup.Rd | 11 +++++++---- man/inversionEstimate_BD_single.Rd | 11 +++++++---- man/nb_cycles_k.Rd | 8 ++++++-- 6 files changed, 46 insertions(+), 12 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 7d1c870f..48f7a7ec 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -56,6 +56,9 @@ export(int_frag) export(inv2UCSCData) export(inversionDistance) export(inversionEstimate_BD) +export(inversionEstimate_BD_many) +export(inversionEstimate_BD_setup) +export(inversionEstimate_BD_single) export(karyotype_index) export(keepLongestPair) export(leftInversionGaps) @@ -67,6 +70,7 @@ export(makeOxfordPlots) export(matchPairs) export(mergeSeqLevels) export(mergeSeqLevels_to_DF) +export(nb_cycles_k) export(orderQuerySeqLevels) export(orderQuerySeqLevels_DF_GR) export(permutationVector) diff --git a/R/inversionEstimate_BD.R b/R/inversionEstimate_BD.R index 3a07525e..699cd2f5 100644 --- a/R/inversionEstimate_BD.R +++ b/R/inversionEstimate_BD.R @@ -1,4 +1,6 @@ +#' Expected number of cycler after k reversals +#' #' Computes the expected number of cycles in a random graph #' after k edges. #' @@ -16,6 +18,9 @@ #' @references Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evolution by double cut-and-join." BMC bioinformatics 16.Suppl 14 (2015): S7. #' #' @author Priscila Biller +#' +#' @export +#' @keywords internal nb_cycles_k <- function(k, N){ # Invalid value (negative number of inversions or not enough markers). @@ -43,6 +48,8 @@ nb_cycles_k <- function(k, N){ } } +#' Inversion Estimate - setup for option 'many' +#' #' The setup step computes the expected number of cycles #' after k edges for various values of k. #' @@ -60,6 +67,9 @@ nb_cycles_k <- function(k, N){ #' or any other genomic region being rearranged. #' #' @author Priscila Biller +#' +#' @export +#' @keywords internal inversionEstimate_BD_setup <- function(n){ # The components of the breakpoint graph for random reversals @@ -103,7 +113,8 @@ inversionEstimate_BD_setup <- function(n){ cyc_all } - +#' Inversion Estimate - option 'single' +#' #' It finds the number of inversions whose expected number of cycles #' is closest to the observed number of cycles. #' @@ -120,6 +131,9 @@ inversionEstimate_BD_setup <- function(n){ #' in the breakpoint graph of two genomes. #' #' @author Priscila Biller +#' +#' @export +#' @keywords internal inversionEstimate_BD_single <- function(n, obs_nb_cycles){ # The components of the breakpoint graph for random reversals @@ -208,6 +222,8 @@ inversionEstimate_BD_single <- function(n, obs_nb_cycles){ list(k_beg=k_beg, k_end=k_end, k_avg=as.integer((k_beg + k_end)/2), nb_cyc_beg=nb_cyc_k_beg, nb_cyc_end=nb_cyc_k_end) } +#' Inversion Estimate - option 'many' +#' #' This function should be preferred when there are many inversion estimates to be computed. #' #' @param n An integer: The number of markers. @@ -218,6 +234,9 @@ inversionEstimate_BD_single <- function(n, obs_nb_cycles){ #' in the breakpoint graph of two genomes. #' #' @author Priscila Biller +#' +#' @export +#' @keywords internal inversionEstimate_BD_many <- function(n, obs_nb_cycles, cyc_all=NA){ # Computes for various values of k the relation: # Number of inversions k <-> Expected number of cycles after k inversions diff --git a/man/inversionEstimate_BD_many.Rd b/man/inversionEstimate_BD_many.Rd index 19f7c785..9bd1a58c 100644 --- a/man/inversionEstimate_BD_many.Rd +++ b/man/inversionEstimate_BD_many.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/inversionEstimate_BD.R \name{inversionEstimate_BD_many} \alias{inversionEstimate_BD_many} -\title{This function should be preferred when there are many inversion estimates to be computed.} +\title{Inversion Estimate - option 'many'} \usage{ inversionEstimate_BD_many(n, obs_nb_cycles, cyc_all = NA) } @@ -20,3 +20,4 @@ This function should be preferred when there are many inversion estimates to be \author{ Priscila Biller } +\keyword{internal} diff --git a/man/inversionEstimate_BD_setup.Rd b/man/inversionEstimate_BD_setup.Rd index 48a8f6a6..5bbc6942 100644 --- a/man/inversionEstimate_BD_setup.Rd +++ b/man/inversionEstimate_BD_setup.Rd @@ -2,8 +2,7 @@ % Please edit documentation in R/inversionEstimate_BD.R \name{inversionEstimate_BD_setup} \alias{inversionEstimate_BD_setup} -\title{The setup step computes the expected number of cycles -after k edges for various values of k.} +\title{Inversion Estimate - setup for option 'many'} \usage{ inversionEstimate_BD_setup(n) } @@ -13,11 +12,14 @@ A marker represents an aligned block, a gene, or any other genomic region being rearranged.} } \description{ +The setup step computes the expected number of cycles +after k edges for various values of k. +} +\details{ The expected number of cycles can take any value between 1 and n (the number of markers). For each of these values, the method computes at least one k. -} -\details{ + Notice that this is the most time-consuming step in the estimate. Since it depends only on the number n of markers, it is recommended that it be computed once @@ -26,3 +28,4 @@ for all permutations of n markers, saving time. \author{ Priscila Biller } +\keyword{internal} diff --git a/man/inversionEstimate_BD_single.Rd b/man/inversionEstimate_BD_single.Rd index cb3b0328..93459384 100644 --- a/man/inversionEstimate_BD_single.Rd +++ b/man/inversionEstimate_BD_single.Rd @@ -2,8 +2,7 @@ % Please edit documentation in R/inversionEstimate_BD.R \name{inversionEstimate_BD_single} \alias{inversionEstimate_BD_single} -\title{It finds the number of inversions whose expected number of cycles -is closest to the observed number of cycles.} +\title{Inversion Estimate - option 'single'} \usage{ inversionEstimate_BD_single(n, obs_nb_cycles) } @@ -16,12 +15,16 @@ or any other genomic region being rearranged.} in the breakpoint graph of two genomes.} } \description{ -This function is an alternative to the function "inversionEstimate_BD_many", -which uses the function "inversionEstimate_BD_setup". +It finds the number of inversions whose expected number of cycles +is closest to the observed number of cycles. } \details{ +This function is an alternative to the function "inversionEstimate_BD_many", +which uses the function "inversionEstimate_BD_setup". + This function should be preferred when there are few inversion estimates to be computed. } \author{ Priscila Biller } +\keyword{internal} diff --git a/man/nb_cycles_k.Rd b/man/nb_cycles_k.Rd index 323aeb80..01c18684 100644 --- a/man/nb_cycles_k.Rd +++ b/man/nb_cycles_k.Rd @@ -2,8 +2,7 @@ % Please edit documentation in R/inversionEstimate_BD.R \name{nb_cycles_k} \alias{nb_cycles_k} -\title{Computes the expected number of cycles in a random graph -after k edges.} +\title{Expected number of cycler after k reversals} \usage{ nb_cycles_k(k, N) } @@ -15,6 +14,10 @@ equivalent to the number of inversions.} equivalent to the number of markers+1.} } \description{ +Computes the expected number of cycles in a random graph +after k edges. +} +\details{ Implements Theorem 3 from Berestycki and Durrett (2006), also described in Equation 6 in Biller et al. (2015). } @@ -28,3 +31,4 @@ Biller, Priscila, Laurent Guéguen, and Eric Tannier. "Moments of genome evoluti \author{ Priscila Biller } +\keyword{internal} From 65e7417bfc4ef851f32996e153e32fbe52ada649 Mon Sep 17 00:00:00 2001 From: pribiller Date: Wed, 19 Aug 2026 14:45:13 +0900 Subject: [PATCH 4/4] Apply proposed suggestions Proposed: - Modifications to the function 'nb_cycles_k' to make it smaller. Other small changes: - Replaces parameter mode="single" to mode="few". - Renames the function 'inversionEstimate_BD_single' to 'inversionEstimate_BD_few'. --- NAMESPACE | 2 +- R/inversionEstimate_BD.R | 51 ++++++++----------- man/inversionEstimate_BD.Rd | 4 +- ..._single.Rd => inversionEstimate_BD_few.Rd} | 8 +-- 4 files changed, 28 insertions(+), 37 deletions(-) rename man/{inversionEstimate_BD_single.Rd => inversionEstimate_BD_few.Rd} (82%) diff --git a/NAMESPACE b/NAMESPACE index 48f7a7ec..61c5a917 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -56,9 +56,9 @@ export(int_frag) export(inv2UCSCData) export(inversionDistance) export(inversionEstimate_BD) +export(inversionEstimate_BD_few) export(inversionEstimate_BD_many) export(inversionEstimate_BD_setup) -export(inversionEstimate_BD_single) export(karyotype_index) export(keepLongestPair) export(leftInversionGaps) diff --git a/R/inversionEstimate_BD.R b/R/inversionEstimate_BD.R index 699cd2f5..e3d227d8 100644 --- a/R/inversionEstimate_BD.R +++ b/R/inversionEstimate_BD.R @@ -22,30 +22,21 @@ #' @export #' @keywords internal nb_cycles_k <- function(k, N){ - - # Invalid value (negative number of inversions or not enough markers). - if ((N < 2) | (k < 0)){ - NA - - # No inversions. Genomes are the same. - } else if(k == 0) { - N - - # At least one inversion. - } else { - - nb_cyc_k <- 0 - # This loop is supposed to go until infinity (sum until infinity), - # but as i grows, the contribution of the i-th term gets smaller and smaller. - # Here, the computation is capped at 10 for efficiency. - # It does not seem to affect much the accuracy. - for(i in 1:10) { - term_base <- 2.0*k/N - term_k <- term_base*exp(-term_base) - nb_cyc_k <- nb_cyc_k + ((1.0 / term_base) * ( i**(i-2) / factorial(i) )* term_k**i) - } - nb_cyc_k*N - } + if ((N < 2) | (k < 0)) return(NA) # Invalid value (negative number of inversions or not enough markers). + if (k == 0) return(N) # No inversions. Genomes are the same + + # At least one inversion. + nb_cyc_k <- 0 + # This loop is supposed to go until infinity (sum until infinity), + # but as i grows, the contribution of the i-th term gets smaller and smaller. + # Here, the computation is capped at 10 for efficiency. + # It does not seem to affect much the accuracy. + for(i in 1:10) { + term_base <- 2.0*k/N + term_k <- term_base*exp(-term_base) + nb_cyc_k <- nb_cyc_k + ((1.0 / term_base) * ( i**(i-2) / factorial(i) )* term_k**i) + } + nb_cyc_k*N } #' Inversion Estimate - setup for option 'many' @@ -113,7 +104,7 @@ inversionEstimate_BD_setup <- function(n){ cyc_all } -#' Inversion Estimate - option 'single' +#' Inversion Estimate - option 'few' #' #' It finds the number of inversions whose expected number of cycles #' is closest to the observed number of cycles. @@ -134,7 +125,7 @@ inversionEstimate_BD_setup <- function(n){ #' #' @export #' @keywords internal -inversionEstimate_BD_single <- function(n, obs_nb_cycles){ +inversionEstimate_BD_few <- function(n, obs_nb_cycles){ # The components of the breakpoint graph for random reversals # on n-1 markers are related to the cycles of random @@ -273,7 +264,7 @@ inversionEstimate_BD_many <- function(n, obs_nb_cycles, cyc_all=NA){ #' #' @param gb A [`GBreaks`] object. #' @param mode Optional parameter that specifies how estimates are computed. Two options are available: -#' 1. `single` : recommended when computing a few values in one function call; +#' 1. `few` : recommended when computing a few values in one function call; #' 2. `many` : use when computing many values. #' @param chrom_sep An optional string that separates the chromosome names in each item of the returned list. #' @@ -315,7 +306,7 @@ inversionEstimate_BD_many <- function(n, obs_nb_cycles, cyc_all=NA){ #' @author Priscila Biller #' #' @export -inversionEstimate_BD <- function(gb, mode="single", chrom_sep="---"){ +inversionEstimate_BD <- function(gb, mode="few", chrom_sep="---"){ # Get matched chromosomes. chrom_pairs <- lapply(split(gb, paste(seqnames(gb),seqnames(gb$query),sep=chrom_sep), drop=TRUE), function(x) {breakpointGraphProperties(x)}) @@ -352,9 +343,9 @@ inversionEstimate_BD <- function(gb, mode="single", chrom_sep="---"){ # Estimate expected number of inversions. chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_many(chr$N, chr$nb_cycles, cyc_all=cyc_all) - # Mode: "single" + # Mode: "few" } else { - chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_single(chr$N, chr$nb_cycles) + chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_few(chr$N, chr$nb_cycles) } } chrom_pairs diff --git a/man/inversionEstimate_BD.Rd b/man/inversionEstimate_BD.Rd index 3920b62a..19fe2f85 100644 --- a/man/inversionEstimate_BD.Rd +++ b/man/inversionEstimate_BD.Rd @@ -4,14 +4,14 @@ \alias{inversionEstimate_BD} \title{Inversion Estimate - Berestycki and Durrett (2006)} \usage{ -inversionEstimate_BD(gb, mode = "single", chrom_sep = "---") +inversionEstimate_BD(gb, mode = "few", chrom_sep = "---") } \arguments{ \item{gb}{A \code{\link{GBreaks}} object.} \item{mode}{Optional parameter that specifies how estimates are computed. Two options are available: \enumerate{ -\item \code{single} : recommended when computing a few values in one function call; +\item \code{few} : recommended when computing a few values in one function call; \item \code{many} : use when computing many values. }} diff --git a/man/inversionEstimate_BD_single.Rd b/man/inversionEstimate_BD_few.Rd similarity index 82% rename from man/inversionEstimate_BD_single.Rd rename to man/inversionEstimate_BD_few.Rd index 93459384..82f59551 100644 --- a/man/inversionEstimate_BD_single.Rd +++ b/man/inversionEstimate_BD_few.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/inversionEstimate_BD.R -\name{inversionEstimate_BD_single} -\alias{inversionEstimate_BD_single} -\title{Inversion Estimate - option 'single'} +\name{inversionEstimate_BD_few} +\alias{inversionEstimate_BD_few} +\title{Inversion Estimate - option 'few'} \usage{ -inversionEstimate_BD_single(n, obs_nb_cycles) +inversionEstimate_BD_few(n, obs_nb_cycles) } \arguments{ \item{n}{An integer: The number of markers.