## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("coMethDMR") ## ----load_comethDMR, message=FALSE-------------------------------------------- library(coMethDMR) ## ----------------------------------------------------------------------------- data(betasChr22_df) betasChr22_df[1:5, 1:5] ## ----------------------------------------------------------------------------- data(pheno_df) head(pheno_df) ## ----------------------------------------------------------------------------- CpGisland_ls <- readRDS( system.file( "extdata", "CpGislandsChr22_ex.rds", package = 'coMethDMR', mustWork = TRUE ) ) ## ----CoMethAllRegions--------------------------------------------------------- system.time( coMeth_ls <- CoMethAllRegions( dnam = betasChr22_df, betaToM = TRUE, # converts beta to m-values method = "pearson", CpGs_ls = CpGisland_ls, arrayType = "450k", returnAllCpGs = FALSE, output = "CpGs", nCores_int = 1 ) ) # ~15 seconds coMeth_ls ## ----message=FALSE------------------------------------------------------------ WriteCorrPlot <- function(beta_mat){ require(corrplot) require(coMethDMR) CpGs_char <- row.names(beta_mat) CpGsOrd_df <- OrderCpGsByLocation( CpGs_char, arrayType = c("450k"), output = "dataframe" ) betaOrdered_mat <- t(beta_mat[CpGsOrd_df$cpg ,]) corr <- cor( betaOrdered_mat, method = "spearman", use = "pairwise.complete.obs" ) corrplot(corr, method = "number", number.cex = 1, tl.cex = 0.7) } # subsetting beta values to include only co-methylated probes betas_df <- subset( betasChr22_df, row.names(betasChr22_df) %in% coMeth_ls[[1]] ) WriteCorrPlot(betas_df) ## ----lmmTestAllRegions_1------------------------------------------------------ out_df <- lmmTestAllRegions( betas = betasChr22_df, region_ls = coMeth_ls, pheno_df, contPheno_char = "stage", covariates_char = NULL, modelType = "randCoef", arrayType = "450k" # generates a log file in the current directory # outLogFile = paste0("lmmLog_", Sys.Date(), ".txt") ) out_df ## ----AnnotateResults, message = FALSE----------------------------------------- system.time( outAnno_df <- AnnotateResults( lmmRes_df = out_df, arrayType = "450k" ) ) # ~12 seconds outAnno_df ## ----CpGsInfoOneRegion-------------------------------------------------------- outCpGs_df <- CpGsInfoOneRegion( regionName_char = "chr22:18268062-18268249", betas_df = betasChr22_df, pheno_df = pheno_df, contPheno_char = "stage", covariates_char = NULL, arrayType = "450k" ) outCpGs_df ## ----message = FALSE---------------------------------------------------------- # library("GenoGAM") # NOTE 2022-03-28: what are we using this package for again? ## ----CpGsInfoOneRegion_2------------------------------------------------------ CpGsInfoOneRegion( regionName_char = "chr22:19709548-19709755", betas_df = betasChr22_df, pheno_df = pheno_df, contPheno_char = "stage", covariates_char = NULL, arrayType = "450k" ) ## ----------------------------------------------------------------------------- gene_ls <- readRDS( system.file( "extdata", "450k_Gene_3_200.rds", package = 'coMethDMR', mustWork = TRUE ) ) ## ----reduce_data-------------------------------------------------------------- Cgi_ls <- readRDS( system.file( "extdata", "CpGislandsChr22_ex.rds", package = 'coMethDMR', mustWork = TRUE ) ) betasChr22_df <- betasChr22_df[rownames(betasChr22_df) %in% unlist(Cgi_ls)[1:20], ] ## ----GetResiduals------------------------------------------------------------- resid_mat <- GetResiduals( dnam = betasChr22_df, # converts to Mvalues for fitting linear model betaToM = TRUE, pheno_df = pheno_df, covariates_char = c("age.brain", "sex", "slide") ) ## ----CoMethAllRegions_2------------------------------------------------------- Cgi_ls <- readRDS( system.file( "extdata", "CpGislandsChr22_ex.rds", package = 'coMethDMR', mustWork = TRUE ) ) coMeth_ls <- CoMethAllRegions( dnam = resid_mat, betaToM = FALSE, method = "pearson", CpGs_ls = Cgi_ls[1], arrayType = "450k", returnAllCpGs = FALSE, output = "CpGs" ) coMeth_ls ## ----CoMethAllRegions_3------------------------------------------------------- coMeth_ls <- CoMethAllRegions( dnam = resid_mat, betaToM = FALSE, CpGs_ls = Cgi_ls[5], arrayType = "450k", returnAllCpGs = FALSE, output = "CpGs" ) coMeth_ls ## ----CoMethAllRegions_4------------------------------------------------------- coMethData_df <- CoMethAllRegions( dnam = resid_mat, betaToM = FALSE, CpGs_ls = Cgi_ls[5], arrayType = "450k", returnAllCpGs = FALSE, output = "dataframe" ) [[1]] coMethData_df ## ----lmmTestAllRegions_2, message = FALSE------------------------------------- lmmTestAllRegions( betas = betasChr22_df, region_ls = coMeth_ls[1], pheno_df, contPheno_char = "stage", covariates_char = "age.brain", modelType = "randCoef", arrayType = "450k" ) ## ----------------------------------------------------------------------------- data(betasChr22_df) ## ----message = FALSE---------------------------------------------------------- # list probes for this gene ARFGAP3_CpGs_char <- c( "cg00079563", "cg01029450", "cg02351223", "cg04527868", "cg09861871", "cg26529516", "cg00539564", "cg05288033", "cg09367092", "cg10648908", "cg14570855", "cg15656623", "cg23778094", "cg27120833" ) # list probes located closely on this gene gene3_200 <- CloseBySingleRegion( CpGs_char = ARFGAP3_CpGs_char, arrayType = "450k", maxGap = 200, minCpGs = 3 ) CpGsOrdered_ls <- lapply( gene3_200, OrderCpGsByLocation, arrayType = "450k", output = "dataframe" ) names(gene3_200) <- lapply(CpGsOrdered_ls, NameRegion) # List of regions gene3_200 ## ----find_and_test_regions---------------------------------------------------- # co-methlyated region within the gene coMeth_ls <- CoMethAllRegions( dnam = betasChr22_df, betaToM = TRUE, method = "pearson", CpGs_ls = gene3_200, arrayType = "450k", returnAllCpGs = FALSE, output = "CpGs" ) coMeth_ls # test the co-methylated regions within the gene results <- lmmTestAllRegions( betas = betasChr22_df, region_ls = coMeth_ls, pheno_df, contPheno_char = "stage", covariates_char = "age.brain", modelType = "randCoef", arrayType = "450k" # generates a log file in the current directory # outLogFile = paste0("lmmLog_", Sys.Date(), ".txt") ) ## ----------------------------------------------------------------------------- AnnotateResults(lmmRes_df = results, arrayType = "450k") ## ----------------------------------------------------------------------------- utils::sessionInfo()