R/Functions.R

In isqg: In Silico Quantitative Genetics

## -*- mode: R -*-
###########################################################################
## This file is part of isqg: R package for in silico quantitative genetics
##
##              Copyright (C) 2018 Fernando H. Toledo CIMMYT
##
## * Filename: Functions.R
##
## * Description: Defines package functionalities (standalone methods)
##
## * Author: Fernando H. Toledo
##
## * Maintainer: Fernando H. Toledo
##
## * Created: Fr Mar 09 2018
##
##   This program is free software; you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##   This program is distributed in the hope that it will be useful, but
##   WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
##   General Public License for more details.
##
##   You should have received a copy of the GNU General Public License
##   along with this program; if not, write to the Free Software Foundation,
##   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
##
##   `` Far better an approximate answer to the right question, which is
##   often vague, than the exact answer to the wrong question, which can
##   always be made precise ''
##                          --John Tukey, Ann. Math. Stat. 33(1):13 1962
##
###########################################################################

##' @export
##' @title Constructor of Instances of the Specie Class
##'
##' @name set_specie
##'
##' @description Constructor of instances of the Specie class given the map of
##'     the genome and optionally a pointer to a C++ function which will drive
##'     the meiosis process.
##'
##' @details By standard the meiosis recombination and de novo genetic variability
##' is generated by means a count-location process \insertCite{karlin1978}{isqg}.
##'
##' @references
##'     \insertAllCited{}
##'
##' @return Objects of R6 class with methods to mimic in silico Genomes.
##'
##' @param data    A data frame with the map of the Genome to be simulates.
##' @param meiosis A pointer to a C++ function of the meiosis process.
##'
##' @examples
##' data(ToyMap)
##' spc_standard <- set_specie(ToyMap)
##'
##' ## generate standard _de novo_ variability
##' spc_standard$gamete(n = 100)
##'
##' \dontrun{
##' ## write your function in C++ and then wrap it as a pointer
##' ## check the examples in extdata
##' ## compile the code
##' Rcpp::sourceCpp(file = system.file("extdata", "Independent.cpp", package = "isqg"),
##'                 rebuild = TRUE)
##'
##' ## define a specie w/ custom meiosis
##' spp_custom <- set_specie(ToyMap, meiosis = indepp())
##'
##' ## check meiosis process
##' spp_custom$gamete(n = 100)
##' }
##'
##' @rdname set_specie
"set_specie" <- function(data, meiosis = NULL) {
  data_order <- data[order(data$chr, data$pos),]    # reorder data
  list_map <- split(data_order$pos, data_order$chr) # chromosomes
  snp <- data_order$snp                             # snp names
  pts <- unname(do.call(c, list_map))               # all loci
  chr <- data_order$chr - 1                         # C++ indices
  pos <- unname(do.call(c, lapply(list_map, function(x) order(x, decreasing = TRUE) - 1)))
  ## finding chromosome limits
  lwr <- findInterval(0:(length(list_map) - 1), chr, left.open = TRUE)
  upr <- findInterval(0:(length(list_map) - 1), chr) - 1
  if (methods::is(meiosis, "externalptr")) { # switch custom/standard
    return(.Cpp_Specie_cus_ctor(list_map, snp, chr, pts, pos, lwr, upr, meiosis))
  } else {
    return(.Cpp_Specie_cus_ctor(list_map, snp, chr, pts, pos, lwr, upr, .Cpp_meiosis_standard()))
  }
}

