Interactive Use of nprcgenekeepr

Direct Edge Style

## edgeStyle is now explicit -- Track 2 flipped the package's own default
## to "rectilinear", so this section's own "Direct Edge Style" labeling
## needs edgeStyle = "direct" spelled out rather than relying on the
## (now different) implicit default.
diagramData <- makePedigreeMatingLayout(demoPed, edgeStyle = "direct")
names(diagramData)

## [1] "nodes"              "edges"              "duplicateToReal"   
## [4] "isolatedIds"        "conditions"         "conditionsNotDrawn"

nrow(diagramData$nodes)

## [1] 56

nrow(diagramData$edges)

## [1] 61

makePedigreeMatingLayout returns a list of four elements: nodes (one row per real animal that has recorded parents, mates, or offspring – shaped by sex, dot = Female, square = Male, star = Hermaphrodite, triangle = Unknown, diamond = Other/Unrecorded, with an HTML hover-tooltip giving ID, sex, generation, sire, and dam – plus one small, unlabeled dot per mating and one extra row per duplicate occurrence, each carrying its animal’s own shape, label, and tooltip with a note that it is a duplicate occurrence), edges (parent-to-mating and mating-to-child edges, plus a dashed edge from each duplicate occurrence back to the real animal), duplicateToReal (a named lookup from each duplicate occurrence’s id back to the real animal it represents), and isolatedIds (the ids of animals with no recorded parents, mates, or offspring, which are left out of the diagram; the function also prints a message naming them).

Each node already carries its own fixed x/y position (computed by makePedigreeMatingLayout itself), so rendering it turns vis.js’s usual physics/hierarchical layout off rather than turning one on. drawPedigreeDiagram draws diagramData the way the Diagram tab does – the tab calls the same function – so it reproduces the Diagram tab’s rendering exactly: the legend, a slash through the symbol of an animal whose status is DECEASED, the colored sections of an animal with more than one affected condition, and a genuinely working Export Diagram (PNG) button – try it below. Give it the edgeStyle the layout was made with.⁵

drawPedigreeDiagram(diagramData, edgeStyle = "direct")

The Select by id dropdown above the diagram – filtered to real animals, excluding mating and duplicate-occurrence nodes – lets you jump straight to a specific animal and dims every node except it and its direct connections – useful for locating one animal in a diagram this size, and just as useful in a real diagram of 704 animals like the one built earlier in this tutorial, where finding one animal by eye is much harder.

Rectilinear Edge Style

The diagram above uses the direct edge style, requested explicitly. Omitting edgeStyle entirely (or passing edgeStyle = "rectilinear" explicitly, as below) draws the same relationships – same demoPed, same animals, same parent/child/mating structure – with strict right angles: a horizontal line directly between mates, a vertical drop to their shared mating dot, a horizontal bar across their children, and a vertical drop into each child, matching the convention the kinship2 R package uses, and matching makePedigreeMatingLayout’s own default (Track 2).

diagramDataRectilinear <- makePedigreeMatingLayout(demoPed,
                                                     edgeStyle = "rectilinear")
names(diagramDataRectilinear)

## [1] "nodes"              "edges"              "duplicateToReal"   
## [4] "isolatedIds"        "conditions"         "conditionsNotDrawn"

nrow(diagramDataRectilinear$nodes)

## [1] 94

nrow(diagramDataRectilinear$edges)

## [1] 99

94 nodes and 99 edges here, versus 56 and 61 for the same 33-animal demoPed under the direct style above – the difference is entirely extra invisible “waypoint” nodes/edges that carry out the right-angle routing (zero size, transparent color, excluded from the Select by id dropdown below just like the mating dots above), not any change to which animals or relationships are shown.⁶ nodes also gains an extra color.border column under rectilinear – visNetwork::visNetwork() reads it directly, so no extra rendering code is needed.

Rendering it is the same call with edgeStyle = "rectilinear", which makes one difference: drawPedigreeDiagram raises the hover highlight’s reach (highlightNearest’s degree) from 1 to 6. Under rectilinear, a plain individual’s nearest graph neighbor is often one of the invisible waypoint nodes above rather than a visible mating dot, so degree = 1 (correct for direct) can highlight nothing visible at all on hover; 6 restores visible feedback for the common case. It is a deliberate compromise, not a complete fix, because the invisible waypoints count as steps: on the bundled 375-animal rhesus pedigree, hovering an animal lights everything that the direct style’s degree = 1 lights for all but 2 animals, whose parents’ mating dot is just out of reach (their family’s connectors were rerouted around crowded lines, which adds waypoints). One mating pair with nine or more full siblings reaches the same limit; the widest pair in the bundled pedigrees has four. A larger degree would cover these cases, but it highlights more of the diagram on every hover, so the simpler, tighter hover was kept. The live Shiny app’s Diagram tab applies this same style-aware degree automatically when you switch the “Diagram Edge Style” toggle.

drawPedigreeDiagram(diagramDataRectilinear, edgeStyle = "rectilinear")

The live Shiny app also caps how large a pedigree it will render as a diagram, and the cap itself is style-aware: 750 animals under direct, 400 under rectilinear – because the rectilinear style’s extra waypoint nodes make an equivalent-size diagram roughly 1.9x heavier to render. Both caps are internal to the Shiny app (R/modPedigree.R); the script-callable functions demonstrated here have no such limit of their own.

Twin/Zygosity Connectors

A recorded pedigree has no way to represent that two siblings are twins, or how confident that twin call is – that information has to come from outside the pedigree, typically a curator-maintained sidecar file. readTwinRelations reads such a file (Excel or delimited text, columns id1, id2, code) with no validation of its own; checkTwinRelations then checks the result against the pedigree – both ids must exist, an MZ twin (monozygotic) or DZ twin (dizygotic) pair must already share both recorded parents, and an MZ twin pair must additionally match in recorded sex, while a UZ twin (zygosity undetermined) pair has no such precondition. The result of checkTwinRelations should be passed to makePedigreeMatingLayout’s optional twinRelations argument; the layout function does not validate it itself. The argument adds a dashed twin-connector edge between the pair, styled by zygosity code – the same connectors shown in the live app’s Diagram tab.

twinPed <- data.frame(
  id = c("F1", "F2", "S1", "S2"),
  sire = c(NA, NA, "F1", "F1"),
  dam = c(NA, NA, "F2", "F2"),
  sex = c("M", "F", "F", "M"),
  stringsAsFactors = FALSE
)
twinPed$gen <- findGeneration(twinPed$id, twinPed$sire, twinPed$dam)
twinRelationsCsv <- tempfile(fileext = ".csv")
writeLines(c("id1,id2,code", "S1,S2,UZ twin"), twinRelationsCsv)

S1 and S2 are full siblings (both recorded parents are F1/F2), so a UZ twin call between them passes validation:

twinRelations <- checkTwinRelations(twinPed, readTwinRelations(twinRelationsCsv))
twinRelations

twinDiagramData <- makePedigreeMatingLayout(twinPed, twinRelations = twinRelations)
twinDiagramData$edges

Compare the edges table above to what makePedigreeMatingLayout returns for the same pedigree without twinRelations: the union/parent -child edges are identical, but only the version above has the extra row connecting S1 directly to S2 (dashed, labeled "?" for the undetermined zygosity code) – that row exists only because a validated twinRelations table was supplied.

Kinship of Identical Twins

Of the three zygosity codes, only MZ twin changes a kinship value: a DZ twin or UZ twin row leaves kinship’s result as it is. Identical twins carry the same genes, so kinship, given a twinRelations table, sets each declared MZ twin pair’s kinship to the kinship an animal has with itself, and carries the change on to the twins’ relatives. The pair above is a female and a male and cannot be an MZ twin pair, so this small pedigree has two full sisters declared identical twins and a daughter of the first of them:

