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..61c5a917 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -55,6 +55,10 @@ export(guessSeqLengths) export(int_frag) export(inv2UCSCData) export(inversionDistance) +export(inversionEstimate_BD) +export(inversionEstimate_BD_few) +export(inversionEstimate_BD_many) +export(inversionEstimate_BD_setup) export(karyotype_index) export(keepLongestPair) export(leftInversionGaps) @@ -66,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 new file mode 100644 index 00000000..e3d227d8 --- /dev/null +++ b/R/inversionEstimate_BD.R @@ -0,0 +1,352 @@ + +#' Expected number of cycler after k reversals +#' +#' 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 +#' +#' @export +#' @keywords internal +nb_cycles_k <- function(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' +#' +#' 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 +#' +#' @export +#' @keywords internal +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 +} + +#' Inversion Estimate - option 'few' +#' +#' 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 +#' +#' @export +#' @keywords internal +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 + # 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) +} + +#' 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. +#' 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 +#' +#' @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 + 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. `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. +#' +#' @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="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)}) + + # 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: "few" + } else { + chrom_pairs[[chrom_pair]]$expinv_BD <- inversionEstimate_BD_few(chr$N, chr$nb_cycles) + } + } + chrom_pairs +} 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..19fe2f85 --- /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 = "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{few} : 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_few.Rd b/man/inversionEstimate_BD_few.Rd new file mode 100644 index 00000000..82f59551 --- /dev/null +++ b/man/inversionEstimate_BD_few.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inversionEstimate_BD.R +\name{inversionEstimate_BD_few} +\alias{inversionEstimate_BD_few} +\title{Inversion Estimate - option 'few'} +\usage{ +inversionEstimate_BD_few(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{ +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/inversionEstimate_BD_many.Rd b/man/inversionEstimate_BD_many.Rd new file mode 100644 index 00000000..9bd1a58c --- /dev/null +++ b/man/inversionEstimate_BD_many.Rd @@ -0,0 +1,23 @@ +% 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{Inversion Estimate - option 'many'} +\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 +} +\keyword{internal} diff --git a/man/inversionEstimate_BD_setup.Rd b/man/inversionEstimate_BD_setup.Rd new file mode 100644 index 00000000..5bbc6942 --- /dev/null +++ b/man/inversionEstimate_BD_setup.Rd @@ -0,0 +1,31 @@ +% 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{Inversion Estimate - setup for option 'many'} +\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 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. + +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 +} +\keyword{internal} 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..01c18684 --- /dev/null +++ b/man/nb_cycles_k.Rd @@ -0,0 +1,34 @@ +% 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{Expected number of cycler after k reversals} +\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{ +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). +} +\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 +} +\keyword{internal} 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()}},