--- title: "Interactive Use of nprcgenekeepr" subtitle: An R Package for the Genetic Management of Colonies output: rmarkdown::html_vignette: df_print: kable fig_width: 7 fig_height: 5 pdf_document: default vignette: > %\VignetteEngine{knitr::rmarkdown_notangle} %\VignetteIndexEntry{Interactive Use of nprcgenekeepr} %\usepackage[UTF-8]{inputenc} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, include = TRUE, results = "markup", cache = FALSE ) pdf.options(useDingbats = TRUE) start_time <- proc.time() ``` ```{r, include = FALSE} library(stringi) # library(rmsutilityr) suppressMessages(library(ggplot2)) library(knitr) library(nprcgenekeepr) set_seed(1L) ``` ## Introduction This tutorial demonstrates the major functions used within the Shiny application provided by the __nprcgenekeepr__ package and provides sufficient insight into those functions that they may be used independently. This tutorial is primarily directed toward someone with experience using R who wants to better understand how the Shiny application works or to perform some actions not directly supported by the Shiny application. Please provide any comments, questions, or bug reports through the GitHub issue tracker at [https://github.com/rmsharp/nprcgenekeepr/issues]( https://github.com/rmsharp/nprcgenekeepr/issues). ## Installation and Help You can install **nprcgenekeepr** from CRAN with the following code. ```{r gh-installation, eval = FALSE, echo=TRUE} install.packages("nprcgenekeepr") ## Use the following code to get the development version from GitHub # install.packages("devtools") # devtools::install_github("rmsharp/nprcgenekeepr") ``` All missing dependencies should be automatically installed. You can get help from the __R__ console with ```{r gh-help, eval = FALSE, echo=TRUE} ?nprcgenekeepr ``` ## Reading in a Pedigree A pedigrees can be imported using either Excel worksheets or text files that contain all of the pedigree information or using either Excel worksheets or text files that contain a list of focal animals with the remainder of the pedigree information is pulled in through the LabKey API. This tutorial will use a pedigree file that can be created using the __makeExamplePedigreeFile__ function as shown below. The function __makeExamplePedigreeFile__ both saves a file and returns the full path name to the saved file, which we are saving into the variable _pedigreeFile_. Note: by default the file is written to the session's temporary directory (`tempdir()`); pass the _file_ argument to choose a different location. ```{r make-example-pedigree, eval = FALSE, include = TRUE} library(nprcgenekeepr) pedigreeFile <- makeExamplePedigreeFile() ``` This writes _examplePedigree.csv_ to the session's temporary directory unless you supply a _file_ argument. ```{r get-file-name, eval = TRUE, echo = FALSE, include = FALSE} pedigreeFile <- "../inst/extdata/examples/ExamplePedigree.csv" ``` You use the file name provided by the __makeExamplePedigreeFile__ function to tell __read.table__ what file to read. ```{r read-in-example-pedigree} breederPedCsv <- read.table(pedigreeFile, sep = ",", header = TRUE, stringsAsFactors = FALSE ) breederPedCsv$fromCenter <- "TRUE" breederPedCsv$fromCenter[ suppressWarnings(sample( which(is.na(breederPedCsv$sire) & is.na(breederPedCsv$dam)), round(0.8 * length(which( is.na(breederPedCsv$sire) & is.na(breederPedCsv$dam) ))) )) ] <- "FALSE" ``` Note the number of rows read. Each row represents an individual within the pedigree. ```{r row-count} nrow(breederPedCsv) ``` The next step is to put the information read from the file into a pedigree object. This is done with the __qcStudbook__ function, which examines the file contents and tests for common pedigree errors. You can see the errors that can be detected by __qcStudbook__ by returning the empty error list with __getEmptyErrorLst()__. We are not showing the output of the function call now because later in this tutorial we will explore errors in more depth. __qcStudbook__ can take five arguments _sb_, _minSireAge_ and _minDamAge_ (in years), _reportChanges_, and _reportErrors_. However, all but _sb_ have default values and only the _sb_ argument is required. It is prudent to ensure that parents are at least of breeding age, which is species and sex specific. Leaving _minSireAge_ and _minDamAge_ blank (their default, `NULL`) looks up the species-specific breeding age for each sex; a number overrides that sex's floor. I have set both to 2 years here.[^1] [^1]: Setting _minDamAge_ to 3.5 and above will cause an error along with the creation of a file _lowParentAge.csv_ in the session's temporary directory (`tempdir()`; the error message prints the path) that will list the parents with low age at the birth of an offspring: two dams in this pedigree were about 3.3 years old at a birth. The sires are all older, so raising _minSireAge_ alone to 3.5 does not trigger the error -- exactly the sex-specific control the two floors provide. ```{r form-pedigree-object} breederPed <- qcStudbook(breederPedCsv, minSireAge = 2L, minDamAge = 2L) ``` If __qcStudbook__ reports an error, change the call by adding the __reportErrors__ argument set to __TRUE__ and examine the returned object. More on this is presented in the __Pedigree Errors__ section. ## Identifying Focal Animals You may want to focus your work on a _focal_ group of animals. This can be done by reading in a list of animal IDs that make up the _focal_ group and use that list to update the pedigree. Alternatively you can created a list of animal IDs based on criteria you have selected. For example, to select living animals at the facility with at least one parent, the following can be used. ```{r select-focal-animals-from-pedigree, results='asis'} focalAnimals <- breederPed$id[!(is.na(breederPed$sire) & is.na(breederPed$dam)) & is.na(breederPed$exit)] print(stri_c( "There are ", length(focalAnimals), " animals in the vector _focalAnimals_." )) ``` As can be seen, these animals have at least one parent and have not left the facility. ```{r short-pedigree-list-of-focal-animals} breederPed[breederPed$id %in% focalAnimals, c("id", "sire", "dam", "exit")][1L:10L, ] ``` We indicate that these are the animals of interest by using the __setPopulation__ function. This function simply sets a column named _population_[^2] to the logical value of __TRUE__ if the row represents an animal in the list and __FALSE__ otherwise. [^2]: The _population_ column is created and added to the pedigree object if it does not already exist. The first line of code below sets the _population_ column and the second counts the number of rows where the value was set to __TRUE__. ```{r set-population-flag} breederPed <- setPopulation(ped = breederPed, ids = focalAnimals) nrow(breederPed[breederPed$population, ]) ``` The IDs used to populate the _population_ flag can be used to trim the pedigree so that it contains only those individuals who are in the ID list or are ancestors of those individuals. ```{r trim-pedigree} trimmedPed <- trimPedigree(focalAnimals, breederPed) nrow(breederPed) nrow(trimmedPed) ``` The __trimPedigree__ function has the ability to remove those ancestors that do not contribute genetic information. Uninformative founders are those individuals who are parents of only one individual and who have no parental information. (_Currently genotypic information is ignored by __trimPedigree___). ```{r removed-uninformatiive-animals} trimmedPedInformative <- trimPedigree(focalAnimals, breederPed, removeUninformative = TRUE ) nrow(trimmedPedInformative) ``` We can find all of the animals that are in the trimmed pedigree but are not focal animals. ```{r get-animals-added-to-focal-animals} nonfocalInTrimmedPed <- trimmedPed$id[!trimmedPed$id %in% focalAnimals] length(nonfocalInTrimmedPed) ``` We can see which of these `r length(nonfocalInTrimmedPed)` are and are not parents. We will first make sure we have all of the parents by getting our list of parents from the entire pedigree. We then demonstrate that they are all in the trimmed pedigree. ```{r find-focal-parents} allFocalParents <- c( breederPed$sire[breederPed$id %in% focalAnimals], breederPed$dam[breederPed$id %in% focalAnimals] ) trimmedFocalParents <- c( trimmedPed$sire[trimmedPed$id %in% focalAnimals], trimmedPed$dam[trimmedPed$id %in% focalAnimals] ) all.equal(allFocalParents, trimmedFocalParents) # Are the IDs the same? ``` However, not all of the animals in the trimmed pedigree are either the focal animals or their parents. They are more distant ancestors as we will show. ```{r find-non-focal-non-focal-parent} notFocalNotParent <- trimmedPed$id[!trimmedPed$id %in% c(focalAnimals, allFocalParents)] length(notFocalNotParent) ``` ```{r find-grandparents, echo = FALSE, include=FALSE, eval = TRUE} allFocalGrandParents <- c( breederPed$sire[breederPed$id %in% allFocalParents], breederPed$dam[breederPed$id %in% allFocalParents] ) ## Not all parents are known so the unknown individuals (NA) are removed. allFocalGrandParents <- allFocalGrandParents[!is.na(allFocalGrandParents)] trimmedFocalGrandParents <- c( trimmedPed$sire[trimmedPed$id %in% allFocalParents], trimmedPed$dam[trimmedPed$id %in% allFocalParents] ) trimmedFocalGrandParents <- trimmedFocalGrandParents[!is.na(trimmedFocalGrandParents)] all.equal(allFocalGrandParents, trimmedFocalGrandParents) ``` Since the trimming process is supposed to retain the focal animals and their ancestors, we will leave it as an exercise for you to demonstrate that at least some of the remaining animals are grandparents of the focal animals. _Hint: there are `r length(trimmedFocalGrandParents)` grandparents of the focal animals in the trimmed pedigree_. ```{r added-animals, echo=FALSE, include=FALSE, eval=TRUE} # get-informative-and-uninformative-added-animals trimmedPedInformative <- trimPedigree(focalAnimals, breederPed, removeUninformative = TRUE ) uninformative <- trimmedPed$id[!trimmedPed$id %in% trimmedPedInformative$id] notFocalInTrimmedPed <- trimmedPed$id[!trimmedPed$id %in% focalAnimals] additionalAnimals <- nrow(trimmedPed) - length(focalAnimals) geneticallyInformative <- nrow(trimPedigree(focalAnimals, breederPed, removeUninformative = TRUE )) - length(focalAnimals) ``` As you can see from the number of rows in the full pedigree (`r nrow(breederPed)`) versus the trimmed pedigree (`r nrow(trimmedPed)`), trimmed pedigrees can be much smaller. Of the additional `r additionalAnimals` animals, `r geneticallyInformative` provide genetic information while the others (`r length(uninformative)`) are genetically uninformative. ```{r animals-no-birth-no-exit, echo = FALSE, eval = TRUE, include=FALSE} unknownBirth <- breederPed$id[is.na(breederPed$birth)] unknownBirthOrExit <- breederPed$id[is.na(breederPed$birth) | !is.na(breederPed$exit)] knownPed <- breederPed[!breederPed$id %in% unknownBirthOrExit, ] otherIds <- knownPed$id[!knownPed$id %in% trimmedPed$id[is.na(trimmedPed$exit)]] ``` As is shown below only `r length(otherIds)` (`r get_and_or_list(otherIds)`) living animals are still in the colony but are not in the trimmed pedigree.[^3] [^3]: Animals with an unknown birth date, and animals that have left the colony, are not considered in this comparison. ```{r show-living-animals-not-in-trimmed-pedigree, results = 'asis'} unknownBirth <- breederPed$id[is.na(breederPed$birth)] knownExit <- breederPed$id[!is.na(breederPed$exit)] unknownBirthKnownExit <- breederPed$id[is.na(breederPed$birth) | !is.na(breederPed$exit)] knownPed <- breederPed[!breederPed$id %in% unknownBirthKnownExit, ] otherIds <- knownPed$id[!knownPed$id %in% trimmedPed$id[is.na(trimmedPed$exit)]] print(stri_c( "The living animals in the pedigree that are not in the trimmed ", "pedigree are ", get_and_or_list(otherIds), "." )) ``` ## Age Sex Pyramid Plot You can examine the population structure using an age-sex pyramid plot with a single function. We will limit our view to just the focal animals and their living relatives. This is appropriate for colony management because in addition to the genetic diversity we seek, we have to remain cognizant of the age and sex distributions within the colonies we manage. ```{r plot-focal-age-sex-pyramid, include = TRUE} #| fig.alt: > #| Age Sex Pyramid Plot. getPyramidPlot(ped = trimmedPed[is.na(trimmedPed$exit), ]) ``` ## Pedigree Diagram The Pedigree Browser's Diagram tab in the Shiny application renders a pedigree as a family-tree diagram. A mate's own mating(s) render as a small connector between the two parents, with a line down to their shared children, rather than two independent lines running straight from each parent -- the same convention traditional pedigree charts use. An animal that mates more than once, or whose lineage loops back on itself (e.g. a consanguineous mating), appears once per mating, with each occurrence after the first joined back to its main occurrence by a dashed line. Both are built on two script-callable functions that work identically outside the Shiny application: __makePedigreeMatingLayout__ prepares a pedigree's node/edge data in this convention, and __drawPedigreeDiagram__ draws it with the __visNetwork__ package (a dependency of __nprcgenekeepr__); the Diagram tab itself calls __drawPedigreeDiagram__. (A simpler function, __makePedigreeDiagramData__, remains available for a plain one-node-per-animal diagram without this convention -- it is no longer what the Shiny app itself uses.) See the _pedigree-diagram_ article for the equivalent point-and-click workflow and screenshots of every feature -- including affected-status shading, name labels, and twin/zygosity connectors, none of which this scripted walkthrough demonstrates. That parent-to-mating-to-child connector has two available routings, selected by __makePedigreeMatingLayout__'s _edgeStyle_ argument: `"rectilinear"` (the default, Track 2, docs/planning/pedigree-diagram- kinship2-fidelity-remediation-plan.md -- strict horizontal/vertical right angles throughout: a horizontal segment between the two parents, a vertical drop to the mating dot, a horizontal sibship bar across the children, and a vertical drop into each one, matching the convention the __kinship2__ R package uses) and `"direct"` (a single straight or lightly sloped line all the way from parent to mating dot to child). The Shiny app's Diagram tab exposes the same choice live as a "Diagram Edge Style" toggle above the diagram; in a script it is just this one argument. Both are demonstrated below -- _direct_ first, _rectilinear_ (matching the tab's own default) afterward. The _trimmedPed_ pedigree built earlier in this tutorial has `r nrow(trimmedPed)` animals -- realistic for the live Shiny application's pan-and-zoom canvas, but too dense for a single static demonstration to usefully show every feature at once. It is also limited to the _Female_, _Male_, and _Unknown_ sex codes that happen to occur in the bundled example data -- the _Hermaphrodite_ and unrecorded-sex shapes never appear in it. So instead we build a small, synthesized pedigree spanning a few generations, deliberately including all five sex codes __makePedigreeMatingLayout__ recognizes -- and, since two of its founders (_M1_ and _F2_) each mate more than once, the diagram's duplicate-node convention as well.[^4] [^4]: This pedigree is entirely synthetic -- constructed for this tutorial, not drawn from any real colony's records. ```{r pedigree-diagram-demo-data} founders <- data.frame( id = c(paste0("F", 1L:4L), paste0("M", 1L:4L)), sire = NA_character_, dam = NA_character_, sex = c(rep("F", 4L), rep("M", 4L)), stringsAsFactors = FALSE ) gen1 <- data.frame( id = paste0("C", 1L:12L), sire = c( "M1", "M1", "M1", "M2", "M2", "M3", "M3", "M3", "M4", "M4", "M1", "M1" ), dam = c( "F1", "F1", "F1", "F2", "F2", "F3", "F3", "F3", "F4", "F4", "F2", "F2" ), sex = c("F", "M", "H", "F", "M", "M", "F", "U", "F", "M", "F", "M"), stringsAsFactors = FALSE ) gen2 <- data.frame( id = paste0("D", 1L:9L), sire = c("C2", "C2", "C5", "C5", "C5", "C6", "C10", "C10", "C12"), dam = c("C4", "C4", "C7", "C7", "C7", "C1", "C9", "C9", "C11"), sex = c("F", "M", "M", "F", "M", "F", "M", "F", "F"), stringsAsFactors = FALSE ) gen3 <- data.frame( id = paste0("E", 1L:4L), sire = c("D2", "D2", "D5", "D7"), dam = c("D4", "D4", "D8", "D6"), sex = c("F", "M", "U", NA_character_), stringsAsFactors = FALSE ) demoPed <- rbind(founders, gen1, gen2, gen3) demoPed$gen <- findGeneration(demoPed$id, demoPed$sire, demoPed$dam) nrow(demoPed) table(demoPed$sex, useNA = "ifany") ``` Sex code __H__ (_C3_) and the unrecorded sex code (_E4_, `NA`) each appear exactly once, alongside code __U__ (_C8_ and _E3_) and the ordinary __F__/ __M__ codes that make up the rest -- enough for every one of the five node shapes to appear in the diagram below at least once, across `r length(unique(demoPed$gen))` generations and `r nrow(demoPed)` animals. ### Direct Edge Style ```{r pedigree-diagram-data} ## 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) nrow(diagramData$nodes) nrow(diagramData$edges) ``` __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.[^5] [^5]: The Diagram tab's click-to-navigate behavior (clicking a node narrows the focal-animal selection to it, resolving a duplicate-occurrence click to the real animal and ignoring mating-node clicks) is Shiny-specific -- it relies on a live Shiny session to receive the click event -- and is therefore not reproduced here. ```{r pedigree-diagram-render} 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 `r nrow(trimmedPed)` 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). ```{r pedigree-diagram-rectilinear-data} diagramDataRectilinear <- makePedigreeMatingLayout(demoPed, edgeStyle = "rectilinear") names(diagramDataRectilinear) nrow(diagramDataRectilinear$nodes) nrow(diagramDataRectilinear$edges) ``` `r nrow(diagramDataRectilinear$nodes)` nodes and `r nrow(diagramDataRectilinear$edges)` edges here, versus `r nrow(diagramData$nodes)` and `r nrow(diagramData$edges)` for the same `r nrow(demoPed)`-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.[^6] _nodes_ also gains an extra `color.border` column under _rectilinear_ -- __visNetwork::visNetwork()__ reads it directly, so no extra rendering code is needed. [^6]: These waypoint nodes are internal plumbing needed for the right-angle routing itself, not new diagram content. 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. ```{r pedigree-diagram-rectilinear-render} 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. ```{r twin-relations-setup} 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: ```{r twin-relations-demo} 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: ```{r twin-kinship} 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")] twinKmat[c("T1", "T2", "D1"), c("T1", "T2", "D1")] ``` The twins' kinship rises from `r plainKmat["T1", "T2"]`, that of full sisters, to `r twinKmat["T1", "T2"]`, and the daughter's kinship with her mother's twin rises from `r plainKmat["D1", "T2"]`, that of a niece and her aunt, to `r twinKmat["D1", "T2"]`, 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. ```{r layout-kinship-matrix} demoKmat <- kinship(demoPed$id, demoPed$sire, demoPed$dam, demoPed$gen) layoutFromPed <- makePedigreeMatingLayout(demoPed) layoutFromMatrix <- makePedigreeMatingLayout(demoPed, kinshipMatrix = demoKmat) identical(layoutFromMatrix$edges, layoutFromPed$edges) unnamedMatrixMessage <- tryCatch( makePedigreeMatingLayout(demoPed, kinshipMatrix = unname(demoKmat)), error = conditionMessage ) unnamedMatrixMessage ``` 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_. ## Genetic Value Analysis Your genetic value analysis must be carefully performed. The next three commands set up the entire pedigree for analysis. The first of these three commands set all of the pedigree members to be part of the population of interest by setting the _population_ column to __TRUE__ for all individuals. ```{r set-entire-pedigree-as-population} ped <- setPopulation(breederPed, NULL) ``` Note that a new pedigree object (__ped__) is being created. ```{r full-pedigree-genetic-value-summary} probands <- ped$id[ped$population] ped <- trimPedigree(probands, ped, removeUninformative = FALSE, addBackParents = FALSE ) ``` The arguments to __reportGV__ are all optional except for _ped_, but you may often want to use non-default values. The principal arguments are described below. Further arguments (__updateProgress__, `breedingTable`, `gestationTable`, `breedingAgeDefault`, `gestationDefault`, `kinshipOverrides`, __twinRelations__, `guCutoff`, `zScoreCutoff` and `axisPriority`) are described in `?reportGV`. - __ped__ Pedigree information in data.frame format - __guIter__ Integer indicating the number of iterations for the gene-drop analysis. Default is 1000 iterations - __guThresh__ Integer indicating the threshold number of animals for defining a unique allele. Default considers an allele "unique" if it is found in only 1 animal. - __pop__ Character vector with animal IDs to consider as the population of interest. The default is NULL. - __byID__ Logical variable of length 1 that is passed through to eventually be used by alleleFreq(), which calculates the count of each allele in the provided vector. If byID is TRUE and ids are provided, the function will only count the unique alleles for an individual (homozygous alleles will be counted as 1). ```{r full-geneticValue} geneticValue <- reportGV(ped, guIter = 50L, guThresh = 3L, byID = TRUE, updateProgress = NULL ) summary(geneticValue) ``` What happens if we limit our analysis to the trimmed pedigree? Remember that the trimmed pedigree still contains all of the genetic information that the full pedigree has for the focal animals. ```{r trimmed-genetic-value-analysis} trimmedGeneticValue <- reportGV(trimmedPed, guIter = 50L, guThresh = 3L, byID = TRUE, updateProgress = NULL ) summary(trimmedGeneticValue) ``` It is clear that limiting your analysis to the animals of interest reduces the effort required to examine the animals of interest. ### 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. ```{r genotypes-from-label-columns} 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))) 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) ``` ### 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. ```{r list-genetic-value-objects} names(trimmedGeneticValue) ``` The _report_ object (an R dataframe) can in-turn be examined. ```{r list-report-object-parts} names(trimmedGeneticValue$report) ## column names nrow(trimmedGeneticValue$report) ## Number of rows ``` 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. ```{r look-at-genetic-value-report} 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. ```{r kinship-and-founders} 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. ```{r genetic-uniqueness-boxplot} #| fig.alt: > #| Genetic Uniqueness Box Plot. 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. ```{r extraction-of-mk-zs} 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. ```{r kinship-x-chromosome} 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. ## Colony Snapshots Over Time A _snapshot_ is one dated row of colony-level numbers: the settings and counts that say what the numbers rest on, the values __reportGV__ reports for the colony as a whole, and summaries of the per-animal kinship and genome uniqueness values. Rows for successive dates make a _history_, which is what the Genetic-Health Trends tab of the Shiny application keeps and plots. The functions in this section do the same work in a script. The chunks below use the example pedigree that ships with the package, __qcPed__, and few gene-drop iterations so that each runs in a moment. Use more for an analysis you will report (the default of __reportGV__ is `r formals(reportGV)$guIter`). ### 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. ```{r snapshot-record} 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 `r ncol(firstSnapshot)` 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. ```{r snapshot-backfill} 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 `r nrow(earlierSnapshots$skipped)` 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. ```{r snapshot-save-read} historyFile <- tempfile(fileext = ".csv") write.csv(snapshotHistory, historyFile, row.names = FALSE) readBack <- readSnapshotHistory(historyFile) class(readBack$snapshotDate) checkedHistory <- checkSnapshotHistory(readBack) class(checkedHistory$snapshotDate) ``` ### 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. ```{r snapshot-compare} mixedRuleMessage <- tryCatch( calcSnapshotDeltas(checkedHistory, from = "2005-01-01", to = "2010-01-01"), error = conditionMessage ) mixedRuleMessage 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. ```{r snapshot-plot} #| fig.height: 4 #| fig.alt: > #| Trend lines of three colony genetic-health measures over the snapshot #| dates. The rebuilt snapshots form a dashed line and the one recorded #| snapshot is a single point. trendPlot <- plotSnapshotTrends(checkedHistory, metrics = c("fe", "fg", "meanGu"), ncol = 3L) trendPlot ``` ## Breeding Group Formation The primary purpose of __nprcgenekeepr__ is to form breeding groups according to our best information regarding maintaining the genetic characteristics we desire and the realities associated with other animal husbandry needs. There are several options you must consider when forming groups using __nprcgenekeepr__, which we will examine using code below. - Animals used to form groups - _high-value_: Randomly select from only high-value animals in genetic value analysis - _all_: Randomly select from all animals in genetic value analysis - _candidates_: Use candidate animals entered below to form groups - Sex ratio - Randomly assign animals without regard to sex - Use a harem structure with one breeding male per group - Specify the sex ratio between 0.5 and 20 (F/M) - Whether or not to pre-populate (seed) groups with animals of your choice - Number of groups to be formed - Whether or not to ignore females at or above the minimum parent age - Number of simulations used to search for the optimal group makeup - Whether or not kinship coefficients are to be included in results You decisions regarding each of the above options are expressed in a call to the function __groupAddAssign__. A complete description of the function and its arguments is available using the code shown below. ```{r groupAddAssign-help-request, eval = FALSE} ?groupAddAssign ``` Below are the descriptions of the principal function parameters extracted from the documentation near the time this tutorial was prepared. The function has further parameters (`maxCandidates`, __exhaustive__, `maxExhaustiveCandidates`, `exhaustiveTimeLimit`, and __ancestryRules__, which is used in the _Ancestry Rules_ section below), and it returns a list with the items __group__, __score__, __candidates__ (up to `maxCandidates` distinct candidate solutions, best first) and, when __withKin__ is __TRUE__, `groupKin`; see `?groupAddAssign` for their descriptions. - __candidates__ Character vector of IDs of the animals available for use in forming the groups. The animals that may be present in currentGroups are not included within candidates. - __currentGroups__ List of character vectors of IDs of animals currently assigned to groups. Defaults to a list with character(0) in each sub-list element (one for each group being formed) assuming no groups are pre-populated. - __kmat__ Numeric matrix of pairwise kinship values. Rows and columns are named with animal IDs. - __ped__ Dataframe that is the 'Pedigree'. It contains pedigree information including the IDs listed in candidates. - __threshold__ Numeric value indicating the minimum kinship level to be considered in group formation. Pairwise kinship below this level will be ignored. The default values is 0.015625. - __ignore__ List of character vectors representing the sex combinations to be ignored. If provided, the vectors in the list specify if pairwise kinship should be ignored between certain sexes. Default is to ignore all pairwise kinship between females. - __minAge__ Integer value indicating the minimum age to consider in group formation. Pairwise kinships involving an animal of this age or younger will be ignored. Default is 1 year. - __iter__ Integer indicating the number of times to perform the random group formation process. Default value is 1000 iterations. - __numGp__ Integer value indicating the number of groups that should be formed from the list of IDs. Default is 1. - __updateProgress__ Function or NULL. If this function is defined, it will be called during each iteration to update a shiny::Progress object. - __harem__ Logical variable when set to TRUE, the formed groups have a single male at least minAge old. - __sexRatio__ Numeric value indicating the ratio of females to males x from 0.5 to 20 by increments of 0.5. - __withKin__ Logical variable when set to TRUE, the kinship matrix for the group is returned along with the group and score. Defaults to not return the kinship matrix. This maintains compatibility with earlier versions. We will use the _trimmedPed_ pedigree in our code. For illustration purposes we are going to keep the numbers of candidates, groups, and iterations fairly small. We will get first some animal IDs to use for our candidates by selecting animals at least 2 years old at the time this pedigree was sampled (01-01-2015) that are still in the colony. ```{r get-candidates} candidates <- trimmedPed$id[trimmedPed$birth < as.Date("2013-01-01") & !is.na(trimmedPed$birth) & is.na(trimmedPed$exit)] table(trimmedPed$sex[trimmedPed$id %in% candidates]) ``` Our candidates are made up of `r table(trimmedPed$sex[trimmedPed$id %in% candidates])[[1]]` females and `r table(trimmedPed$sex[trimmedPed$id %in% candidates])[[2]]` males. The parameters __currentGroups__, __threshold__, __ignore__, __minAge__, __sexRatio__, __withKin__, and __updateProgress__ are allowed to take their default values in the calls below. The default __sexRatio__ of 0 means that no sex ratio is imposed. ### Harems The following group assignments will be forming harem groups. This is done by setting __harem__ to \code{TRUE}. Setting `iter` to 100 or more will increase optimal composition of breeding groups ```{r create-harem-groups} haremGrp <- groupAddAssign( kmat = trimmedGeneticValue[["kinship"]], ped = trimmedPed, candidates = candidates, iter = 10L, numGp = 6L, harem = TRUE ) haremGrp$group ``` 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. ```{r list-males-in-harem-groups} sapply(haremGrp$group, function(ids) { ids[ids %in% trimmedPed$id[trimmedPed$sex == "M"]] }) ``` 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. ```{r harem-count-and-sexratios} lines <- sapply(haremGrp$group, function(ids) { paste0( "Count: ", length(ids), " Sex Ratio: ", round(calculateSexRatio(ids, trimmedPed), 2) ) }) for (line in lines) print(line) ``` Examination of this table shows that of the `r table(trimmedPed$sex[trimmedPed$id %in% candidates])[[1]]` females `r sum(unlist(haremGrp$group[1:6]) %in% trimmedPed$id[trimmedPed$sex == "F"])` 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__). ```{r create-sexratio-groups} sexRatioGrp <- groupAddAssign( kmat = trimmedGeneticValue[["kinship"]], ped = trimmedPed, candidates = candidates, iter = 10L, numGp = 6L, sexRatio = 9.0 ) sexRatioGrp$group ``` As before, the last element of the list is not a breeding group. Here it holds the `r length(sexRatioGrp$group[[7]])` 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. ```{r list-males-in-sexratio-groups} sapply(sexRatioGrp$group, function(ids) { ids[ids %in% trimmedPed$id[trimmedPed$sex == "M"]] }) ``` We can also see the number of animals and the sex ratios created in each group. ```{r list-count-and-sexratios} lines <- sapply(sexRatioGrp$group, function(ids) { paste0( "Count: ", length(ids), " Sex Ratio: ", round(calculateSexRatio(ids, trimmedPed), 2L) ) }) for (line in lines) print(line) ``` Examination of this table shows that of the `r table(trimmedPed$sex[trimmedPed$id %in% candidates])[[1]]` females `r sum(unlist(sexRatioGrp$group[1:6]) %in% trimmedPed$id[trimmedPed$sex == "F"])` are included in the six groups. ## Individual Mate-Pair Analysis __groupAddAssign__, used throughout the previous section, forms whole breeding groups. __reportMatePairs__ instead answers a narrower question -- for a candidate population, which individual sire/dam pairs are even eligible to be considered (opposite sex, both at or above __minAge__, or of unknown age), and what does every available kinship/genetic-value signal say about each pair? It composes the package's existing pair-eligibility machinery into a single report, returning two tables: _pairs_, one row per eligible sire/dam combination, and _excluded_, one row per pair dropped and why. It computes no composite ranking score of its own -- sort or filter the returned columns as your own workflow requires. Marker-based kinship (__markerKinship__'s output) and each parent's genetic-value context (_indivMeanKin_/_gu_, from __reportGV__'s report) can both be supplied as optional enrichment and default to `NA` when not. We reuse the _candidates_ vector already built above, restricted further to a small illustrative slice -- one male and seven females -- so the output stays readable. One of those females is not in the kinship matrix, so six pairs come back below and the _excluded_ table is empty. A full call against all `r length(candidates)` candidates would return one row per opposite-sex, age-eligible pair, which runs into the tens of thousands. ```{r mate-pairs-setup} smallCandidates <- head(candidates, 8L) trimmedPed[trimmedPed$id %in% smallCandidates, c("id", "sex", "age")] ``` ```{r mate-pairs-demo} matePairs <- reportMatePairs(trimmedPed, trimmedGeneticValue[["kinship"]], geneticValues = trimmedGeneticValue, populationIds = smallCandidates) matePairs$pairs ``` With only one male in this slice, every eligible pair shares the same sire, so _sireIndivMeanKin_/_sireGu_ are identical across rows (both come from _trimmedGeneticValue_'s report for that one animal); _damIndivMeanKin_/ _damGu_ vary by dam. _markerKinship_ is `NA` throughout, since no marker genotype matrix was supplied to this call. __exclude__ removes a specific id from consideration and records why in the _excluded_ table -- useful when a curator already knows a particular animal is unavailable this cycle, independent of age or sex eligibility. ```{r mate-pairs-exclude-demo} matePairsExcluded <- reportMatePairs(trimmedPed, trimmedGeneticValue[["kinship"]], geneticValues = trimmedGeneticValue, populationIds = smallCandidates, exclude = matePairs$pairs$damId[1L]) matePairsExcluded$excluded ``` One caveat worth knowing before relying on __minAge__ alone to bound a real-world call: a missing (`NA`) recorded age *passes* the age screen rather than failing it, on the premise that an unknown age should not silently disqualify an otherwise-eligible animal. __populationIds__, not __minAge__, is the reliable way to bound which candidates a call considers. ### 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. ```{r mate-pairs-ancestry-setup} 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 ``` ```{r mate-pairs-ancestry-demo} ancMatePairs <- reportMatePairs(ancPed, ancKmat, ancestryRules = ancestryRules) names(ancMatePairs) ``` 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_. ```{r mate-pairs-ancestry-tables} 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. ```{r mate-pairs-ancestry-coverage} 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. ```{r mate-pairs-ancestry-override} overriddenMatePairs <- reportMatePairs( ancPed, ancKmat, ancestryRules = ancestryRules, overriddenRules = data.frame(ancestry1 = "CHINESE", ancestry2 = "INDIAN")) table(overriddenMatePairs$pairs$ancestryStatus, useNA = "ifany") table(overriddenMatePairs$excluded$reason) ``` 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. ```{r ancestry-violations} 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. ```{r ancestry-violations-overridden} ancViolationsOverridden <- reportAncestryViolations( formedGroups, ancPed, ancestryRules, overriddenRules = data.frame(ancestry1 = "CHINESE", ancestry2 = "INDIAN")) ancViolationsOverridden$violations ``` ## Pedigree Errors As stated earlier you can see which types of errors are detected by _qcStudbook_ by looking at names returned by _getEmptyErrorLst()_ as shown below. ```{r getEmptyErrorLst} names(getEmptyErrorLst()) ``` Each is defined below. ```{r make-errorList-definition-tbl, echo = FALSE, eval=TRUE} errorTypes <- names(getEmptyErrorLst()) errorDescriptions <- c( "Database connection failed: configuration or permissions are invalid", "Columns that must be within the pedigree file are missing.", "Values, which are supposed to be dates, cannot be interpreted as a date.", "Parents were too young on the date of birth of the offspring to have been the parent.", "Individuals listed as female and as a sire.", "Individuals are listed as male and as a dam.", "Individuals who are listed as both a sire and a dam.", "IDs listed more than once.", "IDs (id, sire, or dam) that contain a period ('.'), which is not allowed.", stri_c( "Rows with a placeholder value other than TRUE, FALSE, 1, 0 or blank ", "(the optional column that marks stand-ins for unknown parents)." ), stri_c( "Columns that have been changed to conform to internal naming ", "conventions and what they were changed to." ) ) errorTbl <- data.frame( Error = errorTypes, Definition = errorDescriptions, stringsAsFactors = FALSE ) ``` ```{r print-error-definition-tbl, echo=FALSE} knitr::kable(errorTbl) ``` We are going to use the small imaginary pedigree listed below that has multiple errors to discuss pedigree error detection and reporting. First note the birth dates of ego_id _o4_ (2006-04-13) and the purported sire _s2_ (2006-06-19). Note also the purported birth date of the _d2_ and the birth dates of her offspring. Obviously dates or IDs may be incorrect. This is the pedigree. _(We will discuss the column names shortly.)_ ```{r list-pedOne} knitr::kable(nprcgenekeepr::pedOne) ``` If we try to convert this pedigree file into the standardized studbook format, we are going to get an error message and the creation of a file in the R sessions temporary directory that lists the records that have generated the errors. ```{r summary-pedOne-no-errors, error = TRUE} pedOne <- nprcgenekeepr::pedOne # put it in the local environment ped <- qcStudbook(pedOne, minSireAge = 0.0, minDamAge = 0.0) ``` The contents of _lowParentAge.csv_ is shown below. ```{r read-lowParentAge.csv, echo = FALSE, include = FALSE} lowParentAge <- read.csv(paste0(tempdir(), "/lowParentAge.csv")) ``` ```{r print-lowParentAge.csv, echo = FALSE} knitr::kable(lowParentAge) ``` Examination of the ages of the parents reveals the issues being reported. We can remove the errors by changing the birth dates of _o4_ from 2006-04-13 to 2015-09-16 and of _d2_ from 2015-09-16 to 2006-04-13. ```{r correct-birthdate, error = TRUE} pedOne$birth_date[pedOne$ego_id == "o4"] <- as.Date("2015-09-16") pedOne$birth_date[pedOne$ego_id == "d2"] <- as.Date("2006-04-13") ``` Note the changes made to the column names between the original __pedOne__ pedigree and the pedigree (__ped__) we get from __qcStudbook__. We have chosen to limit the displayed pedigree by selecting the _ego_id_'s and _id_'s in __pedOne__ and __ped__ respectively. ```{r display-corrected-birth-records} ped <- qcStudbook(pedOne, minSireAge = 0.0, minDamAge = 0.0) ped[ped$id %in% c("s2", "d2", "o3", "o4"), ] ``` However, the preferred method of creating the standardized studbook format with __qcStudbook__ is to examine all errors found and correcting them before proceeding. This is done by explicitly requesting that all aspects inconsistent with the studbook format be identified by setting _reportChanges_ and _reportErrors_ to \code{TRUE}. ```{r look-for-errors} errorList <- qcStudbook(pedOne, minSireAge = 0.0, minDamAge = 0.0, reportChanges = TRUE, reportErrors = TRUE ) summary(errorList) ``` We will discuss each of these newly identified errors in a moment, however, let's look at shortening this report, because often you will not be interested in the more trivial changes to the column names made by __qcStudbook__ and in those cases you simply choose not to report changes to the column names as is shown here by setting _reportChanges_ to \code{FALSE}. For this illustration, we are going to bring back the original copy of __pedOne__ to see how the suspicious parents are reported by the __summary__ function. ```{r look-for-errors-only} pedOne <- nprcgenekeepr::pedOne errorList <- qcStudbook(pedOne, minSireAge = 0L, minDamAge = 0L, reportChanges = FALSE, reportErrors = TRUE ) options(width = 90L) summary(errorList) ``` The first two errors mentioned are of particular interest. Currently __qcStudbook__ automatically changes the sex of dams to _F_ (female) and sires to _M_ (male) when __reportErrors__ is set to \code{FALSE}. ## Genetic Loops This feature is not supported within the Shiny application and is not fully implemented. To use the __findLoops__ function run the following code on a pedigree that has a loop in it. We are continuing to use the example pedigree that comes with the software, _ExamplePedigree.csv_, which is held in memory as `breederPed`. ```{r look-at-loops} exampleTree <- createPedTree(breederPed) exampleLoops <- findLoops(exampleTree) ``` You can count the animals checked, the animals that are in loops, and the total of the per-animal loop counts returned by __countLoops__ with the following code. ```{r countLoops} cat("Animals checked:", length(exampleLoops), "\n") cat("Animals in loops:", sum(unlist(exampleLoops)), "\n") nLoops <- countLoops(exampleLoops, exampleTree) cat("Total of per-animal loop counts:", sum(unlist(nLoops[nLoops > 0L])), "\n") ``` You can list the first 10 sets of ids, sires and dams in loops with the following line of code: ```{r listLoops} loopIds <- names(exampleLoops)[unlist(exampleLoops)] breederPed[match(loopIds, breederPed$id), c("id", "sire", "dam")][1L:10L, ] ``` ## Shrinking a Pedigree to a Genotyping Budget __trimPedigree__ cuts a pedigree down around the focal animals you name. __shrinkPedigree__ answers a different question: given which animals are genotyped, which animals can be dropped and still leave a pedigree that is informative within a genotyping-cost budget? It works on this package's pedigree data frame (the _id_, _sire_ and _dam_ columns). The budget, _maxBits_, limits the pedigree's bit size, which is twice the number of non-founders less the number of founders. Three trims are applied in order: - The _unavailable_ trim repeatedly removes animals that are not genotyped and have no offspring in the pedigree, then founder pairs that are not genotyped and have a single child and no other mate (that child becomes a founder), then any founder left without offspring, genotyped or not. - The _non-informative_ trim removes a genotyped animal that is nobody's parent, has both parents known and genotyped, and is not affected, because it adds no genotype information beyond what its parents supply. - The _affected-priority_ trim, while the bit size is still above _maxBits_, removes one genotyped animal that is nobody's parent at a time, choosing the removal that leaves the smallest bit size. It tries animals of unknown affected status first, then unaffected ones, then affected ones. The optional _affected_ vector marks the animals that show the trait you care about; without it every animal counts as unaffected. _genotyped_ is a logical vector with one value per row of the pedigree, `NA` not allowed. Suppose only the animals still living, those with no exit date, are genotyped: ```{r shrink-pedigree} genotyped <- is.na(qcPed$exit) shrinkResult <- shrinkPedigree(qcPed[, c("id", "sire", "dam")], genotyped) unlist(shrinkResult[c("pedSizeOriginal", "pedSizeIntermed", "pedSizeFinal")]) lengths(shrinkResult$idList) ``` The `r shrinkResult$pedSizeOriginal` animals of __qcPed__ shrink to `r shrinkResult$pedSizeIntermed` after the first two trims and `r shrinkResult$pedSizeFinal` at the end. _idList_ groups the removed ids by the trim that took them: _unavail_, _noninform_ and _affected_. _idTrimmed_ lists every removed id in the order of removal, and _bitSize_ records the bit size before any trimming, after the first two trims, and after each round of the third. The default _maxBits_ is `r formals(shrinkPedigree)$maxBits`; a smaller budget asks the third trim to go further. ```{r shrink-pedigree-budget} shrinkResultSmall <- shrinkPedigree(qcPed[, c("id", "sire", "dam")], genotyped, maxBits = 8L) c(default = nrow(shrinkResult$ped), smaller = nrow(shrinkResultSmall$ped)) ``` The smaller budget leaves `r nrow(shrinkResultSmall$ped)` animals in place of `r nrow(shrinkResult$ped)`. ## Marker Genetics Everything up to this point has estimated kinship and diversity from a _recorded_ pedigree -- __kinship__, __reportGV__, and the rest all assume the sire/dam relationships in your pedigree data are correct. The functions in this section instead estimate relatedness, heterozygosity, parentage consistency, and between-center differentiation directly from a marker genotype panel (e.g. a SNP or STR panel), independent of any known pedigree. This is the same functionality behind the Shiny application's __Marker Genetics__ tab -- see the "Marker Genetics" section of the _colony-manager-guide_ article for the equivalent point-and-click workflow and screenshots. Several of the small example genotypes below (kinship, heterozygosity, exclusion, and Fst) match the ones used there, so those numbers are the same ones shown in that article's tables; the candidate-parent likelihood example here differs from the article's. ### 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. ```{r marker-genotype-setup} 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. ```{r marker-genotype-matrix} markerGenotype <- checkMarkerGenotypeFile(markerGenotype) genotypeMatrix <- buildMarkerGenotypeMatrix(markerGenotype) genotypeMatrix ``` ### 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. ```{r marker-kinship-demo} markerKmat <- markerKinship(genotypeMatrix) markerKmat ``` _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_. ```{r marker-heterozygosity-setup} 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) ``` ```{r marker-heterozygosity-demo} markerObservedHeterozygosity(hetMatrix) markerExpectedHeterozygosity(hetMatrix) ``` 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_. ```{r marker-parentage-setup} pedigree <- data.frame(id = c("P", "C", "U"), sire = c(NA, "U", NA), dam = c(NA, "P", NA), stringsAsFactors = FALSE) ``` ```{r marker-parentage-demo} 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. ```{r parentage-likelihood-setup} 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: ```{r parentage-exclusion-recheck} 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: ```{r parentage-likelihood-demo} 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. ```{r cross-center-setup} 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) ``` ```{r cross-center-mapping-demo} 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: ```{r cross-center-mapping-bad-demo} 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. ```{r cross-center-identity-demo} 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. ```{r marker-fst-setup} 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)) ``` ```{r marker-fst-demo} markerFst(matrixA, matrixB) ``` _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. ```{r linkage-genotype-setup} 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: ```{r linkage-genotype-demo} checkLinkageMarkerGenotypeFile(strGenotype) ``` The same data fed to __checkMarkerGenotypeFile__ instead fails, exactly as that function is designed to: ```{r biallelic-genotype-fails-demo, error = TRUE} checkMarkerGenotypeFile(strGenotype) ``` 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. ```{r locus-metadata-demo} 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). ```{r realized-relatedness-variance-demo} 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. ```{r ld-block-setup} 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 )) ``` ```{r ld-block-demo} 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. ```{r ld-block-obfuscate-setup} 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) ``` ```{r ld-block-obfuscate-demo} 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. ```{r obfuscate-id-placeholder-demo} ids <- c("U1234", "A100", "U0007", "B200") isPlaceholder <- c(FALSE, FALSE, TRUE, FALSE) set.seed(1) obfuscateId(ids, size = 6L) set.seed(1) obfuscateId(ids, size = 6L, placeholder = isPlaceholder) ``` 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. ## MHC Haplotype Reporting The MHC Haplotype Reporting tab of the Shiny application has four counterparts for scripts. The input is a wide table with one row for each animal: the animal's id, then its two haplotype designations. __checkMhcHaplotypeFile__ checks that structure (exactly three columns, the first an id, and one row for each animal) and gives the columns fixed names whatever the file called them; the example data set __rhesusGenotypes__ calls the designation columns _first_name_ and _second_name_. A designation is an opaque label: it is never split or compared with MHC region names. A trailing `?` marks an uncertain call, and an empty cell is a missing one. ```{r mhc-check} mhcChecked <- checkMhcHaplotypeFile(rhesusGenotypes) names(rhesusGenotypes) names(mhcChecked) head(mhcChecked, 3L) ``` __mhcHaplotypeFrequency__ summarizes the table. The frequency of a haplotype is its number of copies among the certain calls: twice the animals, less the missing and the uncertain calls. Uncertain calls are left out of the copy counts, the carrier counts and the frequencies, but they are never hidden: _nUncertain_ gives them for each haplotype and the _counts_ table for the file as a whole. A haplotype is flagged rare (_isRare_) when its frequency is at or below _rareFrequencyThreshold_, default `r formals(mhcHaplotypeFrequency)$rareFrequencyThreshold`, or its carrier count is at or below _rareCarrierThreshold_, default `r formals(mhcHaplotypeFrequency)$rareCarrierThreshold`. Both defaults are working filters, not scientific claims. ```{r mhc-frequency} mhcFrequency <- mhcHaplotypeFrequency(mhcChecked) mhcFrequency$counts head(mhcFrequency$summary) ``` In this colony of `r mhcFrequency$counts$nAnimals` animals, `r sum(mhcFrequency$summary$isRare)` of the `r nrow(mhcFrequency$summary)` haplotypes are flagged rare. The smallest frequency, `r round(min(mhcFrequency$summary$frequency), 3L)`, is above the frequency threshold, so here only the carrier count can flag a haplotype. __mhcHaplotypeCarriers__ goes from which haplotypes are rare to which animals carry them: one row for each haplotype and carrying animal. By default it lists only the haplotypes flagged rare; `rareOnly = FALSE` lists them all. The _uncertain_ column is `TRUE` when the animal's only calls of that haplotype are uncertain, a provisional carrier that is listed here but is not counted as a carrier above. ```{r mhc-carriers} rareCarriers <- mhcHaplotypeCarriers(mhcChecked) allCarriers <- mhcHaplotypeCarriers(mhcChecked, rareOnly = FALSE) head(rareCarriers) c(rare = nrow(rareCarriers), all = nrow(allCarriers)) ``` A carrier table names real animals, and a paired MHC type can identify an animal on its own. __obfuscateMhcHaplotypes__ replaces the ids in a carrier table with the aliases that __obfuscatePed__ gives the same animals (its _map_ output) and leaves the haplotype labels and the _uncertain_ column unchanged, because a haplotype name is shared nomenclature and not an animal identifier. An id missing from the map stops the call, so a real id cannot slip through. ```{r mhc-obfuscate} set.seed(1L) obfuscated <- obfuscatePed(rhesusPedigree[, c("id", "sire", "dam", "sex")], map = TRUE) hiddenCarriers <- obfuscateMhcHaplotypes(rareCarriers, obfuscated$map) head(hiddenCarriers) ``` The aliases are drawn at random, so the seed fixes them for this document. ```{r elapsed-time} elapsed_time <- get_elapsed_time_str(start_time) ``` The current date and time is `r Sys.time()`. The processing time for this document was `r elapsed_time`. ```{r session-info} sessionInfo() ```