mzPed <- data.frame(
  id = c("F1", "F2", "T1", "T2", "M3", "D1"),
  sire = c(NA, NA, "F1", "F1", NA, "M3"),
  dam = c(NA, NA, "F2", "F2", NA, "T1"),
  sex = c("M", "F", "F", "F", "M", "F"),
  stringsAsFactors = FALSE
)
mzPed$gen <- findGeneration(mzPed$id, mzPed$sire, mzPed$dam)
mzTwins <- checkTwinRelations(mzPed, data.frame(
  id1 = "T1", id2 = "T2", code = "MZ twin", stringsAsFactors = FALSE
))
plainKmat <- kinship(mzPed$id, mzPed$sire, mzPed$dam, mzPed$gen)
twinKmat <- kinship(mzPed$id, mzPed$sire, mzPed$dam, mzPed$gen,
                    twinRelations = mzTwins)
plainKmat[c("T1", "T2", "D1"), c("T1", "T2", "D1")]

##      T1    T2    D1
## T1 0.50 0.250 0.250
## T2 0.25 0.500 0.125
## D1 0.25 0.125 0.500

twinKmat[c("T1", "T2", "D1"), c("T1", "T2", "D1")]

##      T1   T2   D1
## T1 0.50 0.50 0.25
## T2 0.50 0.50 0.25
## D1 0.25 0.25 0.50

The twins’ kinship rises from 0.25, that of full sisters, to 0.5, and the daughter’s kinship with her mother’s twin rises from 0.125, that of a niece and her aunt, to 0.25, the same as with her own mother. A matrix built this way is the one to pass on as kinshipMatrix, as the next section explains.

Passing a Kinship Matrix You Already Have

makePedigreeMatingLayout marks a mating unit as consanguineous when its sire and dam are related, and it works that out by computing a kinship matrix of the pedigree itself. If you already hold a kinship matrix for the same pedigree, pass it as kinshipMatrix and the layout uses it instead of making its own. The matrix needs the animal ids as its row and column names. It is then the only source of the consanguinity marks: twinRelations goes on drawing the twin connectors and nothing else, so a matrix that should reflect declared monozygotic twins has to be built with kinship’s own twinRelations argument before it is passed in.

demoKmat <- kinship(demoPed$id, demoPed$sire, demoPed$dam, demoPed$gen)
layoutFromPed <- makePedigreeMatingLayout(demoPed)
layoutFromMatrix <- makePedigreeMatingLayout(demoPed, kinshipMatrix = demoKmat)
identical(layoutFromMatrix$edges, layoutFromPed$edges)

## [1] TRUE

unnamedMatrixMessage <- tryCatch(
  makePedigreeMatingLayout(demoPed, kinshipMatrix = unname(demoKmat)),
  error = conditionMessage
)
unnamedMatrixMessage

## [1] "makePedigreeMatingLayout() requires 'kinshipMatrix' to be a matrix (or Matrix) with row and column names set to individual ids, or NULL."

The layout built from the matrix has the same edges as the one that computed its own, and a matrix without ids is refused with a message that names kinshipMatrix.

Using known genotypes

Known genotypes give the gene drop real alleles to start from in place of simulated ones. The Input tab reads them when “Pedigree(s) and genotypes in one file” is chosen and the file has two allele columns (first and second, first_name and second_name, or allele_1 and allele_2), and its QC Summary says which it read. A script gives a pedigree the same integer codes with addGenotype. It numbers each distinct allele label from 10001, in sorted order over both columns, and adds the codes as the integer columns first and second. reportGV finds those two columns by itself, and getGVGenotype shows what it finds.

The next chunk reads the package’s one-file example, a rhesus MHC pedigree in which some animals have allele labels in first_name and second_name. It hands addGenotype the rows that have labels, counts the animals with codes, runs reportGV and counts the genotyped animals at each genome uniqueness value. Each of those animals starts the gene drop from its own two alleles; the alleles of every other animal are simulated or passed down from them.

genotypeFile <- system.file("extdata", "examples",
  "obfuscated_rhesus_mhc_ped_with_genotypes.csv",
  package = "nprcgenekeepr"
)
labelPed <- qcStudbook(getPedigree(genotypeFile),
  reportErrors = FALSE, reportChanges = FALSE
)
labelColumns <- c("first_name", "second_name")
hasLabels <- !is.na(labelPed$first_name) | !is.na(labelPed$second_name)
genotypePed <- addGenotype(
  labelPed[, !(names(labelPed) %in% labelColumns)],
  labelPed[hasLabels, c("id", labelColumns)]
)
c(animals = nrow(genotypePed), with_codes = sum(!is.na(genotypePed$first)))

##    animals with_codes 
##        375         31

genotypeGV <- reportGV(genotypePed,
  guIter = 200L,
  guThresh = 3L,
  byID = TRUE,
  updateProgress = NULL
)
genotypedGu <- genotypeGV$report$gu[genotypeGV$report$id %in%
  genotypePed$id[!is.na(genotypePed$first)]]
table(genotypedGu)

## genotypedGu
##   0  50 100 
##   4  15  12

Detailed look at the Genetic Value Report object

The names of the object within the genetic value report object (trimmedGeneticValue) can be listed as shown in the next line of code.

names(trimmedGeneticValue)

##  [1] "report"          "kinship"         "gu"              "fe"             
##  [5] "fg"              "fgSE"            "neGD"            "neSexRatio"     
##  [9] "neVariance"      "maleFounders"    "femaleFounders"  "nMaleFounders"  
## [13] "nFemaleFounders" "total"

The report object (an R dataframe) can in-turn be examined.

names(trimmedGeneticValue$report) ## column names

##  [1] "id"              "sex"             "age"             "birth"          
##  [5] "exit"            "population"      "origin"          "sire"           
##  [9] "dam"             "indivMeanKin"    "zScores"         "gu"             
## [13] "guSE"            "totalOffspring"  "livingOffspring" "parentage"      
## [17] "flagged"         "value"           "rank"

nrow(trimmedGeneticValue$report) ## Number of rows

## [1] 327

The report is more conveniently used as a separate object. The next section of code rounds some of the numerical values and converts all columns to characters for display as a table where only the first 10 lines are included.

rpt <- trimmedGeneticValue[["report"]]
rpt$indivMeanKin <- round(rpt$indivMeanKin, 5L)
rpt$zScores <- round(rpt$zScores, 2L)
rpt$gu <- round(rpt$gu, 5L)
rpt <- toCharacter(rpt)
names(rpt) <- headerDisplayNames(names(rpt))
knitr::kable(rpt[1L:10L, ]) # needs more work for display purposes.

We start the next lines of code by getting a fresh copy of the genetic value report since we changed all of the numeric values to characters in the last section to print the table. These lines demonstrate one way of extracting the component objects from the genetic value object.

rpt <- trimmedGeneticValue[["report"]]
kmat <- trimmedGeneticValue[["kinship"]]
f <- trimmedGeneticValue[["total"]]
mf <- trimmedGeneticValue[["maleFounders"]]
ff <- trimmedGeneticValue[["femaleFounders"]]
nmf <- trimmedGeneticValue[["nMaleFounders"]]
nff <- trimmedGeneticValue[["nFemaleFounders"]]
fe <- trimmedGeneticValue[["fe"]]
fg <- trimmedGeneticValue[["fg"]]

It is informative to examine the distribution of genetic uniqueness, mean kinship, and z-scores (normalized mean kinship values).

Creation of the boxplot for the genetic uniqueness values is shown below.

gu <- rpt[, "gu"]
guBox <- ggplot(data.frame(gu = gu), aes(x = "", y = gu)) +
  geom_boxplot(
    color = "darkblue",
    fill = "lightblue",
    notch = TRUE,
    outlier.color = "red",
    outlier.shape = 1L
  ) +
  theme_classic() +
  geom_jitter(width = 0.2) +
  coord_flip() +
  ylab("Score") +
  ggtitle("Genetic Uniqueness")
print(guBox)

Extraction of the individual mean kinship values and their corresponding z-scores is shown in the next code chunk.

mk <- rpt[, "indivMeanKin"]
zs <- rpt[, "zScores"]

Creation of boxplots for the mean kinship and z-scores is left as an exercise.

X-Chromosome Kinship

By default kinship gives the autosomal kinship coefficient. With chrtype = "x" it gives X-chromosome kinship instead. A male carries a single X, which comes from his mother alone; a female’s X-linked kinship follows the same average-of-the-parents rule as the autosomes. Because the calculation depends on sex, sex is then required: one value per animal in the order of id, using the package’s "M" and "F" codes (any other value, NA included, counts as unknown). An animal of unknown sex that has recorded parents gets NA kinship with everyone.

