update

DESCRIPTION

@@ -1,15 +1,18 @@
 Package: IMAGE
 Type: Package
-Title: What the Package Does Using Title Case
+Title: Title: High-Powered Detection of Genetic Effects on DNA
+        Methylation Using Integrated Methylation QTL Mapping and
+        Allele-Specific Analysis
 Version: 1.0
-Date: 2019-03-01
-Author: Your Name
-Maintainer: Your Name <[email protected]>
-Description: More details about what the package does. See
-        <http://cran.r-project.org/doc/manuals/r-release/R-exts.html#The-DESCRIPTION-file>
-        for details on how to write this part.
+Date: 2019-07-08
+Author: Yue Fan, Shiquan Sun, Xiang Zhou
+Maintainer: Yue Fan <[email protected]>
+Description: MQTL mapping in bisulfite sequencing studies by fitting a binomial mixed model, incorporating allelic-specific methylation pattern.
+Encoding: UTF-8
 License: GPL (>= 2)
+Depends: R (>= 3.6.0)
 Imports: Rcpp (>= 0.12.19),foreach,doParallel,parallel,Matrix
 LinkingTo: Rcpp, RcppArmadillo
+RoxygenNote: 6.1.1
 NeedsCompilation: yes
-Packaged: 2019-06-21 18:09:38 UTC; yuef
+Packaged: 2019-07-08 20:14:35 UTC; yuef

NAMESPACE

@@ -1,10 +1,19 @@
-useDynLib(IMAGE, .registration=TRUE)
+# Generated by roxygen2: do not edit by hand
+
 export(image)
-importFrom(Rcpp, evalCpp)
-importFrom("stats", "binomial", "glm", "model.frame", "model.matrix",
-             "na.omit", "na.pass", "pchisq", "poisson", "var")
-import(parallel)
-import(foreach)
-import(doParallel)
 import(Matrix)
-
+import(doParallel)
+import(foreach)
+import(parallel)
+importFrom("stats","binomial")
+importFrom("stats","glm")
+importFrom("stats","lm")
+importFrom("stats","model.frame")
+importFrom("stats","model.matrix")
+importFrom("stats","na.omit")
+importFrom("stats","na.pass")
+importFrom("stats","pchisq")
+importFrom("stats","poisson")
+importFrom("stats","var")
+importFrom(Rcpp,evalCpp)
+useDynLib(IMAGE, .registration=TRUE)

R/IMAGE.R

@@ -8,6 +8,36 @@
 #          University of Michigan, Department of Biostatistics
 ########################################################################################################################
 
