Initial commit

.Rbuildignore

@@ -0,0 +1,7 @@
+^segregatr\.Rproj$
+^\.Rproj\.user$
+^notes-private$
+^cran-comments\.md$
+^CRAN-RELEASE$
+^revdep$
+^data-raw$

.gitignore

@@ -0,0 +1,7 @@
+.Rproj.user
+.Rhistory
+.RData
+notes-private
+cran-comments.md
+revdep/
+CRAN-RELEASE

DESCRIPTION

@@ -0,0 +1,22 @@
+Package: segregatr
+Title: Segregation Analysis for Variant Pathogenicity Classification
+Version: 0.1.0
+Authors@R: 
+    person(given = "Magnus Dehli",
+           family = "Vigeland",
+           role = c("aut", "cre"),
+           email = "[email protected]",
+           comment = c(ORCID = "0000-0002-9134-4962"))
+Description: Computation of Bayes factors for evaluating segregation evidence
+    for pathogenicity of a genetic variant. 
+License: GPL-3
+Encoding: UTF-8
+Language: en-GB
+LazyData: true
+Suggests: 
+    testthat (>= 2.1.0),
+Imports: 
+    pedtools,
+    pedprobr
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.1.0

NAMESPACE

@@ -0,0 +1,5 @@
+# Generated by roxygen2: do not edit by hand
+
+export(FLB)
+import(pedtools)
+importFrom(pedprobr,likelihood)

R/FLB.R

@@ -0,0 +1,99 @@
+#' Full-likelihood Bayes factor
+#'
+#' Computes the Bayes factor for co-segregation, as described by Thompson et al
+#' (2003).
+#'
+#' @param x A `ped` object
+#' @param carriers A character vector (or coercible to such), containing the ID
+#'   labels of pedigree members known to carry the variant in question.
+#' @param noncarriers A character vector (or coercible to such), containing the
+#'   ID labels of pedigree members known *not* to carry the variant in question.
+#' @param freq A single number; the allele frequency of the variant in question.
+#' @param affected The affected pedigree members.
+#' @param unknown Pedigree members with unknown affection status.
+#' @param proband The ID label of the proband. This person must also be in both
+#'   `carriers` and `affected`.
+#' @param penetrances Either a numeric vector of length 3 (f0, f1, f2) or a data
+#'   frame with 3 or more columns. Each row contains the penetrance values of a
+#'   liability class. The first three columns are interpreted as penetrance
+#'   values f0, f1, f2 respectively; additional columns are ignored.
+#' @param liabilityClasses A vector of length `pedsize(x)`, containing for each
+#'   pedigree member the row number of `penetrances` which should be used for
+#'   that individual. (If `penetrances` is just a vector, it will be used for
+#'   all classes.) If `liabilityClasses` is NULL (the default), it is set to `1`
+#'   for all individuals.
+#' @param details A logical, indicating if detailed output should be returned.
+#' @param plot A logical.
+#'
+#' @return A positive number. If `details = TRUE`, a list of intermediate
+#'   results is returned.
+#'
+#' @examples
+#' FLB(pedtools::nuclearPed(2), carriers = 3:4, aff = 3:4, unknown = 1:2,
+#'     freq = 0.0001, penetrances = c(0, 1, 1), proband = 3)
+#'
+#' @export
+FLB = function(x, carriers, noncarriers = NULL, freq,
+               affected, unknown = NULL, proband,
+               penetrances, liabilityClasses = NULL,
+               details = FALSE, plot = FALSE) {
+
+  labs = labels(x)
+
+  if(!proband %in% affected)
+    stop2("The proband must be affected")
+  if(!proband %in% carriers)
+    stop2("The proband must be a carrier")
+
+  # Affection status vector, sorted along labs
+  aff = logical(pedsize(x))
+  aff[internalID(x, affected)] = TRUE
+  aff[internalID(x, unknown)] = NA
+
+  # Empty marker (= disease locus)
+  dis = m = mProband = marker(x, afreq = c(a = 1 - freq, b = freq))
+
+  # Full marker
+  genotype(m, carriers) = c("a", "b")
+  genotype(m, noncarriers) = c("a", "a")
+
+  # Proband marker
+  genotype(mProband, proband) = c("a", "b")
+
+  if(plot)
+    plot(x, m, skip.empty = T, shaded = labs[aff], starred = proband)
+
+  # Utility for setting up likelihood under causative hypothesis
+  quickStart = function(locus)
+    startdata_causative(x, marker = locus, aff = aff, penetrances = penetrances,
+                        liabilityClasses = liabilityClasses)
+
+  # Main Bayes factor
+  peelOrder = peelingOrder(x)
+  setup1 = list(informativeNucs = peelOrder, startdata = quickStart(m))
+  numer1 = likelihood(x, m, setup = setup1)
+
+  setupDis = list(informativeNucs = peelOrder, startdata = quickStart(dis))
+  likDis = likelihood(x, dis, setup = setupDis)
+  likM = likelihood(x, m)
+
+  denom1 = likDis * likM
+  BF1 = numer1/denom1
+
+  # Correction factor
+  likMproband = likelihood(x, mProband)
+  numer2 = likDis * likMproband
+
+  setup2 = list(informativeNucs = peelOrder, startdata = quickStart(mProband))
+  denom2 = likelihood(x, mProband, setup = setup2)
+
+  BF2 = numer2/denom2
+
+  FLB = BF1 * BF2
+
+  if(details)
+    return(list(c(FLB = FLB, BF1 = BF1, BF2 = BF2),
+                c(numer1 = numer1, denom1 = denom1, numer2 = numer2, denom2 = denom2),
+                c(likM = likM, likMproband = likMproband, likDis = likDis)))
+  FLB
+}

