sequoia (version 2.3.5)

EstConf: Confidence Probabilities

Description

Estimate confidence probabilities ('backward') and assignment error rates ('forward') per category (genotyped/dummy) by repeatedly simulating genotype data from a reference pedigree using SimGeno, reconstruction a pedigree from this using sequoia, and counting the number of mismatches using PedCompare.

Usage

EstConf(
  Pedigree = NULL,
  LifeHistData = NULL,
  args.sim = list(nSnp = 400, SnpError = 0.001, ParMis = c(0.4, 0.4)),
  args.seq = list(Module = "ped", Err = 0.001, Tassign = 0.5, CalcLLR = FALSE),
  nSim = 10,
  nCores = 1,
  quiet = TRUE
)

Arguments

Pedigree

reference pedigree from which to simulate, dataframe with columns id-dam-sire. Additional columns are ignored.

LifeHistData

dataframe with id, sex (1=female, 2=male, 3=unknown), and birth year.

args.sim

list of arguments to pass to SimGeno, such as nSnp (number of SNPs), SnpError (genotyping error rate) and ParMis (proportion of non-genotyped parents). Set to NULL to use all default values.

args.seq

list of arguments to pass to sequoia, such as Module ('par' or 'ped'), Err (assumed genotyping error rate), and Complex. May include (part of) SeqList, a list of sequoia output (i.e. as a list-within-a-list). Set to NULL to use all default values.

nSim

number of iterations of simulate - reconstruct - compare to perform, i.e. number of simulated datasets.

nCores

number of computer cores to use. If >1, package parallel is used. Set to NULL to use all but one of the available cores, as detected by parallel::detectCores() (using all cores tends to freeze up your computer).

quiet

suppress messages. TRUE runs SimGeno and sequoia quietly, 'very' also suppresses other messages and the iteration counter when nCores=1 (there is no iteration counter when nCores>1).

Value

A list, with elements:

ConfProb

See below

PedErrors

See below

Pedigree.reference

the pedigree from which data was simulated

LifeHistData

Pedigree.inferred

a list with for each iteration the inferred pedigree based on the simulated data

SimSNPd

a list with for each iteration the IDs of the individuals simulated to have been genotyped

PedComp.fwd

Counts from the 'forward' PedCompare, from which PedErrors is calculated

RunParams

a list with the call to EstConf, as well as the default parameter values for SimGeno, and sequoia.

RunTime

sequoia runtime per simulation in seconds, as measured by system.time()['elapsed'].

Dataframe ConfProb has 7 columns:

id.cat, dam.cat, sire.cat

Category of the focal individual, dam, and sire, in the pedigree inferred based on the simulated data. Coded as G=genotyped, D=dummy, X=none

dam.conf

Probability that the dam is correct, given the categories of the assigned dam and sire (ignoring whether or not the sire is correct)

sire.conf

as dam.conf, for the sire

pair.conf

Probability that both dam and sire are correct, given their categories

N

Number of individuals per category-combination, across all nSim iterations

Array PedErrors has three dimensions:

class

• FalseNeg(atives): could have been assigned but was not (individual + parent both genotyped or dummyfiable; P1only in PedCompare).

• FalsePos(itives): no parent in reference pedigree, but one was assigned based on the simulated data (P2only)

• Mismatch: different parents between the pedigrees

cat

Category of individual + parent, as a two-letter code where the first letter indicates the focal individual and the second the parent; G=Genotyped, D=Dummy, T=Total

parent

dam or sire

Assumptions

Because the actual true pedigree is (typically) unknown, the provided reference pedigree is used as a stand-in and assumed to be the true pedigree, with unrelated founders. It is also assumed that the probability to be genotyped is equal for all parents; in each iteration, a new random set of parents (proportion set by ParMis) is mimicked to be non-genotyped. In addition, SNPs are assumed to segregate independently.

Object size

The size in Kb of the returned list can become pretty big, as each of the inferred pedigrees is included. When running EstConf many times for a range of parameter values, it may be prudent to save the required summary statistics for each run rather than the full output.