xPed <- nprcgenekeepr::smallPed
xPed$gen <- findGeneration(xPed$id, xPed$sire, xPed$dam)
autosomeKmat <- kinship(xPed$id, xPed$sire, xPed$dam, xPed$gen)
xKmat <- kinship(xPed$id, xPed$sire, xPed$dam, xPed$gen,
                 chrtype = "x", sex = xPed$sex)
unique(data.frame(sex = xPed$sex, autosome = unname(diag(autosomeKmat)),
                  x = unname(diag(xKmat))))

xPairs <- data.frame(
  first = c("A", "B", "D", "E"), second = c("C", "C", "F", "F"),
  relation = c("father and son", "mother and son", "father and daughter",
               "mother and daughter"))
pairIndex <- as.matrix(xPairs[, c("first", "second")])
xPairs$autosome <- autosomeKmat[pairIndex]
xPairs$x <- xKmat[pairIndex]
xPairs

The first table is each animal’s kinship with itself: on the X a male’s is higher than his autosomal value, because he carries a single X, and a female’s is unchanged. In the second, a father and his son share no X, since a son’s X comes from his mother alone, while a mother and her son, and a father and his daughter, share more of it than they do on the autosomes.

Recording a Snapshot

createColonySnapshot turns one reportGV result into a one-row snapshot. Give it the pedigree the analysis was run on, the result, the membershipRule that says who was analyzed (here "wholePedigree", every animal in the pedigree), and the same guIter and guThresh the analysis used, because the snapshot records them and the result does not carry them. appendColonySnapshot adds the row to a history (NULL starts a new one) and returns the history in date order; it writes nothing. The date is fixed here only so that the printed values are the same each time the document is built; leave it out to date the snapshot today.

set.seed(1L)
snapshotPed <- qcPed
snapshotGV <- reportGV(snapshotPed, guIter = 50L, guThresh = 1L)
firstSnapshot <- createColonySnapshot(snapshotPed, snapshotGV,
  "wholePedigree",
  guIter = 50L, guThresh = 1L, snapshotDate = as.Date("2018-01-01")
)
snapshotHistory <- appendColonySnapshot(NULL, firstSnapshot)
snapshotHistory[, c("snapshotDate", "membershipRule", "nAnimals", "fe", "fg")]

The snapshot has 27 columns; the table shows five of them.

Rebuilding Earlier Dates

A history of one row shows no trend. backfillColonySnapshots rebuilds a snapshot for each date you name from the birth and exit dates already in the pedigree, so you need not wait for seasons of recording. It works out which animals were alive on each date, analyzes them and their ancestors with reportGV, and summarizes the result as createColonySnapshot does. Give it the full pedigree, not one already trimmed to a focal list of animals. The rows are filed under their own membership rule, "reconstructedLiving", so they form a line of their own and never mix with recorded rows. Passing the history as existingHistory keeps its rows and adds the rebuilt ones.

Before the gene drop is spent, previewBackfillSupport gives the counts each date rests on, without the analysis.

set.seed(1L)
earlierDates <- as.Date(c("2005-01-01", "2008-01-01", "2010-01-01"))
backfillPreview <- previewBackfillSupport(snapshotPed, dates = earlierDates)
backfillPreview

earlierSnapshots <- backfillColonySnapshots(snapshotPed,
  dates = earlierDates,
  guIter = 50L, guThresh = 1L, existingHistory = snapshotHistory
)
snapshotHistory <- earlierSnapshots$history
snapshotHistory[, c("snapshotDate", "membershipRule", "nAnimals", "fe", "fg")]

nLiving is the number of animals each date rests on, and pctSireKnown and pctDamKnown are the percentage of them whose sire or dam is recorded. nPlacedByOffspring and nPossiblyAlive count animals with no birth date (?backfillColonySnapshots says how each is placed). Read these counts beside the trend: where few sires are recorded the sex-ratio effective size reads low because of the missing records, not because of the colony.

A rebuilt snapshot looks back with today’s records: a parent identified later counts as known earlier, today’s genotypes are used for every date, and only animals with a birth date, or with an offspring that has one, are placed. Two close dates can differ only by the random sampling every gene drop has.

A date that cannot be analyzed does not stop the others. It is listed in the skipped element with the reason, and one warning says how many dates failed. Here 0 dates were skipped.

Saving, Reading and Checking a History

Neither appendColonySnapshot nor backfillColonySnapshots writes a file. You keep the history yourself, as the application does when you download it, with write.csv and row.names = FALSE. readSnapshotHistory reads the file back (delimited text or Excel) but does not check it, and it leaves the dates as text. checkSnapshotHistory checks the structure and the values, stops with a message that names the problem, and otherwise returns the history with the dates as dates.

historyFile <- tempfile(fileext = ".csv")
write.csv(snapshotHistory, historyFile, row.names = FALSE)
readBack <- readSnapshotHistory(historyFile)
class(readBack$snapshotDate)

## [1] "character"

checkedHistory <- checkSnapshotHistory(readBack)
class(checkedHistory$snapshotDate)

## [1] "Date"

Comparing Two Dates

calcSnapshotDeltas lists, metric by metric, the value on each of two dates and the difference between them. The two dates must belong to one membership-rule series, because rules are never compared with each other. A history that holds both recorded and rebuilt rows therefore needs to be told which series; without membershipRule the function stops and names the rules it finds.

mixedRuleMessage <- tryCatch(
  calcSnapshotDeltas(checkedHistory, from = "2005-01-01", to = "2010-01-01"),
  error = conditionMessage
)
mixedRuleMessage

## [1] "Snapshot history holds several membershipRule values (reconstructedLiving, wholePedigree); pass 'membershipRule' to pick the series to compare."

deltas <- calcSnapshotDeltas(checkedHistory,
  from = "2005-01-01", to = "2010-01-01",
  membershipRule = "reconstructedLiving"
)
deltas[deltas$metric %in% c("nAnimals", "fe", "fg", "neSexRatio"), ]

The comparabilityFlag column is NA when the two snapshots were made with the same settings. When they were not, it names each setting that differs, so a change caused by a different number of gene-drop iterations is not mistaken for a change in the colony.

Plotting the Trend

plotSnapshotTrends draws the history as one ggplot object with a panel for each metric, the date along the bottom and one colored series for each membership rule. The rebuilt series is dashed, and the caption says that the dashed lines were rebuilt and not recorded at the time. Where the history carries a standard error, as for founder genome equivalents and mean genome uniqueness, the panel shows a ribbon for the sampling uncertainty of the gene drop. A snapshot made with different settings from the previous one of its series is drawn with a different point shape and named in the caption. A trend needs at least two snapshots, and metrics chooses the panels (the default draws every metric column). ncol sets how many columns of panels to draw: three keeps these three panels side by side, and a narrow figure needs fewer, so that each panel stays wide enough for its title.

trendPlot <- plotSnapshotTrends(checkedHistory,
                                metrics = c("fe", "fg", "meanGu"), ncol = 3L)
trendPlot

Harems

The following group assignments will be forming harem groups. This is done by setting harem to . Setting iter to 100 or more will increase optimal composition of breeding groups

haremGrp <- groupAddAssign(
  kmat = trimmedGeneticValue[["kinship"]],
  ped = trimmedPed,
  candidates = candidates,
  iter = 10L,
  numGp = 6L,
  harem = TRUE
)
haremGrp$group

