initial commit

.github/workflows/R-CMD-check.yaml

@@ -0,0 +1,50 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: R-CMD-check.yaml
+
+permissions: read-all
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: true
+          build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'

.gitignore

@@ -0,0 +1,49 @@
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+.RDataTmp
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+
+# pkgdown site
+docs/
+
+# translation temp files
+po/*~
+
+# RStudio Connect folder
+rsconnect/

DESCRIPTION

@@ -0,0 +1,25 @@
+Package: geneticTMT
+Version: 0.1.0
+Date: 2024-07-01
+Depends: R (>= 4.0.0)
+Imports: stats
+Suggests: testthat
+Title: Genetic Transmission Mean Test
+Authors@R: c(
+    person(given = "Yushi",
+	   family = "Tang",
+	   role = c("aut", "cre"),
+	   email = "yushi.tang@princeton.edu",
+           comment = c(ORCID = "0000-0002-3809-2129")),
+    person(given = "John D.",
+           family = "Storey",
+	   role = c("aut"),
+	   email = "jstorey@princeton.edu",
+           comment = c(ORCID = "0000-0001-5992-402X"))
+    )
+Maintainer: Yushi Tang <yushi.tang@princeton.edu>
+Description: Functions for inferring causal genotype-phenotype relationships for population-sampled nuclear families.
+License: GPL-3 + file LICENSE
+URL: https://github.com/StoreyLab/geneticTMT
+RoxygenNote: 7.3.1
+Encoding: UTF-8

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2024 Storey Lab
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

NAMESPACE

@@ -0,0 +1,13 @@
+# Generated by roxygen2: do not edit by hand
+
+export(assignment_index)
+export(calc_dtdt)
+export(calc_dtmt)
+export(calc_dtmt_nc)
+export(calc_dtmt_pc)
+export(calc_mu_hat)
+export(calc_varhat_dtmt)
+export(tdtest)
+export(tmtest)
+importFrom(stats,pchisq)
+importFrom(stats,pnorm)

R/assignment_index.R

@@ -0,0 +1,78 @@
+#' Assignment Index
+#'
+#' This function constructs assignment indices based on trio genotypes. 
+#'
+#' @param is_phased Label for unphased/phased genotype. If `FALSE` (default), the
+#'     child's genotype is unphased. Else, the child's genotype is phased.
+#' @param trio_geno A \eqn{3 \times n}{3 x n} matrix of trio genotype if `FALSE`
+#'     for \code{is_phased}. Or a \eqn{4 \times n}{4 x n} matrix of trio genotype if
+#'     `TRUE` for \code{is_phased}. Must contain only values \eqn{(0, 1, 2)}{0, 1, 2}.
+#'     The first row is assumed maternal, the second row is assumed paternal, and
+#'     the third (or together with the fourth) row is assumed child's genotypes.
+#'
+#' @return A named list of two vectors containing indicators for assigning to 
+#'     treatment group 0 or treatment group 1.
+#'     
+#' @note Undesired dimensions of the genotype matrices will result in errors.
+#' 
+#' @examples
+#' trio_example <- matrix(c(1,1,2,0, 1,2,1,0, 1,1,1,0), nrow=3, byrow=TRUE)
+#' list_assign <- assignment_index(trio_example, is_phased=FALSE)
+#' 
+#' @export
+assignment_index <- function(trio_geno, is_phased=FALSE){
+  if(is_phased){
+    # validate input
+    if(nrow(trio_geno)!=4) {
+      stop('undesired dimension for phased trio genotypes')
+    }
+    # assignment index
+    z_m <- trio_geno[1,]
+    z_p <- trio_geno[2,]
+    a_m <- trio_geno[3,]
+    a_p <- trio_geno[4,]
+    w1m <- as.numeric((z_m==1) & (a_m==1))
+    w0m <- as.numeric((z_m==1) & (a_m==0))
+    w1p <- as.numeric((z_p==1) & (a_p==1))
+    w0p <- as.numeric((z_p==1) & (a_p==0))
+    w1 <- w1m + w1p
+    w0 <- w0m + w0p
+    # output
+    list_assign <- list(w1, w0, w1m, w0m, w1p, w0p)
+    names(list_assign) <- c('W1','W0','W1m','W0m','W1p','W0p')
+    return(list_assign)
+  } else {
+    # validate input
+    if(nrow(trio_geno)!=3) {
+      stop('undesired dimension for unphased trio genotypes')
+    }
+    # assignment index
+    z_m <- trio_geno[1,]
+    z_p <- trio_geno[2,]
+    g <- trio_geno[3,]
+    # assign to control
+    w0 <- as.numeric((z_m==2) & (z_p==1) & (g==1)) +
+          as.numeric((z_m==1) & (z_p==2) & (g==1)) +
+          as.numeric((z_m==1) & (z_p==1) & (g==1)) +
+          2 * as.numeric((z_m==1) & (z_p==1) & (g==0)) +
+          as.numeric((z_m==1) & (z_p==0) & (g==0)) +
+          as.numeric((z_m==0) & (z_p==1) & (g==0)) 
+    # assign to treatment
+    w1 <- as.numeric((z_m==2) & (z_p==1) & (g==2)) +
+          as.numeric((z_m==1) & (z_p==2) & (g==2)) +
+          2 * as.numeric((z_m==1) & (z_p==1) & (g==2)) +
+          as.numeric((z_m==1) & (z_p==1) & (g==1)) +
+          as.numeric((z_m==1) & (z_p==0) & (g==1)) +
+          as.numeric((z_m==0) & (z_p==1) & (g==1))
+    # output
+    list_assign <- list(w1, w0)
+    names(list_assign) <- c('W1','W0')
+    return(list_assign)
+  }
+}
+
+
+
+
+
+

R/calc_dtdt.R

@@ -0,0 +1,43 @@
+#' Transmission-Disequilibrium Test Statistic
+#'
+#' This function calculates the TDT statistic given trio genotypes and children's 
+#' dichotomous phenotype.
+#'
+#'
+#' @param Y The n-length vector of the children's dichotomous trait
+#' @param is_phased Label for unphased/phased genotype. If `FALSE` (default), the
+#'     child's genotype is unphased. Else, the child's genotype is phased.
+#' @param trio_geno A \eqn{3 \times n}{3 x n} matrix of trio genotype if `FALSE`
+#'     for is_phased. Or a \eqn{4 \times n}{4 x n} matrix of trio genotype if
+#'     `TRUE` for is_phased. Must contain only values \eqn{(0, 1, 2)}{0, 1, 2}.
+#'     The first row is assumed maternal, the second row is assumed paternal, and
+#'     the third (or together with the fourth) row is assumed child's genotype.
+#'
+#' @return A numeric value of the TDT statistic.
+#' 
+#' @note Inputting phenotypes that are not dichotomous and undesired dimensions
+#'     of the genotype matrices will result in errors.
+#'
+#' @examples 
+#' trio_example <- matrix(c(1,1,2,0, 1,2,1,0, 1,1,1,0), nrow=3, byrow=TRUE) 
+#' y_binary <- c(0,1,1,0)
+#' dtdt <- calc_dtdt(y_binary, trio_example)
+#' 
+#' @export
+calc_dtdt <- function(Y, trio_geno, is_phased=FALSE){
+  if(!all(Y %in% 0:1)){stop("Y must be a binary vector")}
+  list_assign <- assignment_index(trio_geno, is_phased)
+  W0 <- list_assign$W0
+  W1 <- list_assign$W1
+  # sizes of treatment and control
+  N <- sum(W1) + sum(W0)
+  # assign trait value to treatment and control group
+  vec_n1 <- rep(Y, list_assign$W1)
+  vec_n0 <- rep(Y, list_assign$W0)
+  n1 <- sum(vec_n1==1)
+  n0 <- sum(vec_n0==1)
+  ttdt <- ifelse(N>0, 2*(n1-n0)/N, NaN)
+  return(ttdt)
+}
+
+

R/calc_dtmt.R

@@ -0,0 +1,44 @@
+# Transmission-Mean Test Statistic
+#'
+#' This function computes the unbiased TMT estimand given trio genotypes and 
+#' children's phenotype. The phenotype can be either quantitative or dichotomous.
+#'
+#' @param Y The n-length vector of children' phenotype.
+#' @param is_phased Label for unphased/phased genotype. If `FALSE` (default), the
+#'     child's genotype is unphased. Else, the child's genotype is phased.
+#' @param trio_geno A \eqn{3 \times n}{3 x n} matrix of trio genotype if `FALSE`
+#'     for is_phased. Or a \eqn{4 \times n}{4 x n} matrix of trio genotype if
+#'     `TRUE` for is_phased. Must contain only values \eqn{(0, 1, 2)}{0, 1, 2}.
+#'     The first row is assumed maternal, the second row is assumed paternal, and
+#'     the third (or together with the fourth) row is assumed child's genotypes.
+#'
+#' @return A numeric value of the TMT statistic.
+#' 
+#' @note Undesired dimensions of the genotype matrices will result in errors.
+#'
+#' @examples 
+#' trio_example <- matrix(c(1,1,2,0, 1,2,1,0, 1,1,1,0), nrow=3, byrow=TRUE) 
+#' y_quant <- c(11.8,1.7,3.6,12.2)
+#' dtmt <- calc_dtmt(y_quant, trio_example)
+#' 
+#' @export
+calc_dtmt <- function(Y, trio_geno, is_phased=FALSE){
+  # construct assignment vectors
+  list_assign <- assignment_index(trio_geno, is_phased)
+  W0 <- list_assign$W0
+  W1 <- list_assign$W1
+  # estimate mu to center trait value
+  mu_hat <- calc_mu_hat(Y, W0, W1)
+  # sizes of treatment and control
+  N <- sum(W1) + sum(W0)
+  # calculate the tmt estimand
+  dtmt <- ifelse(N>0, 2/N*(sum((Y-mu_hat)*W1) - sum((Y-mu_hat)*W0)), NaN)
+  return(dtmt)
+}
+
+
+
+
+
+
+

R/calc_dtmt_nc.R

@@ -0,0 +1,37 @@
+# Transmission-Mean Test Preliminary Statistic 
+#'
+#' This function computes the preliminary TMT statistic given trio genotypes and 
+#' children's phenotype. The phenotype can be either quantitative or dichotomous.
+#' The preliminary TMT statistic is unbiased and non-centered by the population
+#' mean.
+#'
+#' @param Y The n-length vector of children' phenotype.
+#' @param is_phased Label for unphased/phased genotype. If `FALSE` (default), the
+#'     child's genotype is unphased. Else, the child's genotype is phased.
+#' @param trio_geno A \eqn{3 \times n}{3 x n} matrix of trio genotype if `FALSE`
+#'     for is_phased. Or a \eqn{4 \times n}{4 x n} matrix of trio genotype if
+#'     `TRUE` for is_phased. Must contain only values \eqn{(0, 1, 2)}{0, 1, 2}.
+#'     The first row is assumed maternal, the second row is assumed paternal, and
+#'     the third (or together with the fourth) row is assumed child's genotypes.
+#'
+#' @return A numeric value of the preliminary TMT statistic.
+#' 
+#' @note Undesired dimensions of the genotype matrices will result in errors.
+#'
+#' @examples 
+#' trio_example <- matrix(c(1,1,2,0, 1,2,1,0, 1,1,1,0), nrow=3, byrow=TRUE) 
+#' y_quant <- c(11.8,1.7,3.6,12.2)
+#' dtmt_nc <- calc_dtmt_nc(y_quant, trio_example)
+#' 
+#' @export
+calc_dtmt_nc <- function(Y, trio_geno, is_phased=FALSE){
+  # construct assignment vectors
+  list_assign <- assignment_index(trio_geno, is_phased)
+  W0 <- list_assign$W0
+  W1 <- list_assign$W1
+  # sizes of treatment and control
+  N <- sum(W1) + sum(W0)
+  # calculate the tmt estimand
+  dtmt <- ifelse(N>0, 2/N*(sum(Y*W1) - sum(Y*W0)), NaN)
+  return(dtmt)
+}