Add MZ twin support to adjacency matrix via twinID column

R/buildComponent.R

@@ -19,7 +19,7 @@
 #' @param isChild_method character. The method to use for computing the isChild matrix.  Options are "classic" or "partialparent"
 #' @param adjBeta_method numeric The method to use for computing the building the adjacency_method matrix when using the "beta" build
 #' @param compress logical. If TRUE, use compression when saving the checkpoint files.  Defaults to TRUE.
-#' @param mz_twins logical. If TRUE, treat MZ twins as having an additional parent-child relationship for the purposes of computing the relatedness matrix. Defaults to FALSE.
+#' @param mz_twins logical. If TRUE, merge MZ co-twin columns in the r2 matrix before tcrossprod so that MZ twins are coded with relatedness 1 instead of 0.5. Twin pairs are identified from the \code{twinID} column. When a \code{zygosity} column is also present, only pairs where both members have \code{zygosity == "MZ"} are used; otherwise all \code{twinID} pairs are assumed to be MZ. Defaults to FALSE.
 #' @param ... additional arguments to be passed to \code{\link{ped2com}}
 #' @details The algorithms and methodologies used in this function are further discussed and exemplified in the vignette titled "examplePedigreeFunctions". For more advanced scenarios and detailed explanations, consult this vignette.
 #' @export
@@ -123,9 +123,9 @@ ped2com <- function(ped, component,
     ped <- standardizeColnames(ped, verbose = config$verbose)
   }
 
+  mz_pairs <- NULL
   if (mz_twins == TRUE && "twinID" %in% colnames(ped)) {
-    # TODO
-    # ped <- addMZtwins(ped, verbose = config$verbose)
+    mz_pairs <- findMZtwins(ped, verbose = config$verbose)
   }
 
 
@@ -289,6 +289,22 @@ ped2com <- function(ped, component,
     compress = config$compress
   )
 
+  # --- Step 3b: Temporarily merge MZ twin columns in r2 ---
+  # MZ twins share the same genetic source.  We absorb twin2's column into
+  # twin1's before tcrossprod so all path-traced relatedness flows through a
+  # single source.  After tcrossprod we copy twin1's row/col back to twin2.
+  if (!is.null(mz_pairs) && length(mz_pairs) > 0) {
+    for (pair in mz_pairs) {
+      idx1 <- pair[1]
+      idx2 <- pair[2]
+      r2[, idx1] <- r2[, idx1] + r2[, idx2]
+      r2[, idx2] <- 0
+    }
+    if (config$verbose == TRUE) {
+      message("Merged ", length(mz_pairs), " MZ twin pair column(s) in r2")
+    }
+  }
+
   # --- Step 4: T crossproduct  ---
 
   if (config$resume == TRUE && file.exists(checkpoint_files$tcrossprod_checkpoint) &&
@@ -308,6 +324,20 @@ ped2com <- function(ped, component,
     }
   }
 
+  # --- Step 4b: Restore MZ twins ---
+  # Copy twin1's row/col to twin2 so both twins appear in the final matrix.
+  if (!is.null(mz_pairs) && length(mz_pairs) > 0) {
+    for (pair in mz_pairs) {
+      idx1 <- pair[1]
+      idx2 <- pair[2]
+      r[idx2, ] <- r[idx1, ]
+      r[, idx2] <- r[, idx1]
+    }
+    if (config$verbose == TRUE) {
+      message("Restored ", length(mz_pairs), " MZ twin pair(s) in relatedness matrix")
+    }
+  }
+
   if (config$component %in% c("mitochondrial", "mtdna", "mitochondria")) {
     r@x <- rep(1, length(r@x))
     # Assign 1 to all nonzero elements for mitochondrial component
@@ -343,6 +373,7 @@ ped2add <- function(ped, max_gen = 25, sparse = TRUE, verbose = FALSE,
                     save_rate_parlist = 100000 * save_rate,
                     save_path = "checkpoint/",
                     compress = TRUE,
+                    mz_twins = FALSE,
                     ...) {
   ped2com(
     ped = ped,
@@ -361,6 +392,7 @@ ped2add <- function(ped, max_gen = 25, sparse = TRUE, verbose = FALSE,
     save_rate_parlist = save_rate_parlist,
     save_path = save_path,
     compress = compress,
+    mz_twins = mz_twins,
     ...
   )
 }

R/constructAdjacency.R

@@ -559,3 +559,70 @@ isChild <- function(isChild_method, ped) {
     })
   }
 }
+
+
+#' Find MZ twin pairs in a pedigree
+#'
+#' Identifies MZ twin pairs from the \code{twinID} column and returns their
+#' row indices.  These indices are used later to merge the twins' columns in
+#' the \code{r2} matrix before \code{tcrossprod}, which correctly produces
+#' relatedness 1 between MZ co-twins with no diagonal or downstream artifacts.
+#'
+#' @param ped A pedigree data.frame with columns \code{ID} and \code{twinID}.
+#'   Optionally a \code{zygosity} column; when present only pairs where both
+#'   members have \code{zygosity == "MZ"} are used.
+#' @param verbose logical. If TRUE, print progress messages.
+#' @return A list of length-2 integer vectors \code{c(idx1, idx2)} giving the
+#'   row indices of each MZ pair in the pedigree, or \code{NULL} if none found.
+#' @keywords internal
+findMZtwins <- function(ped, verbose = FALSE) {
+  if (!"twinID" %in% colnames(ped)) {
+    return(NULL)
+  }
+
+  twin_rows <- which(!is.na(ped$twinID))
+
+  # If zygosity column exists, restrict to MZ pairs
+  if ("zygosity" %in% colnames(ped)) {
+    twin_rows <- twin_rows[!is.na(ped$zygosity[twin_rows]) &
+      ped$zygosity[twin_rows] %in% c("mz", "MZ")]
+  }
+
+  if (length(twin_rows) == 0) {
+    return(NULL)
+  }
+
+  processed <- c()
+  pairs <- list()
+
+  for (idx in twin_rows) {
+    twin_id <- ped$ID[idx]
+    co_twin_id <- ped$twinID[idx]
+
+    # Skip if already processed this pair
+    if (twin_id %in% processed || co_twin_id %in% processed) next
+
+    idx1 <- which(ped$ID == twin_id)
+    idx2 <- which(ped$ID == co_twin_id)
+
+    if (length(idx1) != 1 || length(idx2) != 1) next
+
+    # Always put the lower index first for consistency
+    if (idx1 > idx2) {
+      tmp <- idx1
+      idx1 <- idx2
+      idx2 <- tmp
+    }
+
+    processed <- c(processed, twin_id, co_twin_id)
+    pairs[[length(pairs) + 1]] <- c(idx1, idx2)
+
+    if (verbose) {
+      message("MZ twin pair found: ", twin_id, " (row ", idx1,
+        ") and ", co_twin_id, " (row ", idx2, ")")
+    }
+  }
+
+  if (length(pairs) == 0) return(NULL)
+  return(pairs)
+}