## [[1]]
##  [1] "SXSVEH" "G8MCV7" "99BMJW" "5EDIEE" "6X6BG9" "2Z4YLY" "IH1KPA" "0X4W26"
##  [9] "XEC0M5" "3GECJJ" "Q7U139" "6F9FB8" "JLFKV8" "B1WVCN" "BKWE4D" "C18V6I"
## [17] "F45799" "S056D5" "QRZK48" "DHNQ1W" "D33J06" "38K2SR" "TEACA3" "5EDLL7"
## [25] "M9PVG5" "Y6DB6L" "LS184H" "AR5U44" "ZPS15A" "1CZM30"
## 
## [[2]]
##  [1] "9FRCIE" "50D77I" "AW400C" "Q8U9LB" "7B9CA6" "5KWNMZ" "92UG4N" "1KJ2MG"
##  [9] "PI4VHT" "87AQLF" "YLRNIK" "DCJJYS" "B228Q6" "MPIQ4N" "WLMGS1" "Z904TJ"
## [17] "VWC5ZH" "FG0SFA" "W0GUKI" "MFKT9C" "PVY432" "W5WIRP" "2F6J3U" "1FAZ0K"
## [25] "MKY9TK" "MTCAIG" "WJXIH9"
## 
## [[3]]
##  [1] "1E8KD1" "YTJ2UL" "WK89I9" "72LYDE" "465ERA" "FL170P" "CS23RV" "DKIM6U"
##  [9] "DI4AHD" "B134XZ" "SHG3RB" "5W621W" "MX4J7G" "QW2Z3R" "RVHVTZ" "46ZHKN"
## [17] "PU7RSG" "FB5L3N" "CRPXY7" "EX5K0S" "SCFSBF" "6KWVRI" "J3F6PD" "9MG040"
## [25] "7NE2UT" "X694YR" "FJS7RQ" "I5CI33" "I8ABC7" "WNEAS6" "1QVS67" "83HQBN"
## 
## [[4]]
##  [1] "5KFB90" "AFZKBS" "QCA36T" "967Y3D" "N79QXB" "DPXEQE" "8JUUJ9" "QQMBT1"
##  [9] "QWKFBH" "MYUMMX" "D4B0RM" "DRXMW4" "PBAFJF" "QCENKM" "ZQXZYB" "1SPLS8"
## [17] "EZ2F8A" "KX0RJ3" "R5AYJK" "9P0DES" "F7I2ED" "S5H1GC" "Y0TCYX" "YFCIHJ"
## [25] "WKY2SZ" "XYRDKV" "S7IWWA" "5ERY5Z" "1SSCJC" "MH88T6" "EMV4P6" "01QRQ4"
## [33] "AIHJ8Z" "5BPBUI" "KZY6PD" "ESUIAF" "321LLB" "CLSVU6" "BCJJKN" "0XTZQ1"
## [41] "414N7M" "SH3FB7" "7ZNY75"
## 
## [[5]]
##  [1] "80F2MI" "ILVQVB" "0SGJ12" "W6MDVK" "TXZUKC" "AP1YLW" "ZH3YG1" "PYPM1W"
##  [9] "CMMUKU" "G25E3F" "D9P18Y" "GTLA8R" "AR17R5" "1CIRC9" "3DTD2N" "XFWVVX"
## [17] "GCBYDW" "30J3CQ" "KEA4QG" "E5Q33K" "S3EBGZ" "AZ3L0D" "PJ72W1" "13B1QL"
## [25] "LVYYNY" "G58RGY" "1VP3UC" "S222R3" "NK802Y" "ZATMEE" "0IIAEN" "GIIEUD"
## 
## [[6]]
##  [1] "B2CKHA" "HE0SCR" "N5QBWD" "7RA57Q" "3YJIMV" "DH9WJQ" "WI38KZ" "RJ4JPC"
##  [9] "K3TNHP" "CHK1ZX" "TYEWF1" "TQEMY6" "N4NV8B" "Q17CG3" "0HYZ23" "GAS52W"
## [17] "WTE53B" "LYSLPP" "MB6NYQ" "H2J6UA" "5IAFMK" "BS3RLE" "T3QPW5" "3SKITJ"
## [25] "1GF3GM" "NN3GDQ"
## 
## [[7]]
## [1] NA

Notice that the list has one more element than the six groups requested. The last element is not a breeding group. With harem groups it is a lone NA, which is not an animal (the next section shows what it holds when a sex ratio is set).

We can identify and list the males in each group with the following code.

sapply(haremGrp$group, function(ids) {
  ids[ids %in% trimmedPed$id[trimmedPed$sex == "M"]]
})

## [[1]]
## [1] "SXSVEH"
## 
## [[2]]
## [1] "9FRCIE"
## 
## [[3]]
## [1] "1E8KD1"
## 
## [[4]]
## [1] "5KFB90"
## 
## [[5]]
## [1] "80F2MI"
## 
## [[6]]
## [1] "B2CKHA"
## 
## [[7]]
## logical(0)

It is easy to notice that the male is listed first within each breeding group.

We can also see the number of animals and the sex ratios created in each group. Since these are harem groups the sex ratios are determined by the number of animals in the individual groups.

lines <- sapply(haremGrp$group, function(ids) {
  paste0(
    "Count: ", length(ids), " Sex Ratio: ",
    round(calculateSexRatio(ids, trimmedPed), 2)
  )
})
for (line in lines) print(line)

## [1] "Count: 30 Sex Ratio: 29"
## [1] "Count: 27 Sex Ratio: 26"
## [1] "Count: 32 Sex Ratio: 31"
## [1] "Count: 43 Sex Ratio: 42"
## [1] "Count: 32 Sex Ratio: 31"
## [1] "Count: 26 Sex Ratio: 25"
## [1] "Count: 1 Sex Ratio: Inf"

Examination of this table shows that of the 184 females 184 are included in the six groups.

Controlling Sex Ratios

The following group assignments control the sex ratio of the groups formed. This is done by setting sexRatio to the desired number of females per male (here 9, with harem left at its default of FALSE).

sexRatioGrp <- groupAddAssign(
  kmat = trimmedGeneticValue[["kinship"]],
  ped = trimmedPed,
  candidates = candidates,
  iter = 10L,
  numGp = 6L,
  sexRatio = 9.0
)
sexRatioGrp$group