##' @export
##' @title Constructor of a Founder Instances of the Specimen Class
##'
##' @name founder
##'
##' @description Constructor of instances of the Specimen class given the Specie
##'     from which the individual will belong where all loci will equal to the
##'     provided genotype.
##'
##' @return Objects of R6 class with methods to mimic in silico Specimens.
##'
##' @param specie an instance of the R6 class Specie with the genome's parameters.
##' @param code   a length one character vector with one of the genotype codes:
##'     "AA", "Aa", "aA" or "aa".
##'
##' @details Genotypes can be coded as \strong{AA}, \strong{Aa}, \strong{aA} or
##'     \strong{aa}, that meant to represent both homozigous (\strong{AA} and
##'     \strong{aa}) as well as both heterozigous (\strong{Aa} and \strong{aA}).
##'
##' @examples
##' data(ToyMap)
##' spc <- set_specie(ToyMap)
##'
##' ## through standalone function
##' AA <- founder(spc, "AA")
##' aa <- founder(spc, "aa")
##'
##' ## or by the Specie's method
##' Aa <- spc$founder("Aa")
##' aA <- spc$founder("aA")
##'
##' @rdname Founder
'founder' <- function(specie, code) { return(.Cpp_founder_ctor(specie, code)) }

##' @title fletch and melt ith homologous
##'
##' @keywords internal
'.fuse' <- function(i, ref) paste0(sapply(ref, substr, start = i, stop = i), collapse = '')

##' @title decode character vector genotype as homologous strings
##'
##' @keywords internal
'.decode' <- function(vec) {
    if (any(!grepl("^[12]\\s[12]$", vec)))
        stop("Invalid code found")
    binary <- gsub("\\s", "", gsub("2", 0, vec))
    pair <- sapply(1:2, .fuse, ref = binary)
    return(pair)
}

##' @export
##' @title Constructor of a Custom Instances of the Specimen Class
##'
##' @name import
##'
##' @description Constructor of instances of the Specimen class given the Specie
##'     from which the individual will belong where the loci will equal to the
##'     provided genotype from two strings one for each homologous.
##'
##' @return Objects of R6 class with methods to mimic in silico Specimens.
##'
##' @param specie   an instance of the R6 class Specie with the genome's parameters.
##' @param genotype a named character vector with the coded/phased genotypes.
##'
##' @examples
##' data(ToyMap)
##' spc <- set_specie(ToyMap)
##'
##' ## simulating what is very close to your real genotypes
##' Real <- sample(c('2 2', '2 1', '1 2', '1 1'), size = nrow(ToyMap), replace = TRUE)
##' names(Real) <- ToyMap$snp # ensure snp names!
##'
##' ## now you can play _in silico_
##' Virtual <- import(spc, Real)
##' S1 <- Virtual$selfcross(n = 10)
##'
##' @rdname Import
'import' <- function(specie, genotype) {
  snps <- .Cpp_spc_snps(specie)
  if (!all.equal(snps, names(genotype)))
    stop( "Provided genotypes doesn't belong to the provided Specie" )
  homologous <- .decode(genotype)
  return(.Cpp_import_ctor(specie, homologous[1], homologous[2]))
}

##' @export
##' @title Codify Specimens' Genotypes.
##'
##' @name genotype
##'
##' @description Codify Specimens' genotypes instances as numeric codes [-1/0/1]
##'     or as character vector that keeps the phase information.
##'
##' @return A numeric or character matrix with the codified Specimens' genotypes.
##'
##' @param pop   a list with instances of the R6 class Specimen.
##' @param phase logical should the codes keep the phase.
##'
##' @examples
##' data(ToyMap)
##' spc <- set_specie(ToyMap)
##'
##' Aa <- founder(spc, "Aa")
##' aA <- spc$founder("aA")
##'
##' Both <- list(Aa = Aa, aA = aA)
##'
##' ## different ways
##' genotype(Both)               # as numeric
##' genotype(Both, phase = TRUE) # as character
##'
##' @rdname Genotype
'genotype' <- function(pop, phase = FALSE) {
  if (phase) {
    codes <- sapply(pop, .Cpp_Genotype_cod)
    rownames(codes) <- .Cpp_get_snps(pop[[1]])
    return(codes)
  } else {
    codes <- sapply(pop, .Cpp_Genotype_num)
    rownames(codes) <- .Cpp_get_snps(pop[[1]])
    return(codes)
  }
}

## \EOF
###########################################################################