Ibd partitioning

R/AlphaPart.R

@@ -74,7 +74,11 @@
 #' @param colPath Numeric or character, position or name of a column
 #'   holding path information.
 #' @param colBV Numeric or character, position(s) or name(s) of
-#'   column(s) holding genetic Values.
+#'   column(s) holding genetic values.
+#' @param colPaternalBV Numeric or character, position(s) of column(s) 
+#'   holding the paternal genetic values for IBD partitioning.
+#' @param colMaternalBV Numeric or character, position(s) of column(s) 
+#'   holding the maternal genetic values for IBD partitioning..
 #' @param colBy Numeric or character, position or name of a column
 #'   holding group information (see details).
 #'
@@ -164,6 +168,8 @@ AlphaPart <- function(
   colMid = 3,
   colPath = 4,
   colBV = 5:ncol(x),
+  colPaternalBV = NULL,
+  colMaternalBV = NULL,
   colBy = NULL
 ) {
   ## Test if the data is a data.frame
@@ -181,6 +187,15 @@ AlphaPart <- function(
       "arguments 'colId', 'colFid', 'colMid', 'colPath', and 'colBy' must be of length 1"
     )
   }
+
+  if (!is.null(colPaternalBV) | !is.null(colMaternalBV)){
+    if (is.null(colPaternalBV) | is.null(colMaternalBV)){
+      stop("Both 'colPaternalBV' and 'colMaternalBV' must be provided for IBD partitioning.")
+    }
+    gameticPartition <- TRUE
+  } else {
+    gameticPartition <- FALSE
+  }
 
   if (is.null(colBy)) {
     groupSummary <- FALSE
@@ -237,6 +252,24 @@ AlphaPart <- function(
     }
     testN <- NULL # not needed anymore
   }
+
+  if (gameticPartition){
+    if(!is.numeric(colPaternalBV)){
+      testN <- length(colPaternalBV)
+      colPaternalBV <- which(colnames(x) %in% colPaternalBV)
+      if (length(colPaternalBV) != testN) {
+        stop("Identification not valid for 'colPaternalBV' column(s) name", call. = FALSE)
+      }
+    }
+    if(!is.numeric(colMaternalBV)){
+      testN <- length(colMaternalBV)
+      colMaternalBV <- which(colnames(x) %in% colMaternalBV)
+      if (length(colMaternalBV) != testN) {
+        stop("Identification not valid for 'colMaternalBV' column(s) name", call. = FALSE)
+      }
+    }
+    testN <- NULL # not needed anymore
+  }
 
   ## --- Sort and recode pedigree ---
   ## Make sure that identifications are numeric if  recode=FALSE
@@ -253,6 +286,35 @@ AlphaPart <- function(
     stop("colBV columns must be numeric!")
     str(x)
   }
+
+  ## If gametic partitioning, make sure that colPaternalBV and colMaternalBV:
+  ## - are numeric
+  ## - the same length as colBV
+  ## - the columns sum to the colBV columns
+  if (gameticPartition){
+    test <- !sapply(x[, c(colPaternalBV)], is.numeric)
+    if (any(test)){
+      stop("colPaternalBV columns must be numeric!")
+      str(x)
+    }
+    test <- !sapply(x[, c(colMaternalBV)], is.numeric)
+    if (any(test)){
+      stop("colMaternalBV columns must be numeric!")
+      str(x)
+    }
+    test <- length(colPaternalBV) != length(colMaternalBV)
+    if (any(test)){
+      stop("colPaternalBV and colMaternalBV must be of the same length!")
+    }
+    test <- length(colPaternalBV) != length(colBV)
+    if (any(test)){
+      stop("colPaternalBV and colMaternalBV must be of the same length as colBV!")
+    }
+    test <- any(abs(x[, c(colBV)] - x[, c(colPaternalBV)] - x[, c(colMaternalBV)]) > 1e-8)
+    if (any(test)){
+      stop("The sum of colPaternalBV and colMaternalBV must be equal to colBV for each individual!")
+    }
+  }
 
   ## Sort so that parents precede children
   if (sort) {
@@ -307,7 +369,14 @@ AlphaPart <- function(
       }
     }
   }