## [[1]]
##  [1] "PVY432" "E3JP0C" "0HYZ23" "321LLB" "ZATMEE" "KZY6PD" "W5WIRP" "B134XZ"
##  [9] "F45799" "PJ72W1" "1SSCJC" "1GF3GM" "WK89I9" "KZM9RB" "MYUMMX" "QCA36T"
## [17] "MTCAIG" "DI4AHD" "R5AYJK" "DH9WJQ" "X694YR" "BKWE4D" "S056D5" "YLRNIK"
## [25] "KXHGRH" "Y6DB6L" "LS184H" "BS3RLE"
## 
## [[2]]
##  [1] "ESUIAF" "RY1AZM" "G58RGY" "RJ4JPC" "F7I2ED" "6F9FB8" "NK802Y" "3SKITJ"
##  [9] "EMV4P6" "QWKFBH" "6KWVRI" "ZPS15A" "M9PVG5" "WQUN84" "Q17CG3" "BCJJKN"
## [17] "TQEMY6" "CRPXY7" "VWC5ZH" "AR17R5" "YTJ2UL" "CLSVU6" "7ZNY75" "414N7M"
## [25] "A6A1M1" "I8ABC7" "AZ3L0D" "0IIAEN" "S222R3" "6X6BG9" "J3F6PD" "SCFSBF"
## [33] "7NE2UT"
## 
## [[3]]
##  [1] "1CZM30" "QRWYQZ" "CHK1ZX" "S3EBGZ" "FB5L3N" "AFZKBS" "D4B0RM" "3DTD2N"
##  [9] "C18V6I" "2Z4YLY" "01QRQ4" "1QVS67" "GCBYDW" "K7900I" "92UG4N" "MFKT9C"
## [17] "Y0TCYX" "0XTZQ1" "DCJJYS" "30J3CQ" "AW400C" "IH1KPA" "G25E3F" "GIIEUD"
## [25] "SXSVEH" "XFWVVX" "1KJ2MG" "HE0SCR" "MB6NYQ" "Q7U139"
## 
## [[4]]
##  [1] "3GECJJ" "3QHAFI" "5KWNMZ" "DRXMW4" "1FAZ0K" "FL170P" "13B1QL" "TXZUKC"
##  [9] "WI38KZ" "8JUUJ9" "MH88T6" "0SGJ12" "GTLA8R" "JSAP3H" "QCENKM" "TEACA3"
## [17] "72LYDE" "JLFKV8" "WLMGS1" "50D77I" "PU7RSG" "W6MDVK" "5EDIEE" "YFCIHJ"
## [25] "IZDV8K" "QW2Z3R" "465ERA" "38K2SR" "967Y3D" "WTE53B" "RVHVTZ" "87AQLF"
## [33] "AR5U44" "TYEWF1" "Z7NBA2"
## 
## [[5]]
##  [1] "D9P18Y" "XL658N" "LVYYNY" "9P0DES" "MX4J7G" "FJS7RQ" "2F6J3U" "ZQXZYB"
##  [9] "46ZHKN" "EZ2F8A" "KX0RJ3" "B1WVCN" "WNEAS6" "T38W6H" "9MG040" "N79QXB"
## [17] "DKIM6U" "Q8U9LB" "AIHJ8Z" "CMMUKU" "W0GUKI" "ILVQVB" "PI4VHT" "5EDLL7"
## [25] "LN1DLY" "DPXEQE" "S5H1GC" "D33J06" "SH3FB7" "99BMJW" "5ERY5Z"
## 
## [[6]]
##  [1] "7B9CA6" "ER464J" "1SPLS8" "0X4W26" "7RA57Q" "1CIRC9" "XEC0M5" "FG0SFA"
##  [9] "AP1YLW" "GAS52W" "5BPBUI" "NN3GDQ" "83HQBN" "R6HV9A" "CS23RV" "H2J6UA"
## [17] "K3TNHP" "SHG3RB" "MPIQ4N" "PBAFJF" "T3QPW5" "B228Q6" "3YJIMV" "QRZK48"
## [25] "A98D7P" "1VP3UC" "WKY2SZ" "KEA4QG"
## 
## [[7]]
##  [1] "5IAFMK" "HLQ9SY" "B2CKHA" "GDXWJ1" "TR5L57" "XC304E" "1E8KD1" "5KFB90"
##  [9] "5PW7WT" "AEP5EG" "BW10CL" "CFD12A" "CHJ9D2" "FTVE03" "G2GYST" "G8MCV7"
## [17] "IRFJ09" "LMJWTN" "N4NV8B" "N5QBWD" "Q9LWGX" "RNQU14" "09LFE4" "3MMZD4"
## [25] "55BPSE" "5XVTVH" "8IG767" "9FRCIE" "FFGPS4" "FG6L7S" "LYSLPP" "NHWTJ9"
## [33] "P7RBPI" "TBCE78" "YI16QD" "2F1IV1" "4LHK19" "59NYZE" "5IYDXN" "6KLWVC"
## [41] "80F2MI" "8TV4MT" "AZ4D19" "BTTHAJ" "CHSCFG" "DHNQ1W" "EEGLWY" "FX9E4X"
## [49] "G91ZM6" "I5CI33" "J1R2EW" "LDND6J" "MKY9TK" "MQT080" "NSIC4I" "PHB6TE"
## [57] "PYPM1W" "WHQLH5" "XX0GYV" "XYRDKV" "YHHVC7" "YP910X" "Z904TJ" "ZH3YG1"
## [65] "0V4SAC" "0X1RZ9" "3YHBC1" "55VDSQ" "5W621W" "653J82" "6MEP2C" "76DIT4"
## [73] "80KACX" "B2YJJP" "E5Q33K" "EX5K0S" "FLIZQI" "GM371F" "MEUZ85" "PA9F3J"
## [81] "QQMBT1" "S7IWWA" "TJN1AD" "WJXIH9" "WNKKW3" "XY2CK7" "XZH41H" "YDRD81"
## [89] "Z25D52" "ZDRSG0" "3P9BX6" "7D09WH" "DGZLV3" "S63QDN" "ZW2X4N"

As before, the last element of the list is not a breeding group. Here it holds the 95 candidates that were not placed in any of the six groups.

Again we can identify and list the males in each group with the following code.

sapply(sexRatioGrp$group, function(ids) {
  ids[ids %in% trimmedPed$id[trimmedPed$sex == "M"]]
})

## [[1]]
## [1] "E3JP0C" "KZM9RB" "KXHGRH"
## 
## [[2]]
## [1] "RY1AZM" "WQUN84" "A6A1M1"
## 
## [[3]]
## [1] "QRWYQZ" "K7900I" "SXSVEH"
## 
## [[4]]
## [1] "3QHAFI" "JSAP3H" "IZDV8K" "Z7NBA2"
## 
## [[5]]
## [1] "XL658N" "T38W6H" "LN1DLY"
## 
## [[6]]
## [1] "ER464J" "R6HV9A" "A98D7P"
## 
## [[7]]
##  [1] "HLQ9SY" "B2CKHA" "GDXWJ1" "TR5L57" "XC304E" "1E8KD1" "5KFB90" "5PW7WT"
##  [9] "AEP5EG" "BW10CL" "CFD12A" "CHJ9D2" "FTVE03" "G2GYST" "IRFJ09" "LMJWTN"
## [17] "Q9LWGX" "RNQU14" "09LFE4" "3MMZD4" "55BPSE" "5XVTVH" "8IG767" "9FRCIE"
## [25] "FFGPS4" "FG6L7S" "NHWTJ9" "P7RBPI" "TBCE78" "YI16QD" "2F1IV1" "4LHK19"
## [33] "59NYZE" "5IYDXN" "6KLWVC" "80F2MI" "8TV4MT" "AZ4D19" "BTTHAJ" "CHSCFG"
## [41] "EEGLWY" "FX9E4X" "G91ZM6" "J1R2EW" "LDND6J" "MQT080" "NSIC4I" "PHB6TE"
## [49] "WHQLH5" "XX0GYV" "YHHVC7" "YP910X" "0V4SAC" "0X1RZ9" "3YHBC1" "55VDSQ"
## [57] "653J82" "6MEP2C" "76DIT4" "80KACX" "B2YJJP" "FLIZQI" "GM371F" "MEUZ85"
## [65] "PA9F3J" "TJN1AD" "WNKKW3" "XY2CK7" "XZH41H" "YDRD81" "Z25D52" "ZDRSG0"
## [73] "3P9BX6" "7D09WH" "DGZLV3" "S63QDN" "ZW2X4N"

We can also see the number of animals and the sex ratios created in each group.

lines <- sapply(sexRatioGrp$group, function(ids) {
  paste0(
    "Count: ", length(ids), " Sex Ratio: ",
    round(calculateSexRatio(ids, trimmedPed), 2L)
  )
})
for (line in lines) print(line)

## [1] "Count: 28 Sex Ratio: 8.33"
## [1] "Count: 33 Sex Ratio: 10"
## [1] "Count: 30 Sex Ratio: 9"
## [1] "Count: 35 Sex Ratio: 7.75"
## [1] "Count: 31 Sex Ratio: 9.33"
## [1] "Count: 28 Sex Ratio: 8.33"
## [1] "Count: 95 Sex Ratio: 0.23"

Examination of this table shows that of the 184 females 166 are included in the six groups.

Ancestry Rules for Mate Pairs

Some centers do not want certain ancestry combinations paired (for example, Indian-origin with Chinese-origin rhesus macaques). reportMatePairs accepts the same center-configurable rules table the Breeding Groups ancestry guardrails use, through ancestryRules. A rule matches a pair whichever animal is the sire and whichever the dam. A block rule moves the matching pair out of pairs and into excluded with the reason "ancestry rule" – here “block” means “not listed as eligible”, not “never placed together”. A flag rule keeps the pair in pairs and annotates it. The pedigree must carry an ancestry column, which the example pedigree below does.

ancPed <- qcStudbook(
  read.csv(system.file("extdata", "examples", "example_ancestry_pedigree.csv",
                       package = "nprcgenekeepr"),
           stringsAsFactors = FALSE, na.strings = c("", "NA")),
  minSireAge = 2, minDamAge = 2, reportChanges = FALSE, reportErrors = FALSE)
ancPed$gen <- findGeneration(ancPed$id, ancPed$sire, ancPed$dam)
ancKmat <- kinship(ancPed$id, ancPed$sire, ancPed$dam, ancPed$gen)
ancestryRules <- checkAncestryRules(readAncestryRules(
  system.file("extdata", "examples", "example_ancestry_rules.csv",
              package = "nprcgenekeepr")))
ancestryRules

ancMatePairs <- reportMatePairs(ancPed, ancKmat, ancestryRules = ancestryRules)
names(ancMatePairs)

## [1] "pairs"            "excluded"         "ancestryCoverage"