R/segregatr-package.R

@@ -0,0 +1,11 @@
+#' segregatr: Segregation Analysis for Identifying Pathogenic Variants
+#'
+#' Quantifying the evidence of co-segregation in the evaluation of putative
+#' pathogenic variants
+#'
+#' @docType package
+#' @import pedtools
+#' @importFrom pedprobr likelihood
+#'
+#' @name segregatr
+NULL

R/startdata_causative.R

@@ -0,0 +1,68 @@
+startdata_causative = function(x, marker, aff, penetrances, liabilityClasses = NULL) {
+
+  # Penetrances
+  if(is.null(liabilityClasses))
+    liabilityClasses = rep_len(1, pedsize(x))
+  if(is.null(dim(penetrances)))
+    penetrances = as.data.frame(as.list(penetrances))
+
+  if(!all(liabilityClasses %in% 1:nrow(penetrances)))
+    stop2("Illegal liability class: ", setdiff(liabilityClasses, 1:nrow(penetrances)))
+
+  glist = pedprobr:::.buildGenolist(x, marker, eliminate = 1)
+
+  if (attr(glist, "impossible")) {
+    dat = structure(list(), impossible = TRUE)
+    return(dat)
+  }
+
+  FOU = founders(x, internal = TRUE)
+
+  # Founder inbreeding: A vector of length pedsize(x), with NA's at nonfounders
+  # Enables quick look-up e.g. FOU_INB[i].
+  FOU_INB = rep(NA_real_, pedsize(x))
+  FOU_INB[FOU] = founderInbreeding(x)
+
+  # Allele frequencies (used in HW below)
+  afr = afreq(marker)
+
+  # Initialise `impossible` attribute
+  impossible = FALSE
+
+  dat = lapply(1:pedsize(x), function(i) {
+    h = glist[[i]]
+
+    # Penetrance values
+    liab = liabilityClasses[i]
+    penet = as.numeric(penetrances[liab, 1:3])
+
+    # Affection status priors
+    nmut = h[1, ] + h[2, ] - 2
+    affi = aff[i]
+    if(is.na(affi))
+      affpriors = rep_len(1, length(nmut))
+    else if(affi)
+      affpriors = penet[nmut + 1]
+    else
+      affpriors = 1 - penet[nmut + 1]
+
+    # Total genotype priors
+    prob = as.numeric(affpriors)
+    if (i %in% FOU)
+      prob = prob * pedprobr::HWprob(h[1, ], h[2, ], afr, f = FOU_INB[i])
+
+    # Remove impossible entries
+    zer = prob == 0
+    if (any(zer)) {
+      h = h[, !zer, drop = F]
+      prob = prob[!zer]
+      if (all(zer))
+        impossible = TRUE
+    }
+
+    list(hap = h, prob = prob)
+  })
+
+  attr(dat, "impossible") = impossible
+  dat
+}

R/utils.R

@@ -0,0 +1,21 @@
+stop2 = function(...) {
+  a = lapply(list(...), toString)
+  a = append(a, list(call. = FALSE))
+  do.call(stop, a)
+}
+
+# Test that input is a positive (or similar) integer.
+isCount = function(x, minimum = 1) {
+  isTRUE(length(x) == 1 &&
+           (is.integer(x) || (is.numeric(x) && x == as.integer(x))) &&
+           x >= minimum)
+}
+
+.mysetdiff = function(x, y) {
+  unique.default(x[match(x, y, 0L) == 0L])
+}
+
+# Fast intersection. NB: assumes no duplicates!
+.myintersect = function (x, y) {
+  y[match(x, y, 0L)]
+}

inst/extdata/RanolaPedigree549.ped

@@ -0,0 +1,23 @@
+id fid mid sex
+1 0 0 1
+2 0 0 2
+3 1 2 2
+4 1 2 1
+5 1 2 2
+6 1 2 1
+7 0 0 1
+8 7 3 1
+9 7 3 1
+10 7 3 2
+11 7 3 1
+12 7 3 2
+13 0 0 2
+14 4 13 1
+15 4 13 2
+16 0 0 1
+17 16 5 1
+18 16 5 1
+19 16 5 2
+20 0 0 2
+21 6 20 1
+22 6 20 2

inst/extdata/ThompsonPedigree2.ped

@@ -0,0 +1,9 @@
+id fid mid sex
+1 0 0 1
+2 0 0 2
+3 1 2 1
+4 1 2 2
+5 1 2 2
+6 0 0 2
+7 3 6 2
+8 3 6 1

segregatr.Rproj

@@ -0,0 +1,21 @@
+Version: 1.0
+
+RestoreWorkspace: No
+SaveWorkspace: No
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: knitr
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(segregatr)
+
+test_check("segregatr")