R/thetaCorrection.R

Defines functions likTheta thetaCorr

# v: vector of (integer) alleles; typically all founder alleles
# afr: output of afreq(m)
# theta: correction
thetaCorr = function(v, afr, theta) {
  p = sapply(seq_along(v), function(j) {
    a = v[j]
    bj = sum(v[seq_len(j - 1)] == a)
    (bj * theta + (1 - theta) * afr[a]) / (1 + (j - 2)*theta)
  })
  p
}


likTheta = function(x, m, theta, peeler, peelOrder) {

  #---------------------
  # Mostly for debug (peeler and peelOrder will not be missing in practice)
  if(missing(peeler)) {
    if(!isXmarker(m))
      peeler = function(dat, sub) .peel_M_AUT(dat, sub, mutmat = mutmod(m))
    else
      peeler = function(dat, sub) .peel_M_X(dat, sub, SEX = x$SEX, mutmat = mutmod(m))
  }

  if(missing(peelOrder)) {
    # NB: Don't use `informativeSubnucs()` here: pedigree cannot be trimmed (at least not upwards).
    peelOrder = peelingOrder(x)
  }

  #----------------------

  # Freqs
  afr = afreq(m)

  # Starting point: all possible genotypes for each
  # NB: no symmetry reduction for founders
  glist = .buildGenolist(x, m, eliminate = 10, foundersUnordered = FALSE)

  # Number of possible genotypes for each
  nGeno = sapply(glist, function(g) length(g$pat))

  # Founders (internal ID)
  FOU = c(founders(x, internal = TRUE), attr(peelOrder, "treatAsFounder"))
  NONFOU = seq_len(pedsize(x))[-FOU]

  # List of all combinations of founder genotypes
  fouGrid = fastGrid(lapply(nGeno[FOU], seq_len), as.list = TRUE)

  # Useful later
  glistFOU = glist[FOU]

  # Initialise likelihood
  likel = 0

  # Loop over combos
  for (gIdx in fouGrid) { # gIdx is a vector of indices, of length same as FOU

    # Collect all founder alleles in that combo
    als = unlist(lapply(seq_along(FOU), function(i) {
      g = glistFOU[[i]]
      idx = gIdx[i]
      c(g$pat[idx], g$mat[idx])
    }))

    # Theta-corrected probs (when als are sampled sequentially)
    probs = thetaCorr(als, afr, theta)

    # Create startdata with single geno for each founder
    # TODO: better to use prod(probs) in first and 1,1,.. for rest?!
    gli = glist
    for(i in seq_along(FOU)) {
      idx = FOU[i]
      gli[[idx]] = list(pat = als[2*i-1], mat = als[2*i],
                        prob = probs[2*i-1]*probs[2*i])
    }
    for(j in NONFOU)
      gli[[j]]$prob = rep(1, nGeno[j])

    attr(gli, "impossible") = any(probs == 0)

    li = peelingProcess(x, m, startdata=gli, peeler, peelOrder)
    likel = likel + li
  }

  likel
}

Try the pedprobr package in your browser

Any scripts or data that you put into this service are public.

pedprobr documentation built on April 14, 2023, 12:31 a.m.