With rules supplied, pairs gains three columns: ancestryRule (the matched rule as a sorted "LEVEL-LEVEL" string), ancestrySeverity ("flag" or "block") and ancestryStatus ("violation", or "overridden" for an overridden block rule). All three are NA for a pair no rule matches. Here the only rows with a rule are the three flagged pairs; the five blocked pairs are in excluded.

ancMatePairs$pairs[!is.na(ancMatePairs$pairs$ancestryRule),
                   c("sireId", "damId", "ancestryRule", "ancestrySeverity",
                     "ancestryStatus")]

ancMatePairs$excluded

The third element, ancestryCoverage, lists every standardized ancestry level with the number of distinct animals at that level (n) and whether any rule names it (covered). A level with animals but covered = FALSE – here JAPANESE – is one no rule speaks to, which is worth a second look before relying on the rules as complete.

ancMatePairs$ancestryCoverage

A curator can override a block rule for one call with overriddenRules, a data frame with ancestry1 and ancestry2 columns (unordered, case-insensitive). The blocked pairs then return to pairs, marked "overridden", rather than silently reappearing.

overriddenMatePairs <- reportMatePairs(
  ancPed, ancKmat, ancestryRules = ancestryRules,
  overriddenRules = data.frame(ancestry1 = "CHINESE", ancestry2 = "INDIAN"))
table(overriddenMatePairs$pairs$ancestryStatus, useNA = "ifany")

## 
## overridden  violation       <NA> 
##          3          3         17

table(overriddenMatePairs$excluded$reason)

## 
## ancestry rule 
##             2

Naming a rule that is not in ancestryRules, a flag rule, or the same rule twice is an error, as is giving overriddenRules without ancestryRules. A valid rules table with zero rules is accepted: the result has the ancestry columns, with nothing moved.

Ancestry Rules for Formed Groups

reportMatePairs applies the rules to pairs. reportAncestryViolations applies the same rules to groups that already exist, for example the group element returned by groupAddAssign without its final element, which holds the animals left unused. Give it a list of animal-id vectors, one per group, the pedigree with its ancestry column, and the rules table. It reports every pair inside a group that matches a rule, block and flag rules alike, and a coverage table: for each standardized ancestry level, how many of the grouped animals sit at that level (n) and whether any rule names it (covered). An id that is not in the pedigree stops the function.

formedGroups <- list(c("C1", "I1", "J1"), c("I2", "U1", "H1"))
ancViolations <- reportAncestryViolations(formedGroups, ancPed, ancestryRules)
ancViolations$violations

ancViolations$coverage

Each row of violations gives the group (its position in the list), the two animals and their levels, the matched rule, its severity and its status. The first group also holds a Japanese animal. No rule names that level (covered is FALSE), so that animal appears in no violation; the coverage table is there so that this silence is visible.

A rule named in overriddenRules (a data frame that gives the two ancestry levels of the rule, in either order and in any case) still reports its pairs, with the status "overridden" in place of "violation", so nothing is dropped silently. A row that matches no rule in the rules table is an error, so a typo cannot override nothing.

ancViolationsOverridden <- reportAncestryViolations(
  formedGroups, ancPed, ancestryRules,
  overriddenRules = data.frame(ancestry1 = "CHINESE", ancestry2 = "INDIAN"))
ancViolationsOverridden$violations

Preparing a Marker Genotype File

A marker genotype file is long-format: one row per id x locus, with columns id, locus, allele1, and allele2. (This is a different shape from the single-locus first_name/second_name genotype format used earlier by addGenotype/geneDrop – the two are unrelated.) You would normally read this from a CSV file with read.csv, the same way the pedigree file was read in at the start of this tutorial; here we build a small example directly, for a parent (P), her offspring (C), and an unrelated founder (U), genotyped at 10 biallelic loci.

markerGenotype <- data.frame(
  id = c(rep("P", 10L), rep("C", 10L), rep("U", 10L)),
  locus = rep(paste0("L", 1L:10L), 3L),
  allele1 = c(
    "A", "A", "A", "B", "A", "A", "B", "A", "A", "A",
    "A", "A", "B", "A", "A", "A", "B", "A", "A", "A",
    "A", "B", "A", "A", "A", "A", "A", "B", "A", "B"
  ),
  allele2 = c(
    "A", "B", "B", "B", "A", "B", "B", "B", "A", "B",
    "B", "A", "B", "B", "A", "B", "B", "B", "B", "A",
    "B", "B", "A", "B", "B", "A", "B", "B", "B", "B"
  ),
  stringsAsFactors = FALSE
)

checkMarkerGenotypeFile validates the column shape, rejects any locus with more than two distinct alleles (the estimators below all require biallelic markers), and stops on duplicate id x locus rows. buildMarkerGenotypeMatrix then pivots the checked long-format table into the wide id x locus matrix the rest of this section’s functions consume, one cell per individual/locus combination (e.g. "A/B"), or NA where that individual has no genotype call at that locus.

markerGenotype <- checkMarkerGenotypeFile(markerGenotype)
genotypeMatrix <- buildMarkerGenotypeMatrix(markerGenotype)
genotypeMatrix

##   L1    L2    L3    L4    L5    L6    L7    L8    L9    L10  
## P "A/A" "A/B" "A/B" "B/B" "A/A" "A/B" "B/B" "A/B" "A/A" "A/B"
## C "A/B" "A/A" "B/B" "A/B" "A/A" "A/B" "B/B" "A/B" "A/B" "A/A"
## U "A/B" "B/B" "A/A" "A/B" "A/B" "A/A" "A/B" "B/B" "A/B" "B/B"

Marker-Based Kinship

markerKinship estimates pairwise kinship directly from genotypeMatrix, using the “KING-robust” estimator (Manichaikul et al. 2010) – the same estimator implemented by KING, PLINK2, and SNPRelate::snpgdsIBDKING. It never looks at recorded parentage, so it is useful as an independent cross-check on the pedigree-based kinship computed by kinship/reportGV above.

markerKmat <- markerKinship(genotypeMatrix)
markerKmat

##     P    C    U
## P 0.5  0.2  0.0
## C 0.2  0.5 -0.3
## U 0.0 -0.3  0.5

P and C are the true parent/offspring pair in this example, and their marker-based kinship (0.2) is clearly higher than either pair involving the unrelated founder U. The KING-robust estimator is not bounded below by zero – C and U’s negative estimate here is informative (more divergent ancestry than the reference panel), not an error.

Heterozygosity Diagnostic

markerObservedHeterozygosity computes, per animal, the fraction of its genotyped loci at which it is heterozygous (Ho). markerExpectedHeterozygosity computes Nei’s (1973) gene diversity (He) at each locus from the population’s allele frequencies, plus the unweighted mean across loci as a population-wide summary – the standard observed-vs-expected heterozygosity diagnostic. Both take the same genotypeMatrix shape as markerKinship. This example uses three animals (X, Y, Z) genotyped at 4 loci, with Y missing a call at L4.

hetGenotype <- data.frame(
  id = c(rep("X", 4L), rep("Y", 3L), rep("Z", 4L)),
  locus = c(paste0("L", 1L:4L), paste0("L", c(1L, 2L, 3L)), paste0("L", 1L:4L)),
  allele1 = c("A", "A", "A", "A", "A", "A", "B", "B", "A", "A", "A"),
  allele2 = c("A", "B", "B", "B", "B", "A", "B", "B", "A", "B", "A"),
  stringsAsFactors = FALSE
)
hetGenotype <- checkMarkerGenotypeFile(hetGenotype)
hetMatrix <- buildMarkerGenotypeMatrix(hetGenotype)

markerObservedHeterozygosity(hetMatrix)

##         X         Y         Z 
## 0.7500000 0.3333333 0.2500000

markerExpectedHeterozygosity(hetMatrix)

## $perLocus
##        L1        L2        L3        L4 
## 0.5000000 0.2777778 0.4444444 0.3750000 
## 
## $meanHe
## [1] 0.3993056

Note Y’s observed heterozygosity is computed over its own 3 genotyped loci (1/3), not over all 4 – a missing call at one locus does not lower an animal’s own Ho, since the denominator is always that animal’s own non-missing loci, not the full panel.

Parentage Verification (Mendelian Exclusion)