Details

The confidence probability is taken as the number of correct (matching) assignments, divided by all assignments made in the observed (inferred-from-simulated) pedigree. In contrast, the false negative & false positive assignment rates are proportions of the number of parents in the true (reference) pedigree. Each rate is calculated separatedly for dams & sires, and separately for each category (Genotyped/Dummy(fiable)/X (none)) of individual, parent and co-parent.

This function does not know which individuals in the actual Pedigree are genotyped, so the confidence probabilities need to be added to the Pedigree as shown in the example at the bottom.

A confidence of \(1\) means all assignments on simulated data were correct for that category-combination. It should be interpreted as (and perhaps modified to) \(> 1 - 1/N\), where sample size N is given in the last column of the ConfProb and PedErrors dataframes in the output. The same applies for a false negative/positive rate of \(0\) (i.e. to be interpreted as \(< 1/N\)).

See Also

SimGeno, sequoia, PedCompare.

Examples

# NOT RUN {
data(Ped_HSg5, LH_HSg5, package="sequoia")

## Example A: parentage assignment only
conf.A <- EstConf(Pedigree = Ped_HSg5, LifeHistData = LH_HSg5,
   args.sim = list(nSnp = 100, SnpError = 5e-3, ParMis=c(0.2, 0.5)),
   args.seq = list(Module="par", Err=1e-3, Tassign=0.5), nSim = 3)

# parent-pair confidence, per category:
conf.A$ConfProb

# calculate (correct) assignment rates (ignores co-parent)
1 - apply(conf.A$PedErrors, c(1,3), sum, na.rm=TRUE)

## Example B: with sibship clustering, based on sequoia inferred pedigree
RealGenotypes <- SimGeno(Ped = Ped_HSg5, nSnp = 100,
                         ParMis=c(0.19,0.53), SnpError = 6e-3)
SeqOUT <- sequoia(GenoM = RealGenotypes,
                  LifeHistData = LH_HSg5,
                  Err=5e-3, Module="ped",
                  quiet=TRUE, Plot=FALSE)

conf.B <- EstConf(Pedigree = SeqOUT$Pedigree,
              LifeHistData = LH_HSg5,
               args.sim = list(nSnp = 100, SnpError = 5e-3,
                               ParMis=c(0.2, 0.5)),
              args.seq = list(Err=5e-3, Module="ped"),
              nSim = 2, nCores=2)
conf.B$ConfProb

Ped.withConf <- getAssignCat(Pedigree = SeqOUT$Pedigree,
                             SNPd = rownames(RealGenotypes))
Ped.withConf <- merge(Ped.withConf, conf.B$ConfProb, all.x=TRUE, sort=FALSE)
Ped.withConf <- Ped.withConf[, c("id","dam","sire", "dam.conf", "sire.conf",
                                 "id.cat", "dam.cat", "sire.cat")]
head(Ped.withConf[Ped.withConf$dam.cat=="G", ])
head(Ped.withConf[Ped.withConf$dam.cat=="D", ])


## P(actual FS | inferred as FS) etc.
PairL <- list()
for (i in 1:length(conf.A$Pedigree.inferred)) {  # nSim
  cat(i, "\t")
  PairL[[i]] <- ComparePairs(conf.A$Pedigree.reference,
                             conf.A$Pedigree.inferred[[i]],
                             GenBack=1, patmat=TRUE, ExcludeDummies = TRUE,
                             Return="Counts")
}
# P(actual relationship (Ped1) | inferred relationship (Ped2))
PairA <- plyr::laply(PairL, function(M) sweep(M, 2, colSums(M), "/"))
PairRel.prop <- apply(PairA, 2:3, mean, na.rm=TRUE)  # mean across simulations
round(PairRel.prop, 2)
#' # or: P(inferred relationship | actual relationship)
PairA2 <- plyr::laply(PairL, function(M) sweep(M, 1, rowSums(M), "/"))
# }
# NOT RUN {
# }