library(dplyr)cohort <- "MtSinai"
data.dir <- file.path("DATASETS/",cohort,"/")
data.dir.table <- "DATASETS/Summary_Table/"
data.dir.raw <- file.path(data.dir,"/step1_download/")
data.dir.clinical.filter <- file.path(data.dir,"/step2_clinical_available_filtering/")
data.dir.probes.qc <- file.path(data.dir,"/step3_probesQC_filtering/")
data.dir.probes.normalization <- file.path(data.dir,"/step4_normalization/")
data.dir.pca <- file.path(data.dir,"/step5_pca_filtering/")
data.dir.neuron <- file.path(data.dir,"/step6_neuron_comp/")
data.dir.single.cpg.pval <- file.path(data.dir,"/step7_single_cpg_pval/")
data.dir.residuals <- file.path(data.dir,"/step10_residuals/")
data.dir.median <- file.path(data.dir,"/step11_median/")
for(p in grep("dir",ls(),value = T)) dir.create(get(p),recursive = TRUE,showWarnings = FALSE)Required R library GEOquery can be installed as following:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GEOquery")Data from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE80970
library(GEOquery)
dir.create(data.dir,recursive = TRUE,showWarnings = FALSE)
test <- getGEO(GEO = "GSE80970",destdir = data.dir.raw)
assayData <- exprs(test$GSE80970_series_matrix.txt.gz)
saveRDS(assayData, paste0(data.dir.raw, "GSE80970_assay.RDS"))
write.csv(assayData, paste0(data.dir.raw, "GSE80970_assay.csv"))
phenoData <- pData(test$GSE80970_series_matrix.txt.gz)
saveRDS(phenoData, paste0(data.dir.raw, "GSE80970_pheno.RDS"))
write.csv(phenoData, paste0(data.dir.raw, "GSE80970_pheno.csv"))Subset Files:
subset methylation data to match subjects in pheno data (n = 142)
Input: GSE80970_pheno.csv, GSE80970_assay.RDS
Output: pheno_df.RDS, beta_mat.RDS
### Read in phenotype data
phenoRaw_df <- readr::read_csv( paste0(data.dir.raw, "GSE80970_pheno.csv"),
col_types = readr::cols())
colnames(phenoRaw_df) <- gsub(" |[[:punct:]]",".",colnames(phenoRaw_df))
dim(phenoRaw_df)## [1] 286 46### Subset rows and columns
phenoPFC_df <- phenoRaw_df %>%
dplyr::filter(tissue.ch1 == "prefrontal cortex") %>%
dplyr::select(c("geo.accession",
"donor.id.ch1",
"sentrix.id.ch1",
"age..yr..ch1",
"gender.ch1",
"braak.stage.ch1"))
### Rename vars
colnames(phenoPFC_df) <- c(
"sample", "subject.id", "sentrix_id", "age.brain", "sex", "stage"
)
### Get slide from sentrix_id
# e.g. "7796806144_R01C01"(sentrix_id) -- "7796806144"(slide) and "R01C01"(array)
phenoPFC_df$slide <- gsub("_[[:alnum:]]*$","",phenoPFC_df$sentrix_id)
### Order final pheno_df
pheno_df <- phenoPFC_df[, c(
"sample", "subject.id", "sentrix_id", "slide", "age.brain", "sex", "stage"
)] ##### 2. Subset methylation data ###############################################
### Read in methylation data
beta_mat <- readRDS(paste0(data.dir.raw, "GSE80970_assay.RDS"))
### Subset methylation data based on pheno data
beta_mat <- beta_mat[, pheno_df$sample]
all(colnames(beta_mat) == pheno_df$sample)## [1] TRUEnb.probes <- nrow(beta_mat)
nb.samples <- ncol(beta_mat)
nb.samples.with.clinical <- ncol(beta_mat)##### 3. Output datasets #######################################################
## phenotype dataset
saveRDS(pheno_df, paste0(data.dir.clinical.filter, "pheno_df.RDS"))
## methylation beta values dataset
saveRDS(beta_mat, paste0(data.dir.clinical.filter, "beta142_mat.RDS"))Input: beta_mat.RDS Output: beta_CG_XY_SNPfiltered_mat.RDS
##### 1. keep on probes with start with "cg" ###################################
beta142_mat <- readRDS(paste0(data.dir.clinical.filter, "beta142_mat.RDS")) #dim: 485577 142
dim(beta142_mat)## [1] 485577 142beta142cg_mat <- beta142_mat[grep("cg",rownames(beta142_mat)),]
dim(beta142cg_mat)## [1] 482421 142nb.cg.probes <- nrow(beta142cg_mat)##### 2. drop probes that are on X/Y ###########################################
##### 3. drop probes where SNP with MAF >= 0.01 in the last 5 bp of the probe ##
library(DMRcate)beta142Filtered_mat <- rmSNPandCH(
object = beta142cg_mat,
dist = 5,
mafcut = 0.01,
and = TRUE,
rmcrosshyb = FALSE,
rmXY = TRUE
) #dim: 437713 142
dim(beta142Filtered_mat)## [1] 450793 142nb.probes.cg.dmrcate <- nrow(beta142Filtered_mat)##### 4. Output datasets #######################################################
saveRDS(
beta142Filtered_mat,
paste0(data.dir.probes.qc, "beta142_CG_XY_SNPfiltered_mat.RDS")
)Input:
Output:
beta_mat <- readRDS(paste0(data.dir.probes.qc, "beta142_CG_XY_SNPfiltered_mat.RDS"))
pheno_df <- readRDS(paste0(data.dir.clinical.filter, "pheno_df.RDS"))library(lumi)
betaQN <- lumiN(x.lumi = beta_mat, method = "quantile")##### 5. BMIQ ##################################################################
library(wateRmelon)
library(RPMM)
library(sesame)
library(sesameData)### Order annotation in the same order as beta matrix
annotType <- sesameDataGet("HM450.hg19.manifest")
annotType$designTypeNumeric <- ifelse(annotType$designType == "I",1,2)
### Density plot for type I and type II probes
library(sm)
betaQNCompleteCol1 <- betaQN[complete.cases(betaQN[,1]), ]
annotTypeCompleteCol1 <- annotType[row.names(betaQNCompleteCol1), ]
sm.density.compare(
betaQNCompleteCol1[,1],
annotTypeCompleteCol1$designTypeNumeric
)
type12 <- annotType$designTypeNumeric[match(rownames(betaQN),names(annotType))]### BMIQ
set.seed (946)
betaQN_BMIQ <- apply(
betaQN, 2,
function(x){
norm_ls <- BMIQ(x, design.v = type12, plots = FALSE)
return (norm_ls$nbeta)
}
)
saveRDS(betaQN_BMIQ, paste0(data.dir.probes.normalization, "MtSinai_QNBMIQ.RDS"))Description:
Input:
Output:
# plotPCA and OrderDataBySd functions
devtools::source_gist("https://gist.github.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e")## Sourcing https://gist.githubusercontent.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e/raw/a14424662da343c1301b7b2f03210d28d16ae05c/functions.R## SHA-1 hash of file is ef6f39dc4e5eddb5ca1c6e5af321e75ff06e9362##### 0.Import datasets ########################################################
beta_mat <- readRDS(paste0(data.dir.probes.normalization, "MtSinai_QNBMIQ.RDS")) #dim: 437713 142
pheno_df <- readRDS(paste0(data.dir.clinical.filter, "pheno_df.RDS")) #dim: 142 7
identical(colnames(beta_mat), pheno_df$sample)## [1] TRUE### transform to m values
mvalue_mat <- log2(beta_mat / (1 - beta_mat)) #dim: 437713 142
both_df <- subset(pheno_df, pheno_df$sample %in% colnames(beta_mat)) #dim: 142 7
# both_df$stage <- as.numeric(as.character(both_df$stage))
both_df$stage3 <- both_df$stage
both_df$stage3[both_df$stage <= 2] <- '0-2'
both_df$stage3[both_df$stage > 2 & both_df$stage < 5] <- '3-4'
both_df$stage3[both_df$stage >= 5] <- '5-6'
##### 1.Order matrix by most variable probes on top ############################
betaOrd_mat <- OrderDataBySd(beta_mat) #dim: 391482 142
mOrd_mat <- OrderDataBySd(mvalue_mat) #dim: 391482 142
betaOrd_matPlot <- betaOrd_mat[, both_df$sample] #dim: 391482 142
mOrd_matPlot <- mOrd_mat[, both_df$sample] #dim: 391482 142
identical(both_df$sample, colnames(betaOrd_matPlot))## [1] TRUEidentical(both_df$sample, colnames(mOrd_matPlot))## [1] TRUEexpSorted_mat = betaOrd_mat #dim: 391482 142
pca <- prcomp(t(expSorted_mat[1:50000, ]),
center = TRUE,
scale = TRUE)
d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2])
meanPC1 <- mean (d$PC1)
sdPC1 <- sd (d$PC1)
meanPC2 <- mean (d$PC2)
sdPC2 <- sd (d$PC2)
out3sdPC1_1 <- meanPC1 - 3 * sdPC1
out3sdPC1_2 <- meanPC1 + 3 * sdPC1
out3sdPC2_1 <- meanPC2 - 3 * sdPC2
out3sdPC2_2 <- meanPC2 + 3 * sdPC2
d$outlier_PC1[d$PC1 >= out3sdPC1_1 & d$PC1 <= out3sdPC1_2] <- 0
d$outlier_PC1[d$PC1 < out3sdPC1_1 | d$PC1 > out3sdPC1_2] <- 1
d$outlier_PC2[d$PC2 >= out3sdPC2_1 & d$PC2 <= out3sdPC2_2] <- 0
d$outlier_PC2[d$PC2 < out3sdPC2_1 | d$PC2 > out3sdPC2_2] <- 1
readr::write_csv(d, paste0(data.dir.pca, "MtSinai_PCs_usingBetas.csv"))##### 3.Filter samples by PCA, save files ######################################
noOutliers <- d[which(d$outlier_PC1 == 0 & d$outlier_PC2 == 0), ]
betaQN_BMIQ_PCfiltered <- beta_mat[, rownames(noOutliers)] #dim: 437713 141
saveRDS(betaQN_BMIQ_PCfiltered, paste0(data.dir.pca, "MtSinai_QNBMIQ_PCfiltered.RDS"))
pheno141_df <- pheno_df[pheno_df$sample %in% rownames(noOutliers),] #dim: 141 7
saveRDS(pheno141_df, paste0(data.dir.pca, "pheno141_df.RDS"))nb.samples.with.clinical.after.pca <- ncol(betaQN_BMIQ_PCfiltered)
dim(betaQN_BMIQ_PCfiltered)## [1] 434610 141dim(pheno141_df)## [1] 141 7pheno141_df %>%
DT::datatable(filter = 'top',
style = "bootstrap",
extensions = 'Buttons',
options = list(scrollX = TRUE,
dom = 'Bfrtip',
buttons = I('colvis'),
keys = TRUE,
pageLength = 10),
rownames = FALSE,
caption = "Samples metadata")df.samples <- data.frame("Number of samples" = c(nb.samples,
nb.samples.with.clinical,
nb.samples.with.clinical.after.pca),
"Description" = c("total number of samples",
"samples with clinical data",
"Samples after PCA"),
"Difference" = c("-",
nb.samples.with.clinical - nb.samples ,
nb.samples.with.clinical.after.pca - nb.samples.with.clinical)
)
df.samples # Create summary table
df.probes <- data.frame("Number of probes" = c(nb.probes,
nb.cg.probes,
nb.probes.cg.dmrcate),
"Description" = c("total number of probes in raw data",
"only probes that start with cg",
"DMRcate"),
"Difference" = c("-",
nb.cg.probes - nb.probes ,
nb.probes.cg.dmrcate - nb.cg.probes)
)
df.probessave(df.samples,df.probes,file = file.path(data.dir.table, "MtSiani_table.rda"))Data from https://www.tandfonline.com/doi/full/10.4161/epi.23924
objects <- load("../../CET/CETS_Image.RData")
objects## [1] "a" "a1" "b" "b1"
## [5] "brain" "c" "dat" "deconvolute"
## [9] "dilution" "estProportion" "f" "g2"
## [13] "getF" "getReference" "i" "idx"
## [17] "modelIdx" "p" "pdBrain" "prop"
## [21] "r2" "refProfile" "res" "rowttests"
## [25] "variables"idx <- list(
controlNeuron = pdBrain$celltype == "N" & pdBrain$diag == "Control" & pdBrain$ethnicity == "Caucasian",
controlGlia = pdBrain$celltype == "G" & pdBrain$diag == "Control" & pdBrain$ethnicity == "Caucasian"
)
refProfile <- getReference(brain, idx)
##### 2. Estimate proportions of neurons in PFC samples ########################
### Limit to 10,000 cpgs in the refProfile dataset
pfc <- readRDS(paste0(data.dir.pca, "MtSinai_QNBMIQ_PCfiltered.RDS")) #dim: 433656 59
selected <- rownames(pfc) %in% rownames(refProfile)
pfc.refcpgs <- pfc[selected, ] #dim: 9530 59
### Estimate proportion of neurons
prop <- data.frame(estProportion(pfc.refcpgs, profile = refProfile))
colnames(prop) <- "prop.neuron"
##### 3. Merge pfc.refcpgs with phenotype file #################################
pheno <- readRDS(paste0(data.dir.pca, "pheno141_df.RDS"))
pheno_final <- merge(
pheno,
prop,
by.x = "sample",
by.y = "row.names"
)
saveRDS(pheno_final, paste0(data.dir.neuron, "pheno141_withNeuronProp_df.RDS"))Input:
Output:
beta_mat <- readRDS(paste0(data.dir.pca, "MtSinai_QNBMIQ_PCfiltered.RDS")) #dim:433656 59
pheno_df <- readRDS(paste0(data.dir.neuron, "pheno141_withNeuronProp_df.RDS")) #dim:59 7### Compute M values
mval_mat <- log2(beta_mat / (1 - beta_mat))
pheno_df$Sample <- pheno_df$sample
identical(pheno_df$Sample, colnames(mval_mat))## [1] TRUEpheno_df$sex <- as.factor(pheno_df$sex)
pheno_df$slide <- as.factor(pheno_df$slide)
# If rosmap cohort, don't forget batch effect
str(pheno_df)## 'data.frame': 141 obs. of 9 variables:
## $ sample : chr "GSM2139163" "GSM2139164" "GSM2139165" "GSM2139166" ...
## $ subject.id : chr "MS66" "MS546" "MS709" "MS1549" ...
## $ sentrix_id : chr "7796806144_R01C01" "7796806144_R01C02" "7796806144_R02C01" "7796806144_R03C01" ...
## $ slide : Factor w/ 13 levels "7796806144","7796814008",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ age.brain : num 89 80 94 78 86 84 85 83 78 86 ...
## $ sex : Factor w/ 2 levels "female","male": 2 1 1 2 1 1 2 1 2 1 ...
## $ stage : num 5 2 3 4 3 4 5 1 2 6 ...
## $ prop.neuron: num 0.243 0.202 0.289 0.304 0.201 ...
## $ Sample : chr "GSM2139163" "GSM2139164" "GSM2139165" "GSM2139166" ...is(pheno_df$stage,"numeric")## [1] TRUEis(pheno_df$age.brain,"numeric")## [1] TRUEis(pheno_df$prop.neuron,"numeric")## [1] TRUEpredictors_char <- "stage"
covariates_char <- c("age.brain", "sex", "prop.neuron", "slide")
doParallel::registerDoParallel(cores = parallel::detectCores()/2)
devtools::source_gist("https://gist.github.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e")
results_ordered_df <- plyr::adply(mval_mat,1, function(row){
sumOneRegion_df <- data.frame(t(row))
result <- TestSingleRegion(
predictors_char = predictors_char,
covariates_char = covariates_char,
pheno_df = pheno_df,
sumOneRegion_df = sumOneRegion_df
)
result
}, .progress = "time",.parallel = T,.id = "cpg")
colnames(results_ordered_df)[1] <- "cpg"
identical(row.names(mval_mat), results_ordered_df$cpg %>% as.character())
results_ordered_df$fdr <- p.adjust(
results_ordered_df$pValue,
method = "fdr"
)
write.csv(
results_ordered_df,
paste0(data.dir.single.cpg.pval, "MtSinai_single_cpg_pVal_df.csv"),
row.names = FALSE
)results_ordered_df <- readr::read_csv(paste0(data.dir.single.cpg.pval, "MtSinai_single_cpg_pVal_df.csv"))## Parsed with column specification:
## cols(
## cpg = col_character(),
## Estimate = col_double(),
## StdErr = col_double(),
## pValue = col_double(),
## fdr = col_double()
## )results_ordered_dfInput:
Output:
##### 1. Import datasets #######################################################
beta_mat <- readRDS(grep("QNBMIQ",dir(data.dir.pca,full.names = T),ignore.case = T,value = T))
pheno_df <- readRDS(dir(data.dir.neuron,full.names = T))
### Compute M values
mvalue_mat <- log2(beta_mat / (1 - beta_mat))
### Reorder samples based on pheno_df
mvalue_mat <- mvalue_mat[, pheno_df$sample]
identical(colnames(mvalue_mat), pheno_df$sample)
### Take residuals
lmF <- function(mval){
fitE <- lm(
as.numeric(mval) ~ age.brain + sex + prop.neuron + as.character(slide), #add batch if rosmap
data = pheno_df,
na.action = na.exclude
)
residuals (fitE)
}
doParallel::registerDoParallel(cores = detectCores()/2)
resid <- plyr::adply(mvalue_mat,1,.fun = lmF,.progress = "time",.parallel = TRUE)
rownames(resid) <- resid[,1]
resid[,1] <- NULL
colnames(resid) <- colnames(mvalue_mat)
dim(resid)
# [1] 433656 22
dim(mvalue_mat)
# [1] 433656 22
### Save dataset
saveRDS(
resid,
paste0(data.dir.residuals,
"MtSinai_QNBMIQ_PCfiltered_mvalResiduals.RDS")
)
##### 4. Find co-methylated regions ############################################
### Import datasets
mvalue_residuals_mat <- readRDS(
paste0(data.dir.residuals,
"MtSinai_QNBMIQ_PCfiltered_mvalResiduals.RDS")
)
### Call in functions
library(coMethDMR)
library(BiocParallel)
probes.cluster.all <- coMethDMR::getPredefinedCluster(arrayType = "450k",
clusterType = "regions")
ncores <- 10
### Find co-methylated clusters
coMeth_ls <- CoMethAllRegions(
dnam = mvalue_residuals_mat,
betaToM = FALSE,
CpGs_ls = probes.cluster.all,
arrayType = "450k",
rDropThresh_num = 0.4,
minPairwiseCorr = NULL,
method = "spearman",
returnAllCpGs = TRUE,
output = "all",
nCores_int = ncores,
progressbar = FALSE
)
saveRDS(
coMeth_ls,
paste0(data.dir.residuals,
"MtSinai_residuals_cometh_input_ls.RDS")
)Input:
Output:
### Import datasets
beta_mat <- readRDS(dir(data.dir.pca, pattern = "QNBMIQ", full.names = TRUE))
pheno_df <- readRDS(dir(data.dir.neuron,full.names = T))
mval_mat <- log2(beta_mat / (1 - beta_mat)) %>% as.matrix()
coMeth_ls <- readRDS(
paste0(data.dir.residuals,
"MtSinai_residuals_cometh_input_ls.RDS")
)
### Create info dataset
input_cometh <- data.frame(
inputRegion = coMeth_ls$inputRegion_chr,
nCoMethRegion = coMeth_ls$nCoMethRegions_num,
coMethRegion = names(coMeth_ls$coMeth_ls),
nCpGs = unlist(lapply(coMeth_ls$coMeth_ls, length), use.names = FALSE),
stringsAsFactors = FALSE
)
input_cometh_nodup <- input_cometh[
!duplicated(input_cometh$coMethRegion),
]
colnames(input_cometh_nodup) <- c(
paste0(cohort, "_inputRegion"),
paste0(cohort, "_nCoMethRegion"),
paste0(cohort, "_coMethRegion"),
paste0(cohort, "_nCpGs")
)
saveRDS(
input_cometh_nodup,
paste0(data.dir.median, cohort, "_info_df.rds")
)
### Take median of probes in each cluster for each sample
filename <- paste0(paste0(data.dir.median,cohort, "_mediansMval_df.rds"))
library(robustbase)
mval_mat <- mval_mat[rownames(mval_mat) %in% unlist(coMeth_ls$coMeth_ls),]
if(!file.exists(filename)){
medianMval.df <- plyr::ldply(
coMeth_ls$coMeth_ls[!duplicated(names(coMeth_ls$coMeth_ls))],
function(probes){
colMedians(mval_mat[as.character(probes),], na.rm = TRUE)
},
.progress = "time"
)
medianMval.df$.id <- NULL
colnames(medianMval.df) <- colnames(mval_mat)
saveRDS(medianMval.df, file = filename)
} else {
medianMval.df <- readRDS(filename)
}### Import datasets
info_df <- readRDS(dir(data.dir.median,pattern = "info", full.names = TRUE))
mediansMval_df <- readRDS(dir(data.dir.median,pattern = "mediansMval", full.names = TRUE))
pheno_df <- readRDS(dir(data.dir.neuron, full.names = TRUE))
### Check variables before fitting model
pheno_df$Sample <- pheno_df$sample
identical(pheno_df$Sample, colnames(mediansMval_df))## [1] TRUEpheno_df$sex <- as.factor(pheno_df$sex)
pheno_df$slide <- as.factor(pheno_df$slide)
# If rosmap cohort, don't forget batch effect
str(pheno_df)## 'data.frame': 141 obs. of 9 variables:
## $ sample : chr "GSM2139163" "GSM2139164" "GSM2139165" "GSM2139166" ...
## $ subject.id : chr "MS66" "MS546" "MS709" "MS1549" ...
## $ sentrix_id : chr "7796806144_R01C01" "7796806144_R01C02" "7796806144_R02C01" "7796806144_R03C01" ...
## $ slide : Factor w/ 13 levels "7796806144","7796814008",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ age.brain : num 89 80 94 78 86 84 85 83 78 86 ...
## $ sex : Factor w/ 2 levels "female","male": 2 1 1 2 1 1 2 1 2 1 ...
## $ stage : num 5 2 3 4 3 4 5 1 2 6 ...
## $ prop.neuron: num 0.243 0.202 0.289 0.304 0.201 ...
## $ Sample : chr "GSM2139163" "GSM2139164" "GSM2139165" "GSM2139166" ...predictors_char <- "stage"
covariates_char <- c("age.brain", "sex", "prop.neuron", "slide")
devtools::source_gist("https://gist.github.com/tiagochst/d3a7b1639acf603916c315d23b1efb3e")
res_df <- TestAllRegions_noInfo(
predictors_char = predictors_char,
covariates_char = covariates_char,
pheno_df = pheno_df,
summarizedRegions_df = mediansMval_df
)
colnames(res_df) <- c(
paste0(cohort, "_estimate"),
paste0(cohort, "_se"),
paste0(cohort, "_pVal"),
paste0(cohort, "_fdr")
)
res_withInfo_df <- cbind(info_df, res_df)
saveRDS(
res_withInfo_df,
paste0(data.dir.median,cohort, "_linear_df.rds")
)file <- dir(data.dir.median,pattern = paste0(".*linear_df"),
recursive = T,
full.names = TRUE,
ignore.case = T)
file## [1] "DATASETS//MtSinai////step11_median//MtSinai_linear_df.rds"res_withInfo_df <- readRDS(file)
dim(res_withInfo_df)## [1] 38769 8res_withInfo_dfdir(path = data.dir,recursive = T,pattern = ".rda|.csv|.RDS")## [1] "step1_download/GSE80970_assay.csv"
## [2] "step1_download/GSE80970_assay.RDS"
## [3] "step1_download/GSE80970_pheno.csv"
## [4] "step1_download/GSE80970_pheno.RDS"
## [5] "step10_residuals/MtSinai_QNBMIQ_PCfiltered_mvalResiduals.RDS"
## [6] "step10_residuals/MtSinai_residuals_cometh_input_ls.RDS"
## [7] "step2_clinical_available_filtering/beta142_mat.RDS"
## [8] "step2_clinical_available_filtering/pheno_df.RDS"
## [9] "step3_probesQC_filtering/beta142_CG_XY_SNPfiltered_mat.RDS"
## [10] "step4_normalization/MtSinai_QNBMIQ.RDS"
## [11] "step4_normalization/qtestPerm.rda"
## [12] "step5_pca_filtering/MtSinai_PCs_usingBetas.csv"
## [13] "step5_pca_filtering/MtSinai_QNBMIQ_PCfiltered.RDS"
## [14] "step5_pca_filtering/pheno141_df.RDS"
## [15] "step6_neuron_comp/pheno141_withNeuronProp_df.RDS"
## [16] "step7_single_cpg_pval/MtSinai_single_cpg_pVal_df.csv"devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R Under development (unstable) (2020-02-25 r77857)
## os macOS Catalina 10.15.3
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2020-04-26
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib
## acepack 1.4.1 2016-10-29 [1]
## affy 1.65.1 2019-11-06 [1]
## affyio 1.57.0 2019-11-06 [1]
## annotate 1.65.1 2020-01-27 [1]
## AnnotationDbi * 1.49.1 2020-01-25 [1]
## AnnotationFilter 1.11.0 2019-11-06 [1]
## AnnotationHub * 2.19.11 2020-04-16 [1]
## askpass 1.1 2019-01-13 [1]
## assertthat 0.2.1 2019-03-21 [1]
## backports 1.1.6 2020-04-05 [1]
## base64 2.0 2016-05-10 [1]
## base64enc 0.1-3 2015-07-28 [1]
## beanplot 1.2 2014-09-19 [1]
## Biobase * 2.47.3 2020-03-16 [1]
## BiocFileCache * 1.11.6 2020-04-16 [1]
## BiocGenerics * 0.33.3 2020-03-23 [1]
## BiocManager 1.30.10 2019-11-16 [1]
## BiocParallel 1.21.2 2019-12-21 [1]
## BiocVersion 3.11.1 2019-11-13 [1]
## biomaRt 2.43.5 2020-04-02 [1]
## Biostrings * 2.55.7 2020-03-24 [1]
## biovizBase 1.35.1 2019-12-03 [1]
## bit 1.1-15.2 2020-02-10 [1]
## bit64 0.9-7 2017-05-08 [1]
## bitops 1.0-6 2013-08-17 [1]
## blob 1.2.1 2020-01-20 [1]
## BSgenome 1.55.4 2020-03-19 [1]
## bsseq 1.23.2 2020-04-06 [1]
## bumphunter * 1.29.0 2019-11-07 [1]
## callr 3.4.3 2020-03-28 [1]
## cellranger 1.1.0 2016-07-27 [1]
## checkmate 2.0.0 2020-02-06 [1]
## cli 2.0.2 2020-02-28 [1]
## cluster * 2.1.0 2019-06-19 [1]
## codetools 0.2-16 2018-12-24 [1]
## colorspace 1.4-1 2019-03-18 [1]
## crayon 1.3.4 2017-09-16 [1]
## crosstalk 1.1.0.1 2020-03-13 [1]
## curl 4.3 2019-12-02 [1]
## data.table 1.12.9 2020-02-26 [1]
## DBI 1.1.0 2019-12-15 [1]
## dbplyr * 1.4.3 2020-04-19 [1]
## DelayedArray * 0.13.12 2020-04-10 [1]
## DelayedMatrixStats 1.9.1 2020-03-30 [1]
## desc 1.2.0 2018-05-01 [1]
## devtools 2.3.0 2020-04-10 [1]
## dichromat 2.0-0 2013-01-24 [1]
## digest 0.6.25 2020-02-23 [1]
## DMRcate * 2.1.9 2020-03-28 [1]
## DMRcatedata * 2.5.0 2019-10-31 [1]
## DNAcopy 1.61.0 2019-11-06 [1]
## doRNG 1.8.2 2020-01-27 [1]
## dplyr * 0.8.99.9002 2020-04-02 [1]
## DSS 2.35.1 2020-04-14 [1]
## DT 0.13 2020-03-23 [1]
## edgeR 3.29.1 2020-02-26 [1]
## ellipsis 0.3.0 2019-09-20 [1]
## ensembldb 2.11.4 2020-04-17 [1]
## evaluate 0.14 2019-05-28 [1]
## ExperimentHub * 1.13.7 2020-04-16 [1]
## fansi 0.4.1 2020-01-08 [1]
## fastmap 1.0.1 2019-10-08 [1]
## FDb.InfiniumMethylation.hg19 * 2.2.0 2020-04-09 [1]
## foreach * 1.5.0 2020-03-30 [1]
## foreign 0.8-78 2020-04-13 [1]
## Formula 1.2-3 2018-05-03 [1]
## fs 1.4.1 2020-04-04 [1]
## genefilter 1.69.0 2019-11-06 [1]
## generics 0.0.2 2018-11-29 [1]
## GenomeInfoDb * 1.23.17 2020-04-13 [1]
## GenomeInfoDbData 1.2.3 2020-04-20 [1]
## GenomicAlignments 1.23.2 2020-03-24 [1]
## GenomicFeatures * 1.39.7 2020-03-19 [1]
## GenomicRanges * 1.39.3 2020-04-08 [1]
## GEOquery 2.55.1 2019-11-18 [1]
## ggplot2 * 3.3.0 2020-03-05 [1]
## glue 1.4.0 2020-04-03 [1]
## gridExtra 2.3 2017-09-09 [1]
## gtable 0.3.0 2019-03-25 [1]
## gtools 3.8.2 2020-03-31 [1]
## Gviz 1.31.12 2020-03-05 [1]
## HDF5Array 1.15.18 2020-04-10 [1]
## Hmisc 4.4-0 2020-03-23 [1]
## hms 0.5.3 2020-01-08 [1]
## htmlTable 1.13.3 2019-12-04 [1]
## htmltools 0.4.0 2019-10-04 [1]
## htmlwidgets 1.5.1 2019-10-08 [1]
## httpuv 1.5.2 2019-09-11 [1]
## httr 1.4.1 2019-08-05 [1]
## IlluminaHumanMethylation450kanno.ilmn12.hg19 * 0.6.0 2020-03-24 [1]
## IlluminaHumanMethylationEPICanno.ilm10b4.hg19 0.6.0 2020-04-09 [1]
## illuminaio * 0.29.0 2019-11-06 [1]
## interactiveDisplayBase 1.25.0 2019-11-06 [1]
## IRanges * 2.21.8 2020-03-25 [1]
## iterators * 1.0.12 2019-07-26 [1]
## jpeg 0.1-8.1 2019-10-24 [1]
## jsonlite 1.6.1 2020-02-02 [1]
## KernSmooth 2.23-16 2019-10-15 [1]
## knitr 1.28 2020-02-06 [1]
## later 1.0.0 2019-10-04 [1]
## lattice 0.20-41 2020-04-02 [1]
## latticeExtra 0.6-29 2019-12-19 [1]
## lazyeval 0.2.2 2019-03-15 [1]
## lifecycle 0.2.0 2020-03-06 [1]
## limma * 3.43.8 2020-04-14 [1]
## locfit * 1.5-9.4 2020-03-25 [1]
## lumi * 2.39.3 2020-03-27 [1]
## magrittr 1.5 2014-11-22 [1]
## MASS 7.3-51.5 2019-12-20 [1]
## Matrix 1.2-18 2019-11-27 [1]
## matrixStats * 0.56.0 2020-03-13 [1]
## mclust 5.4.6 2020-04-11 [1]
## memoise 1.1.0 2017-04-21 [1]
## methylumi * 2.33.0 2019-11-06 [1]
## mgcv 1.8-31 2019-11-09 [1]
## mime 0.9 2020-02-04 [1]
## minfi * 1.33.1 2020-03-05 [1]
## missMethyl 1.21.4 2020-01-28 [1]
## multtest 2.43.1 2020-03-12 [1]
## munsell 0.5.0 2018-06-12 [1]
## nleqslv 3.3.2 2018-05-17 [1]
## nlme 3.1-147 2020-04-13 [1]
## nnet 7.3-13 2020-02-25 [1]
## nor1mix 1.3-0 2019-06-13 [1]
## openssl 1.4.1 2019-07-18 [1]
## org.Hs.eg.db * 3.10.0 2020-03-30 [1]
## permute 0.9-5 2019-03-12 [1]
## pillar 1.4.3 2019-12-20 [1]
## pkgbuild 1.0.6 2019-10-09 [1]
## pkgconfig 2.0.3 2019-09-22 [1]
## pkgload 1.0.2 2018-10-29 [1]
## plyr 1.8.6 2020-03-03 [1]
## png 0.1-7 2013-12-03 [1]
## preprocessCore 1.49.2 2020-02-01 [1]
## prettyunits 1.1.1 2020-01-24 [1]
## processx 3.4.2 2020-02-09 [1]
## progress 1.2.2 2019-05-16 [1]
## promises 1.1.0 2019-10-04 [1]
## ProtGenerics 1.19.3 2019-12-25 [1]
## ps 1.3.2 2020-02-13 [1]
## purrr 0.3.4 2020-04-17 [1]
## quadprog 1.5-8 2019-11-20 [1]
## R.methodsS3 1.8.0 2020-02-14 [1]
## R.oo 1.23.0 2019-11-03 [1]
## R.utils 2.9.2 2019-12-08 [1]
## R6 2.4.1 2019-11-12 [1]
## randomForest 4.6-14 2018-03-25 [1]
## rappdirs 0.3.1 2016-03-28 [1]
## RColorBrewer 1.1-2 2014-12-07 [1]
## Rcpp 1.0.4.6 2020-04-09 [1]
## RCurl 1.98-1.2 2020-04-18 [1]
## readr 1.3.1 2018-12-21 [1]
## readxl 1.3.1 2019-03-13 [1]
## remotes 2.1.1 2020-02-15 [1]
## reshape 0.8.8 2018-10-23 [1]
## reshape2 * 1.4.4 2020-04-09 [1]
## rhdf5 2.31.10 2020-04-02 [1]
## Rhdf5lib 1.9.3 2020-04-15 [1]
## rlang 0.4.5.9000 2020-03-20 [1]
## rmarkdown 2.1 2020-01-20 [1]
## rngtools 1.5 2020-01-23 [1]
## ROC * 1.63.0 2019-11-06 [1]
## rpart 4.1-15 2019-04-12 [1]
## RPMM * 1.25 2017-02-28 [1]
## rprojroot 1.3-2 2018-01-03 [1]
## Rsamtools 2.3.7 2020-03-18 [1]
## RSQLite 2.2.0 2020-01-07 [1]
## rstudioapi 0.11 2020-02-07 [1]
## rtracklayer 1.47.0 2019-11-06 [1]
## S4Vectors * 0.25.15 2020-04-04 [1]
## scales * 1.1.0 2019-11-18 [1]
## scrime 1.3.5 2018-12-01 [1]
## sesame * 1.5.3 2020-03-03 [1]
## sesameData * 1.5.0 2019-10-31 [1]
## sessioninfo 1.1.1 2018-11-05 [1]
## shiny 1.4.0.2 2020-03-13 [1]
## siggenes 1.61.0 2019-11-06 [1]
## statmod 1.4.34 2020-02-17 [1]
## stringi 1.4.6 2020-02-17 [1]
## stringr 1.4.0 2019-02-10 [1]
## SummarizedExperiment * 1.17.5 2020-03-27 [1]
## survival 3.1-12 2020-04-10 [1]
## testthat 2.3.2 2020-03-02 [1]
## tibble 3.0.1 2020-04-20 [1]
## tidyr 1.0.2 2020-01-24 [1]
## tidyselect 1.0.0 2020-01-27 [1]
## TxDb.Hsapiens.UCSC.hg19.knownGene * 3.2.2 2020-04-02 [1]
## usethis 1.6.0 2020-04-09 [1]
## VariantAnnotation 1.33.4 2020-04-09 [1]
## vctrs 0.2.99.9010 2020-04-02 [1]
## wateRmelon * 1.31.0 2019-11-06 [1]
## wheatmap 0.1.0 2018-03-15 [1]
## withr 2.1.2 2018-03-15 [1]
## xfun 0.13 2020-04-13 [1]
## XML 3.99-0.3 2020-01-20 [1]
## xml2 1.3.1 2020-04-09 [1]
## xtable 1.8-4 2019-04-21 [1]
## XVector * 0.27.2 2020-03-24 [1]
## yaml 2.2.1 2020-02-01 [1]
## zlibbioc 1.33.1 2020-01-24 [1]
## source
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## local
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## Github (tidyverse/dplyr@affb977)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## Github (Bioconductor/GenomicRanges@70e6e69)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## Github (r-lib/rlang@a90b04b)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
## Github (r-lib/vctrs@fd24927)
## Bioconductor
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## CRAN (R 4.0.0)
## Bioconductor
## CRAN (R 4.0.0)
## Bioconductor
##
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library