markerParentageExclusion cross-references a pedigree’s recorded dam/sire against each animal’s marker genotype and flags a recorded parent whose genotype evidence contradicts simple Mendelian inheritance – directly useful for catching the kind of dam/sire misidentification a paper pedigree alone cannot reveal. A locus counts as a conflict only when the animal and the candidate parent are each homozygous for a different allele (“opposite homozygotes”); by default (maxExclusions = 2) a recorded parent is only flagged once 3 or more such conflicts accumulate, since a single mismatching locus can arise from ordinary genotyping error even for a true parent.

We reuse the P/C/U genotypeMatrix from above, with a pedigree in which C’s recorded dam (P) is correct, but C’s recorded sire has been (deliberately, for this example) misrecorded as the unrelated U.

pedigree <- data.frame(id = c("P", "C", "U"), sire = c(NA, "U", NA),
                        dam = c(NA, "P", NA), stringsAsFactors = FALSE)

markerParentageExclusion(genotypeMatrix, pedigree)

C’s true dam (P) has zero exclusions and is not flagged; C’s falsely recorded sire (U) has 3 exclusions, exceeding the default maxExclusions = 2, and is correctly flagged.

Candidate-Parent Likelihood Ranking

Once markerParentageExclusion has flagged a recorded parent as Mendelian-inconsistent with its offspring’s genotype, markerParentageLikelihood picks up the investigation: it ranks candidate replacement parents using a CERVUS-style multilocus likelihood-ratio (LOD) score (Meagher & Thompson 1986; Marshall, Slate, Kruuk & Pemberton 1998), the same approach independently validated in a real captive macaque colony by de Groot et al. (2025). It is deliberately report-only – it never rewrites pedigree\(sire_/_pedigree\)dam; a flagged parent is a curator decision, not something to silently overwrite. Candidates default to getPotentialParents’s own demographically -eligible list (breeding age, gestation window, proven-breeder preference) for the flagged offspring/role, or can be supplied explicitly, as below.

likelihoodGenotype <- data.frame(
  id = c("O", "O", "D", "D", "C1", "C1", "C2", "C2"),
  locus = c("L1", "L2", "L1", "L2", "L1", "L2", "L1", "L2"),
  allele1 = c("A", "A", "A", "A", "A", "A", "B", "B"),
  allele2 = c("A", "B", "B", "A", "A", "B", "B", "B"),
  stringsAsFactors = FALSE
)
likelihoodMatrix <- buildMarkerGenotypeMatrix(likelihoodGenotype)
likelihoodPed <- data.frame(id = "O", sire = "C2", dam = "D",
                             stringsAsFactors = FALSE)

O is the offspring, recorded dam D, recorded sire C2; C1 is an unrecorded candidate. markerParentageExclusion first shows why C2 is worth a second look:

markerParentageExclusion(likelihoodMatrix, likelihoodPed)

C2’s single exclusion (of 2 loci) falls within the default maxExclusions = 2 tolerance, so it is not flagged here – but 1 conflicting locus is still worth ranking against alternatives:

markerParentageLikelihood(likelihoodMatrix, likelihoodPed, id = "O",
                           role = "sire", candidates = c("C1", "C2"),
                           minLoci = 1L)

C1 – genetically identical to O at both loci – decisively outranks the recorded sire C2 (LOD 0.47 vs -Inf). Notice C2’s excluded column is FALSE (the same tolerant test as above) even though its LOD is -Inf: a single opposite-homozygote locus (C2 is B/B at L1, where O is A/A) is a zero-tolerance Mendelian impossibility under the likelihood model, with no allowance for ordinary genotyping error the way the exclusion count is. delta for C1 is Inf, not a large finite number – the gap down to a -Inf-ranked neighbor is reported as infinite, by design.

Validating a Cross-Center Mapping

When an animal transfers between centers that use independent id namespaces, the receiving center often has no way to know the animal’s real parents and records it as an artificial founder – losing its actual lineage even though the originating center’s records still have it. resolveCrossCenterIds (next) fixes this by merging two centers’ pedigrees given a curator-confirmed id mapping, but it stop()s on the first problem it finds in that mapping. checkCrossCenterMapping shares the same validation logic but never stops early – it returns every problem found as a row in a data.frame, so a curator sees everything wrong with a mapping in one call rather than a stop-fix-rerun cycle, one problem at a time.

In this example, T1 (in Center A’s pedigree, pedA) and X9 (in Center B’s pedigree, pedB) are the same physical animal – Center B recorded it as a founder (both parents NA) because it never knew T1’s real parents, which Center A does have.

pedA <- data.frame(
  id = c("P1", "P2", "T1"), sire = c(NA, NA, "P1"), dam = c(NA, NA, "P2"),
  stringsAsFactors = FALSE
)
pedB <- data.frame(
  id = c("X9", "O1"), sire = c(NA, "X9"), dam = c(NA, NA),
  stringsAsFactors = FALSE
)
mapping <- data.frame(idA = "T1", idB = "X9", stringsAsFactors = FALSE)

checkCrossCenterMapping(pedA, pedB, mapping)

Zero rows means the mapping is clean – resolveCrossCenterIds can be called on the same three inputs without error, as it is below. A mapping with real problems (here, idB "X9" mapped twice, and an idA value "NOPE" that doesn’t exist in pedA) instead reports both problems at once, rather than only the first:

badMapping <- data.frame(idA = c("T1", "NOPE"), idB = c("X9", "X9"),
                          stringsAsFactors = FALSE)
checkCrossCenterMapping(pedA, pedB, badMapping)

One subtlety: a genuine undeclared id collision between the two centers (an id present in both pedA and pedB but never named in mapping) is only ever reported once the existence and uniqueness problems above are already clean – a collision check on a mapping that doesn’t resolve to real pedigree rows would be meaningless. The same second check also reports a pedB animal whose id mapping already gives to a different pedA animal (a "mappedIdClash" row), and a mapped pair whose recorded parents conflict. A clean result from this function is what licenses trusting the absence of a collision, not merely the absence of a "collision" row in a dirty one.

Cross-Center Identity Linking

With the mapping validated above, resolveCrossCenterIds merges the two centers’ pedigrees into one, preferring whichever side actually recorded a parent.

resolveCrossCenterIds(pedA, pedB, mapping)

The merged pedigree has a single row for T1 with its real parents (P1, P2) intact, and O1 – recorded at Center B as X9’s offspring – now correctly points to T1 as its sire.

Cross-Center Differentiation (Fst)

markerFst takes a different, population-level view: rather than linking individual animals across centers, it estimates Hudson’s Fst (Hudson, Slatkin & Maddison 1992, as given in closed form by Bhatia et al. 2013) – how differentiated two centers’ colonies are in allele frequency, at each locus genotyped by both, plus a single pooled summary across loci. The pooled value is a ratio of summed numerators and denominators, not a mean of the per-locus ratios, which Bhatia et al. show is materially biased.

centerAGenotype <- data.frame(
  id = c(rep("CA1", 2L), rep("CA2", 2L), rep("CA3", 2L), rep("CA4", 2L)),
  locus = rep(c("L1", "L2"), 4L),
  allele1 = c("A", "A", "A", "A", "A", "A", "B", "A"),
  allele2 = c("A", "A", "A", "B", "B", "B", "B", "A"),
  stringsAsFactors = FALSE
)
centerBGenotype <- data.frame(
  id = c(rep("CB1", 2L), rep("CB2", 2L), rep("CB3", 2L),
         rep("CB4", 2L), rep("CB5", 2L), rep("CB6", 2L)),
  locus = rep(c("L1", "L2"), 6L),
  allele1 = c("A", "B", "B", "A", "A", "B", "B", "B", "A", "B", "B", "A"),
  allele2 = c("B", "B", "B", "B", "B", "B", "B", "B", "A", "B", "B", "B"),
  stringsAsFactors = FALSE
)
matrixA <- buildMarkerGenotypeMatrix(checkMarkerGenotypeFile(centerAGenotype))
matrixB <- buildMarkerGenotypeMatrix(checkMarkerGenotypeFile(centerBGenotype))

markerFst(matrixA, matrixB)

## $perLocus
##         L1         L2 
## 0.05794206 0.45129870 
## 
## $pooledFst
## [1] 0.2749664

