Оптимизировать алгоритм, который рассчитывает ожидаемое богатство аллелей на разных уровнях подвыборки

Question

Я реализовал алгоритм для расчета богатства аллелей на основе формул 1 - 3, представленных в:

• Подсчет аллелей с разрежением: частные аллели и иерархические образцы выборки - Стивен Т. Калиновский, Ссылка на PDF

и та же формула из:

• Различающиеся тенденции между гетерозиготностью и аллельным богатством во время послеледниковой колонизации в европейском буке - B. Comps, Ссылка на бумагу

и нужна помощь в полной его оптимизации.

Богатство аллелей является мерой генетического разнообразия, но на него влияет размер выборки. Эта формула позволяет нам оценить ожидаемое богатство аллелей при меньших размерах выборки без повторной выборки. Я использую его для оценки богатства аллелей при всех возможных размерах подвыборки, что позволяет мне строить кривые разрежения.

В качестве первого шага в алгоритме я вычисляю вероятность не наблюдения аллеля на каждом уровне подсчета на каждом уровне подвыборки, чтобы составить справочную таблицу для расчета фактических вероятностей. Я рассчитываю как можно меньше значений (я думаю), используя тот факт, что многие значения равны 1 - (предыдущие вычисленные значения). Это все еще самая медленная часть, но я думаю, что получил масштабирование при n * log (n). В основном я хочу знать, существует ли более эффективный способ создания векторов и объединения их в таблицу (фрейм данных).

Второй шаг - использование таблицы соответствия для расчета ожидаемого богатства аллелей на каждом уровне подвыборки. Я изменил эту часть на более быструю реализацию и обновил код ниже.

Вот код, который есть на данный момент. Вы можете найти его на моем GitHub, DiDeoxy / PGDA: calc_allele_richness.R .

Для запуска кода вы можете установить пакет с: devtools::install_github("https://github.com/DiDeoxy/PGDA") и использовать его с library(PGDA)

#' Calculate mean allele richness at all sampling levels
#'
#' Calculates the mean allele richness across all markers for a sample at all
#' sampling levels. Missing data is not allowed. Based on the formula presented
#' in https://www.genetics.org/content/157/1/389
#'
#' @param pop a data frame with individuals in columns and markers in rows,
#' there must be atleast two individuals
#' @param allele_coding the coding used for indicating the different alleles
#' @param num_cores the number of cores to use, must be 1 on windows, can use 
#' detectCores() of the parallel package on linux
#'
#' @importFrom magrittr %>%
#' @importFrom parallel mclapply
#' @importFrom scrime rowTables
#' 
#' @return a table of expected allele richness for each marker at each 
#' subsampling level with markers in rows and sampling levels in columns
#'
#' @export

allele_richness <- function (pop, allele_coding = 1:2, num_cores = 1) {
  # the total number of alleles observed at each marker
  n <- ncol(pop)
  # probs contains the probability of not observing allele i at each 
  # sub-sampling level (n - k) for each possible count of allele i with
  # allele count in rows and k in columns
  #
  # for each subsampling level
  probs <- mclapply(0:(n - 1), function (k) {
    # a vector for containng the probs of not observing allele i at each count
    # level at each subsampling (n - k) level
    inter <- rep(0, n)
    # if n - k <= 1 then the prob of not observing allele i is 0 at all count
    # levels, the smaller k is compared to n the more levels will have probs of
    # not observing allel i greater than 0
    if (n - k > 1) {
      # probs of not observing allele i are linear decreasing, therefore the top
      # half and bottom half are 1 - mirrors, we can use this fact to skip a lot
      # of computation
      temp <- lapply(1:floor((n - k) / 2), function (n_i) {
        (n - n_i - k) / (n - k)
      }) %>% do.call(c, .)
      # concatenating the calced probs with their 1 - mirror, if n - k is odd
      # the middle value will equal 0.5 which we do not need to mirror
      temp <- c(temp, rev(1 - temp[which(temp != 0.5)]))
      inter[1:length(temp)] <- temp
      inter
    } else {
      inter
    }
  }, mc.cores = num_cores) %>% do.call(cbind, .)
  # creates a data frame containg the counts of each allele for each marker
  marker_allele_counts <- rowTables(pop, allele_coding)
  # we calcuate the mean allele richness across all markers at each subsampling
  # level (n - k) by calculating the product of not observing each allele at
  # each sub-sampling level then taking the sum of these for each marker and 
  # then taking the mean across all markers
  #
  # for each marker
  mclapply(1:nrow(marker_allele_counts), function (marker) {
    (1 - lapply(1:length(marker_allele_counts[marker, ]), function (allele) {
      # for each allele, calc the probability of not observing the allele at
      # each sub-sampling level
      cumprod(probs[marker_allele_counts[[marker, allele]], ])
    # rbind the probabilities for each allele at each sub-smapling level,
    # subtract from one to turn them into probabilities of observing the allele,
    # and sum the alleles together
    }) %>% do.call(rbind, .)) %>% colSums()
  # return a table with markers in rows and sub-sampling levels in columns
  }, mc.cores = num_cores) %>% do.call(rbind, .)
}

Спасибо за любую помощь, которую вы можете оказать, я впервые запрограммировал что-то подобное.

Приветствия

дидезоксите.