-  y <- cbind(y, as.matrix(x[, colBV]))
+
+  if (!gameticPartition){
+    y <- cbind(y, as.matrix(x[, colBV]))
+    nGP <- 1 # Number of genetic partitions: total
+  } else {
+    y <- cbind(y, as.matrix(x[, colBV]), as.matrix(x[, colPaternalBV]), as.matrix(x[, colMaternalBV]))
+    nGP <- 3 # Number of genetic partitions: total, paternal, maternal
+  }
 
   ## Test if father and mother codes precede children code -
   ## computational engine needs this
@@ -430,6 +499,7 @@ AlphaPart <- function(
       nI = nI,
       nP = nP,
       nT = nT,
+      nGP = nGP,
       ped = y,
       P = as.integer(P),
       Px = as.integer(cumsum(c(0, rep(nP, nT - 1))))
@@ -448,14 +518,24 @@ AlphaPart <- function(
       nG = nG,
       ped = y,
       P = as.integer(P),
-      Px = as.integer(cumsum(c(0, rep(nP, nT - 1))))
+      Px = as.integer(cumsum(c(0, rep(nP, nT - 1)))),
+      g = g
     )
   }
 
   ## Assign nice column names
-  colnames(tmp$pa) <- paste(lT, "_pa", sep = "")
-  colnames(tmp$ms) <- paste(lT, "_ms", sep = "")
-  colnames(tmp$xa) <- c(t(outer(lT, lP, paste, sep = "_")))
+  if (!gameticPartition){
+    colnames(tmp$pa) <- paste(lT, "_pa", sep = "")
+    colnames(tmp$ms) <- paste(lT, "_ms", sep = "")
+    colnames(tmp$xa) <- c(t(outer(lT, lP, paste, sep = "_"))) 
+  } else {
+    #lPT <- colnames(x[, c(colBV, colPaternalBV, colMaternalBV), drop=FALSE])
+    lP <- c(lP, paste0(lP, "_paternal"), paste0(lP, "_maternal"))
+    lPT <- unlist(lapply(c("","_paternal", "_maternal"), function(s) paste0(lT, s)))
+    colnames(tmp$pa) <- paste(lPT, "_pa", sep="")
+    colnames(tmp$ms)  <- paste(lPT, "_ms", sep="")
+    colnames(tmp$xa) <- c(t(outer(lT, lP, paste, sep = "_")))
+  }
 
   ## --- Massage results ---
 