L2 shows substantially more differentiation between the two centers than L1 does in this example; the pooled value summarizes both loci into one number for a quick between-center comparison.

Multiallelic Marker Panels and Locus Metadata

Everything above assumed a biallelic panel – checkMarkerGenotypeFile rejects any locus with more than two distinct alleles, since the KING-robust kinship estimator requires it. Real microsatellite/STR panels are routinely multiallelic, though, so the package provides a second, more permissive validator for the linkage-aware metrics below: checkLinkageMarkerGenotypeFile keeps the same structural checks (column shape, no duplicate id x locus rows) but drops the biallelic restriction.

strGenotype <- data.frame(
  id = c("W", "X", "Y", "Z"),
  locus = c("L1", "L1", "L1", "L1"),
  allele1 = c("A", "A", "A", "A"),
  allele2 = c("B", "C", "A", "D"),
  stringsAsFactors = FALSE
)

This locus has 4 distinct alleles (A, B, C, D) across the four animals – checkLinkageMarkerGenotypeFile accepts it:

checkLinkageMarkerGenotypeFile(strGenotype)

The same data fed to checkMarkerGenotypeFile instead fails, exactly as that function is designed to:

checkMarkerGenotypeFile(strGenotype)

## Error in `checkMarkerGenotypeFile()`:
## ! Marker genotype file has one or more loci with more than two distinct alleles (the KING-robust estimator requires biallelic markers): L1.

The metrics below also need per-locus genomic position information, which real curated panels rarely have complete records for. checkLocusMetadata validates a locus/chrom/pos (optionally cM) sidecar table and classifies each locus’s own coverage as "full" (chromosome and position both known), "partial" (one of the two), or "none" (neither) – downstream functions then decide for themselves which loci they can use, rather than requiring complete metadata up front.

locusMetadata <- data.frame(
  locus = c("L1", "L2", "L3"),
  chrom = c("1", "1", NA),
  pos   = c(1000000, NA, NA),
  stringsAsFactors = FALSE
)
checkLocusMetadata(locusMetadata)

L1 has both chromosome and position and is "full"; L2 has only a chromosome ("partial"); L3 has neither ("none"). computeGenomicROH (the F_ROH tab’s estimator) uses only "full"-coverage loci and drops the rest with a named warning; markerLdBlock below needs only chrom and tolerates a missing pos/cM entirely.

Realized Relatedness Variance

Pedigree-expected relatedness (kinship doubled) is only an average – actual IBD sharing between two relatives varies around that expectation because of Mendelian sampling and linkage. markerRealizedRelatednessVariance computes a closed-form estimate of that variance (Hill & Weir 2011) for Parent-Offspring, Full-Sibling, and Half-Sibling pairs, and lists only those pairs. Despite the “marker” name shared with the rest of this section, it needs no genotype data at all – only a kinship matrix, the pedigree, and the species’ chromosome count and total autosomal genetic map length (in Morgans).

smallPed <- nprcgenekeepr::smallPed
smallPed$gen <- findGeneration(smallPed$id, smallPed$sire, smallPed$dam)
smallPedKmat <- kinship(smallPed$id, smallPed$sire, smallPed$dam, smallPed$gen,
                         sparse = FALSE)
## Rhesus macaque autosome count/map length are used here only as an
## example scale -- supply values appropriate to your own species.
rrv <- markerRealizedRelatednessVariance(smallPedKmat, smallPed, nChr = 20L,
                                          mapLength = 28)
rrv

Parent-offspring pairs have varR of exactly zero, not merely small – biologically, exactly one allele at every locus is always IBD between a parent and its offspring, so there is no variance to estimate around the R = 0.5 expectation. Full-sibling variance is roughly double half -sibling variance in this example (two independently segregating parents vs. one), a useful sanity check when interpreting the estimate for your own pedigree.

Linkage-Disequilibrium Blocks

markerLdBlock takes a different, more classical view of marker relationships: pairwise linkage disequilibrium (LD) between loci sharing a chromosome, using Hedrick’s (1987) frequency-weighted D’ and a multiallelic r² generalization, from phase frequencies estimated by an EM algorithm (generalizing Excoffier & Slatkin 1995) since this package’s genotype matrix is unphased. It is explicitly the secondary, exploratory statistic in a pair with markerRealizedRelatednessVariance above – classical LD theory assumes random mating, which a pedigreed colony violates, so every row of its output carries a fixed caveat saying so.

ldGenotype <- checkLinkageMarkerGenotypeFile(data.frame(
  id = rep(c("W", "X", "Y", "Z"), 2),
  locus = rep(c("L1", "L2"), each = 4),
  allele1 = c("A", "A", "A", "A", "M", "M", "N", "N"),
  allele2 = c("A", "B", "A", "B", "M", "N", "M", "N"),
  stringsAsFactors = FALSE
))
ldGenotypeMatrix <- buildMarkerGenotypeMatrix(ldGenotype)
ldLocusMetadata <- checkLocusMetadata(data.frame(
  locus = c("L1", "L2"), chrom = c("1", "1"), pos = c(NA, NA),
  stringsAsFactors = FALSE
))

markerLdBlock(ldGenotypeMatrix, ldLocusMetadata)

D’ = 1 (the maximum possible) in this toy example – only X is a double heterozygote at both loci (the only source of phase ambiguity here), and the EM estimator resolves it to a fully-consistent phase assignment. Note r² (1/3) does not move in lockstep with D’ here – the two statistics capture different aspects of the same association, and needn’t agree.

De-identifying LD-Block Results

markerLdBlock’s own output is a locus-pair-level population statistic with no per-individual ids in it – except idsUsed, populated only when founderIds restricts the computation to a named subset. That is the one place a real id can leak, and obfuscateLdBlocks closes it: it remaps idsUsed through the same alias vector obfuscatePed’s map output produces, mirroring obfuscateTwinRelations’s existing pattern from elsewhere in the package’s de-identification family. A row whose idsUsed contains an id absent from the map stop()s loudly rather than silently leaking it.

ldBlockResult <- data.frame(
  locus1 = "L1", locus2 = "L2", chrom = "1", Dprime = 0.5, r2 = 0.3,
  nUsed = 2L, idsUsed = "F1,F2", caveat = "Descriptive statistic only.",
  stringsAsFactors = FALSE
)
obfuscatedTwinPed <- obfuscatePed(twinPed[, c("id", "sire", "dam", "sex", "gen")],
                                   map = TRUE)

obfuscateLdBlocks(ldBlockResult, obfuscatedTwinPed$map)

The real ids F1/F2 (reusing the twin-pedigree ids from the Pedigree Diagram section above) no longer appear in idsUsed – every other column (the statistics themselves) is untouched. The exact alias strings are randomly generated on every call and are not reproduced here; what matters is that no real id from the input survives into the de-identified table.

Aliasing Ids with a Known Placeholder Status

obfuscateId gives each id a random alias. By default it reads an id’s shape to decide whether it stands for an unknown parent: an id like "U0007" gets a placeholder-shaped alias (the placeholder prefix plus random characters) and every other id gets a plain one. That guess is wrong for a real animal that happens to be named "U1234". qcStudbook records the truth in a placeholder column, and obfuscateId’s placeholder argument takes that information, one logical value per id: TRUE for a placeholder and FALSE for a real animal. NA (or leaving the argument out) reads that id by its shape.

The aliases are random; the seed is fixed here only so the printed values are reproducible.

ids <- c("U1234", "A100", "U0007", "B200")
isPlaceholder <- c(FALSE, FALSE, TRUE, FALSE)
set.seed(1)
obfuscateId(ids, size = 6L)

##    U1234     A100    U0007     B200 
## "UJNV6H" "68YXCH" "UGZN1S" "09N27H"

set.seed(1)
obfuscateId(ids, size = 6L, placeholder = isPlaceholder)

##    U1234     A100    U0007     B200 
## "JNV6H6" "8YXCHG" "UZN1S0" "9N27HX"

The first call reads U1234 by its shape and gives it a placeholder-shaped alias, so the real animal looks like an unknown parent after de-identification. The second call marks it real, so it gets a plain alias (random, so it could begin with U by chance; with this seed it does not), while U0007, marked TRUE, keeps a placeholder-shaped one. The argument must be logical and the same length as id, or the call stops.