+#' High-powered detection of genetic effects on DNA methylation using integrated methylation QTL mapping and allele-specific analysis.
+#' 
+#' mQTL mapping in bisulfite sequencing studies by fitting a binomial mixed model, incorporating allelic-specific methylation pattern.
+#' 
+#' IMAGE properly accounts for the count nature of bisulfite sequencing data and incorporates allele-specific methylation patterns from heterozygous individuals to enable more powerful mQTL discovery. 
+#' Binomial mixed models (BMM) are fitted using the penalized quasi-likelihood (PQL) method proposed by Breslow and Clayton (1993).
+#'
+#' @param geno a data list containing the genotype data.
+#' @param data a data list containing the methylation data.
+#' @param K a known kinship matrix. This matrix should be a positive semi-definite matrix with dimensions equal to the samplie size.
+#' @param Covariates a matrix containing the covariates subject to adjustment (Default = NULL).
+#' @param numCore a positive integer specifying the number of cores for parallel computing (default = 1).
+#' @param fit.maxiter a positive integer specifying the maximum number of iterations when fitting the generalized linear mixed model (default = 500).
+#' @param fit.tol a positive number specifying tolerance, the difference threshold for parameter estimates below which iterations should be stopped (default = 1e-5).
+#' @param verbose a logical switch for printing detailed information (parameter estimates in each iteration) for testing and debugging purpose (default = TRUE).
+#' @return \item{loc}{ordinal number of SNP-CpG pair being analyzed}
+#' @return \item{numIDV}{number of observations of SNP-CpG pair being analyzed}
+#' @return \item{beta}{the fixed effect parameter estimate for the predictor of interest.}
+#' @return \item{se_beta}{the standard deviation of fixed effect.}
+#' @return \item{pvalue}{P value for the fixed effect, based on the wald test.}
+#' @return \item{h2}{heritability of the transformed rate.}
+#' @return \item{sigma2}{total variance component.}
+#' @return \item{converged}{a logical indicator for convergence.}
+#' @references Breslow, N.E. and Clayton, D.G. (1993) Approximate Inference in Generalized Linear Mixed Models. Journal of the American Statistical Association 88, 9-25.
+#' @author Yue Fan, Shiquan Sun, Xiang Zhou
+#' @examples 
+#' data(exampledata)
+#' res=image(exampledata$geno,exampledata$data,exampledata$K)
+
+
 
 image <- function(geno,data,K,Covariates=NULL,numCore=1,fit.maxiter=500,fit.tol=1e-5,verbose=TRUE) {
 
@@ -25,7 +55,7 @@ image <- function(geno,data,K,Covariates=NULL,numCore=1,fit.maxiter=500,fit.tol=
     numCov     <- dim(Covariates)[2]
     Covariates <- as.matrix(Covariates)
   }
-
+  iVar=0
   resBMM <- foreach(iVar=1:m, .combine=rbind, .errorhandling = 'remove')%dopar%{
 
     res <- data.frame()
@@ -81,7 +111,7 @@ image <- function(geno,data,K,Covariates=NULL,numCore=1,fit.maxiter=500,fit.tol=
       n2 <- 0
       CountData <- c(data$y[homo, iVar])
       LibSize <- c(data$r[homo, iVar])
-      genotypes <- c(pheno1[homo])
+      genotypes <- c(geno1[homo])
       if(numCov>0)
       {
         covariates=Covariates[homo,]
@@ -258,7 +288,7 @@ data_modify <- function(genotypes,ratio,LibSize,CountData) {
 ##########################################################
 #           	   PQLseq FIT FUNCTION					 #
 ##########################################################
-#' @export
+
 PQLseq.fit <- function(model0, RelatednessMatrix, method = "REML", method.optim = "AI", maxiter = 500, tol = 1e-5, verbose = FALSE) {
 
   names(RelatednessMatrix) <- paste("kins", 1:length(RelatednessMatrix), sep="")
@@ -282,7 +312,7 @@ PQLseq.fit <- function(model0, RelatednessMatrix, method = "REML", method.optim
 ##########################################################
 #       PQLseq FIT AVERAGE INFORMATION FUNCTION			 #
 ##########################################################
-#' @export
+
 PQLseq.AI <- function(model0, RelatednessMatrix, tau = rep(0, length(RelatednessMatrix)+1), fixtau = rep(0, length(RelatednessMatrix)+1), maxiter = 500, tol = 1e-5, verbose = FALSE) {
 
   if(model0$family$family %in% c("binomial")){
@@ -475,15 +505,15 @@ PQLseq.AI <- function(model0, RelatednessMatrix, tau = rep(0, length(Relatedness
 ##########################################################
 #           	   PQLseq2 FIT FUNCTION					 #
 ##########################################################
-#' @export
+
 PQLseq2.fit <- function(model0, RelatednessMatrix, method = "REML", method.optim = "AI", maxiter = 500, tol = 1e-5, verbose = FALSE) {
 
   names(RelatednessMatrix) <- paste("kins", 1:length(RelatednessMatrix), sep="")
   # if((method.optim == "AI")&(!sum(model0$fitted.values<1e-5))) {
   if(method.optim == "AI") {
     fixtau.old <- rep(0, length(RelatednessMatrix)+1)
     # to use average information method to fit alternative model
-    model1 <- PQLseq.AI2(model0, RelatednessMatrix, maxiter = maxiter, tol = tol, verbose = verbose)
+    model1 <- PQLseq2.AI(model0, RelatednessMatrix, maxiter = maxiter, tol = tol, verbose = verbose)
     fixtau.new <- 1*(model1$theta < 1.01 * tol)
 
 
@@ -497,7 +527,7 @@ PQLseq2.fit <- function(model0, RelatednessMatrix, method = "REML", method.optim
     while(any(fixtau.new != fixtau.old)) {
       fixtau.old <- fixtau.new
       # to use average information method to fit alternative model
-      model1 <- PQLseq.AI2(model0, RelatednessMatrix, fixtau = fixtau.old, maxiter = maxiter, tol = tol, verbose = verbose)
+      model1 <- PQLseq2.AI(model0, RelatednessMatrix, fixtau = fixtau.old, maxiter = maxiter, tol = tol, verbose = verbose)
 
 
 
@@ -513,7 +543,7 @@ PQLseq2.fit <- function(model0, RelatednessMatrix, method = "REML", method.optim
 ##########################################################
 #       PQLseq2 FIT AVERAGE INFORMATION FUNCTION			 #
 ##########################################################
-#' @export
+
 PQLseq2.AI <- function(model0, RelatednessMatrix, tau = rep(0, length(RelatednessMatrix)+1), fixtau = rep(0, length(RelatednessMatrix)+1), maxiter = 500, tol = 1e-5, verbose = FALSE) {
 
   n=nrow(RelatednessMatrix[[1]])

R/IMAGE_package.R

@@ -0,0 +1,10 @@
+#' @useDynLib IMAGE, .registration=TRUE
+#' @importFrom Rcpp evalCpp
+#' @importFrom "stats" "binomial" "glm" "model.frame" "model.matrix" "na.omit" "na.pass" "pchisq" "poisson" "var" "lm"
+#' @import parallel
+#' @import foreach
+#' @import doParallel
+#' @import Matrix
+#' @export image
+NULL
+

R/exampledata.R

@@ -0,0 +1,15 @@
+#' Example dataset
+#'
+#' A simulated example dataset of IMAGE.
+#' The variables are as follows:
+#'
+#' @format Contains the following objects:
+#' \describe{
+#'   \item{geno}{a data list containing 2 haplotypes.}
+#'   \item{data}{a data list containing the methylated read counts and total read counts for each allele.}
+#'   \item{K}{a genetic relationship matrix.}
+#' }
+#' @examples 
+#' data(exampledata)
+#' res=image(exampledata$geno,exampledata$data,exampledata$K)
+"exampledata"

man/IMAGE-package.Rd

@@ -3,30 +3,13 @@
 \alias{IMAGE}
 \docType{package}
 \title{\packageTitle{IMAGE}}
-\description{\packageDescription{IMAGE}}
-\details{
-  The DESCRIPTION file: \packageDESCRIPTION{IMAGE}
-  \packageIndices{IMAGE}
+\description{\packageDescription{IMAGE}mQTL mapping in bisulfite sequencing studies by fitting a binomial mixed model, incorporating allelic-specific methylation pattern.}
 
-  This section should provide a more detailed overview of how to use the
-  package, including the most important functions.
-}
 \author{
-  \packageAuthor{IMAGE}
-
-  Maintainer: \packageMaintainer{IMAGE}
-}
-\references{
-  This optional section can contain literature or other references for
-  background information.
+  Yue Fan, Shiquan Sun,  Xiang Zhou
+  Maintainer: Yue Fan <yuef@umich.edu>
 }
-% Optionally other standard keywords, one per line,
-% from the file KEYWORDS in the R documentation.
+
 \keyword{package}
-\seealso{
-  Optional links to other man pages
-}
-\examples{
-  ## Optional simple examples of the most important functions
-  ## Use \dontrun{} around code to be shown but not executed
-}
+
+