@@ -466,34 +546,64 @@ AlphaPart <- function(
   colP <- colnames(tmp$pa)
   colM <- colnames(tmp$ms)
   colX <- colnames(tmp$xa)
-
-  for (j in 1:nT) {
-    ## j <- 1
-    Py <- seq(t + 1, t + nP)
-    ret[[j]] <- cbind(tmp$pa[-1, j], tmp$ms[-1, j], tmp$xa[-1, Py])
-    colnames(ret[[j]]) <- c(colP[j], colM[j], colX[Py])
-    t <- max(Py)
+
+  if (!gameticPartition) {
+    for (j in 1:nT) {
+      ## j <- 1
+      Py <- seq(t + 1, t + nP)
+      ret[[j]] <- cbind(tmp$pa[-1, j], tmp$ms[-1, j], tmp$xa[-1, Py])
+      colnames(ret[[j]]) <- c(colP[j], colM[j], colX[Py])
+      t <- max(Py)
+    }
+  } else {
+    t <- 0
+    for (j in 1:nT) {
+      ## j <- 1
+      Py <- seq(t+1, j*nGP*nP)
+      pams <- c(j, j+nT, j+2*nT)
+      ret[[j]] <- cbind(tmp$pa[-1, pams], tmp$ms[-1, pams], tmp$xa[-1, Py])
+      colnames(ret[[j]]) <- c(colP[pams], colM[pams], colX[Py])
+      t <- max(Py)
+    }
   }
-  tmp <- NULL # not needed anymore
+
+  tmp <- NULL # not needed anymore 
 
   ## Add initial data
 
   if (!groupSummary) {
-    for (i in 1:nT) {
-      ## Hassle in order to get all columns and to be able to work with
-      ##   numeric or character column "names"
-      colX <- colX2 <- colnames(x)
-      names(colX) <- colX
-      names(colX2) <- colX2
-      ## ... put current agv in the last column in original data
-      colX <- c(colX[!(colX %in% colX[colBV[i]])], colX[colBV[i]])
-      ## ... remove other traits
-      colX <- colX[
-        !(colX %in%
-          colX2[(colX2 %in% colX2[colBV]) & !(colX2 %in% colX2[colBV[i]])])
-      ]
-      ret[[i]] <- cbind(x[, colX], as.data.frame(ret[[i]]))
-      rownames(ret[[i]]) <- NULL
+    if (!gameticPartition) {
+      for (i in 1:nT) {
+        ## Hassle in order to get all columns and to be able to work with
+        ##   numeric or character column "names"
+        colX <- colX2 <- colnames(x)
+        names(colX) <- colX
+        names(colX2) <- colX2
+        ## ... put current agv in the last column in original data
+        colX <- c(colX[!(colX %in% colX[colBV[i]])], colX[colBV[i]])
+        ## ... remove other traits
+        colX <- colX[
+          !(colX %in%
+              colX2[(colX2 %in% colX2[colBV]) & !(colX2 %in% colX2[colBV[i]])])
+        ]
+        ret[[i]] <- cbind(x[, colX], as.data.frame(ret[[i]]))
+        rownames(ret[[i]]) <- NULL
+      }
+    } else {
+      for (i in 1:nT) {
+        ## Hassle in order to get all columns and to be able to work with
+        ##   numeric or character column "names"
+        colX <- colX2 <- colnames(x)
+        names(colX) <- colX; names(colX2) <- colX2
+        ## ... put current agv in the last three columns in original data
+        colX <- c(colX[!(colX %in% c(colX[colBV[i]], colX[colPaternalBV[i]], colX[colMaternalBV[i]]))],
+                  colX[colBV[i]], colX[colPaternalBV[i]], colX[colMaternalBV[i]])
+        ## ... remove other traits
+        colX <- colX[!(colX %in% colX2[(colX2 %in% c(colX2[colBV], colX2[colPaternalBV], colX2[colMaternalBV])) & !
+                                         (colX2 %in% c(colX2[colBV[i]], colX2[colPaternalBV[i]], colX2[colMaternalBV[i]]))])]
+        ret[[i]] <- cbind(x[, colX], as.data.frame(ret[[i]]))
+        rownames(ret[[i]]) <- NULL
+      }
     }
   }
 
@@ -507,6 +617,7 @@ AlphaPart <- function(
     lP = lP,
     nT = nT,
     lT = lT,
+    gameticPartition = gameticPartition,
     warn = NULL
   )
   ## names(ret)[nT+1] <- "info"

R/RcppExports.R

@@ -1,8 +1,8 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-AlphaPartDrop <- function(c1, c2, nI, nP, nT, ped, P, Px) {
-    .Call(`_AlphaPart_AlphaPartDrop`, c1, c2, nI, nP, nT, ped, P, Px)
+AlphaPartDrop <- function(c1, c2, nI, nP, nT, nGP, ped, P, Px) {
+    .Call(`_AlphaPart_AlphaPartDrop`, c1, c2, nI, nP, nT, nGP, ped, P, Px)
 }
 
 AlphaPartDropGroup <- function(c1, c2, nI, nP, nT, nG, ped, P, Px, g) {

R/methods.R

@@ -262,6 +262,7 @@ summary.AlphaPart <- function(
     nT = nT,
     lT = lT,
     by = by,
+    gameticPartition = object$info$gameticPartition,
     warn = object$info$warn,
     labelSum = labelSum
   )
@@ -526,6 +527,8 @@ plot.summaryAlphaPart <-
     nP <- x$info$nP + x$info$nCov
     ret <- vector(mode = "list", length = nT)
     names(ret) <- x$info$lT
+
+    ifelse(x$info$gameticPartition, nGP <- 3, nGP <- 1)
 
     ## Axis labels
     if (!is.null(xlab) && length(xlab) > 1)
@@ -613,8 +616,8 @@ plot.summaryAlphaPart <-
     }
     ## Line type
     if (is.null(lineTypeList)) {
-      if (length(lineType) < nP) {
-        lineType <- c(1, rep(x = lineType, times = nP))
+      if (length(lineType) < nP*nGP) {
+        lineType <- c(1, rep(x = lineType, times = nP*nGP))
       } else {
         lineType <- c(1, lineType)
       }
@@ -643,7 +646,7 @@ plot.summaryAlphaPart <-
         if (sortValue) {
           nC <- ncol(tmp0)
           pathStat <- sapply(
-            X = tmp0[, (nC - nP + 1):nC],
+            X = tmp0[, (nC - (nP*nGP) + 1):nC],
             FUN = sortValueFUN,
             na.rm = TRUE
           )

src/AlphaPartDrop.cpp

@@ -2,17 +2,17 @@
 using namespace Rcpp;
 
 // [[Rcpp::export]]
-List AlphaPartDrop(double c1, double c2, int nI, int nP, int nT,
+List AlphaPartDrop(double c1, double c2, int nI, int nP, int nT, int nGP,
                    NumericMatrix ped, IntegerVector P, IntegerVector Px) {
   // --- Temp ---
 
-  int i, j, t, p;
+  int i, j, t, p, pt=0, mt=0, k;
 
   // --- Outputs ---
 
-  NumericMatrix pa(nI+1, nT);    // Parent average
-  NumericMatrix ms(nI+1, nT);    // Mendelian sampling
-  NumericMatrix xa(nI+1, nP*nT); // Partitions
+  NumericMatrix pa(nI+1, nT*nGP);    // Parent average
+  NumericMatrix ms(nI+1, nT*nGP);    // Mendelian sampling
+  NumericMatrix xa(nI+1, nP*nT*nGP); // Partitions
   // NOTE: Rcpp::NumericMatrix is filled by 0s by default
 
   // TODO: Maybe we want an algorithm that works on one trait at a time to save on memory?
@@ -22,7 +22,6 @@ List AlphaPartDrop(double c1, double c2, int nI, int nP, int nT,
   //       https://github.com/AlphaGenes/AlphaPart/issues/13
 
   // --- Compute ---
-
   for(i = 1; i < nI+1; i++) {
     for(t = 0; t < nT; t++) {
       // Parent average (PA)
@@ -31,18 +30,41 @@ List AlphaPartDrop(double c1, double c2, int nI, int nP, int nT,
 
       // Mendelian sampling (MS)
       ms(i, t) = ped(i, 3+t) - pa(i, t);
+
+      // Gametic partition
+      if (nGP ==3) {
+        pt = t + nT; // paternal trait index
+        mt = t + 2*nT; // maternal trait index
+        pa(i, pt) = c1 * ped(ped(i,1), 3+pt) + c1 * ped(ped(i,1), 3+mt);
+        pa(i, mt) = c2 * ped(ped(i,2), 3+pt) + c2 * ped(ped(i,2), 3+mt);
+        ms(i, mt) = ped(i, 3+mt) - pa(i, mt);
+        ms(i, pt) = ped(i, 3+pt) - pa(i, pt);
+      }
 
       // Parts
 
       // ... for the MS part
-      j = Px[t] + P[i];
+      j = nGP*Px[t] + P[i];
       xa(i, j) = ms(i, t);
+
+      if (nGP == 3) {
+        j = nGP*Px[t] + P[i] + nP;
+        xa(i, j) = ms(i, pt);
+        j = nGP*Px[t] + P[i] + nP*2;
+        xa(i, j) = ms(i, mt);
+      }
 
       // ... for the PA parts
       for(p = 0; p < nP; p++) {
-        j = Px[t] + p;
+        j = nGP*Px[t] + p;
         xa(i, j) += c1 * xa(ped(i, 1), j) +
                     c2 * xa(ped(i, 2), j);
+        if (nGP == 3) {
+          j = nGP*Px[t] + p + nP;
+          k = nGP*Px[t] + p + nP*2;
+          xa(i, j) += c1 * xa(ped(i, 1), j) + c1 * xa(ped(i, 1), k);
+          xa(i, k) += c2 * xa(ped(i, 2), j) + c2 * xa(ped(i, 2), k);
+        }
       }
     }
   }