vignettes/vignette_analysis.R
In cozygene/TCA: Tensor Composition Analysis
library("TCA")
library("ggplot2")
library("ggpubr")
library("pracma")
library("matrixStats")
prep_data <- function(data_path){
file_name1 <- paste(data_path,.Platform$file.sep,"hannum.chr22.RData",sep="")
file_name2 <- paste(data_path,.Platform$file.sep,"liu.cd4.chr22.RData",sep="")
file_name3 <- paste(data_path,.Platform$file.sep,"paquette.chr22.RData",sep="")
file_names <- c(file_name1, file_name2, file_name3)
if(any(!file.exists(file_name1, file_name2, file_name3))){
library("GEOquery")
library("data.table")
library("EpiDISH")
}
if (!file.exists(file_name1)){
# Download the Hannum et al. data
gse <- GEOquery::getGEO("GSE40279", destdir = data_path, GSEMatrix = TRUE)
# Extract methylation data
X.hannum <- Biobase::exprs(gse[[1]])
# covariates; note that there's also 'ethnicity' covariate in the data, however, it is perfectly captured by the plateinformation
plate.hannum <- as.numeric(as.factor((Biobase::pData(gse[[1]])[["plate:ch1"]])))
cov.hannum <- cbind(as.numeric((Biobase::pData(gse[[1]])[["age (y):ch1"]])), as.numeric(as.factor((Biobase::pData(gse[[1]])[["gender:ch1"]]))), indicator_vars(plate.hannum))
rownames(cov.hannum) <- colnames(X.hannum)
colnames(cov.hannum) <- c("age", "gender", "plate1", "plate2", "plate3", "plate4", "plate5", "plate6", "plate7", "plate8")
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.hannum <- low_var_pcs(X.hannum, rank = 10, p = 1000)
cov.hannum <- cbind(cov.hannum, low_var_pcs.hannum)
# estimate cell-type proportions usign a reference-based approach
ref <- as.matrix(EpiDISH::centDHSbloodDMC.m[,c("Neutro","Eosino","CD4T","CD8T","Mono","B","NK")])
W.hannum <- EpiDISH::epidish(X.hannum, ref)$estF
# merge Neutro and Eosino
W.hannum.gran <- W.hannum[,"Neutro",drop=F] + W.hannum[,"Eosino",drop=F]
colnames(W.hannum.gran) <- "Gran"
W.hannum <- cbind(W.hannum.gran,W.hannum[,c("CD4T","CD8T","Mono","B","NK")])
# remove polymorphic sites, cross-reactive sites, and non-autosomal sites according to Chen et al.; in addition, remove low variance sites, as those are unlikely to demonstrate biological signal.
X.hannum.processed <- filter_data(X.hannum)
# keep only sites on chromosome 22 (so that tutorial can run quickly)
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- map[rownames(X.hannum.processed),1]
chr22sites <- which(chrs == 22)
X.hannum.processed.22 <- X.hannum.processed[chr22sites,]
hannum <- list(X = X.hannum.processed.22, cov = cov.hannum, W = W.hannum)
save(hannum, file = file_name1, compress = "bzip2")
rm(hannum, gse, X.hannum, X.hannum.processed)
file.remove(paste(data_path,.Platform$file.sep,"GPL13534.soft",sep=""))
file.remove(paste(data_path,.Platform$file.sep,"GSE40279_series_matrix.txt.gz",sep=""))
}
if (!file.exists(file_name2)){
# Download the Liu et al. CD4 data
# covariates
gse <- GEOquery::getGEO("GSE56581", destdir = data_path, GSEMatrix = TRUE)
liu.cd4.age <- as.matrix(as.numeric(Biobase::pData(gse[[1]])[["age (yrs):ch1"]]))
colnames(liu.cd4.age) <- "age"
# methylation data
liu.cd4.methfile <- paste(data_path,.Platform$file.sep,"GSE56581_methylome_normalized.txt",sep="")
download.file("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE56nnn/GSE56581/suppl/GSE56581_methylome_normalized.txt.gz", paste(liu.cd4.methfile,".gz",sep=""))
gunzip(paste(liu.cd4.methfile,".gz",sep=""))
X.liu.cd4 <- fread(liu.cd4.methfile)
X.liu.cd4.rownames <- X.liu.cd4$ID_REF
X.liu.cd4 <- as.matrix(X.liu.cd4[,2:ncol(X.liu.cd4)])
rownames(X.liu.cd4) <- X.liu.cd4.rownames
X.liu.cd4 <- X.liu.cd4[,seq(1,ncol(X.liu.cd4),2)]
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.liu.cd4 <- low_var_pcs(X.liu.cd4, rank = 10)
cov.liu.cd4 <- cbind(liu.cd4.age, low_var_pcs.liu.cd4)
# keep only sites on chromosome 22 (so that tutorial can run quickly)
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- map[rownames(X.liu.cd4),1]
chr22sites <- which(chrs == 22)
X.liu.cd4.22 <- X.liu.cd4[chr22sites,]
liu.cd4 <- list(X = X.liu.cd4.22, cov = cov.liu.cd4)
save(liu.cd4, file = file_name2)
rm(liu.cd4, gse, X.liu.cd4.22, X.liu.cd4)
file.remove(paste(data_path,.Platform$file.sep,"GSE56581_methylome_normalized.txt",sep=""))
file.remove(paste(data_path,.Platform$file.sep,"GSE56581_series_matrix.txt.gz",sep=""))
}
if (!file.exists(file_name3)){
# Paquette et al. (Epigenetics 2016)
gse <- GEOquery::getGEO("GSE75248", destdir = data_path, GSEMatrix = TRUE)
# methylation data
X.paquette <- Biobase::exprs(gse[[1]])
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.paquette <- low_var_pcs(X.paquette, rank = 10, p = 1000)
# remove polymorphic sites, cross-reactive sites, and non-autosomal sites according to Chen et al.; in addition, remove low variance sites, as those are unlikely to demonstrate biological signal.
X.paquette.processed <- filter_data(X.paquette)
# keep only sites on chromosome 22 (so that tutorial can run quickly); Paquette et al. reported some signal in chromosome 16
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- as.character(map[rownames(X.paquette.processed),1])
chr22sites <- which(chrs == "22")
X.paquette.processed.22 <- X.paquette.processed[chr22sites,]
# Extract covaraites
gestational_age <-as.numeric(unlist(lapply(Biobase::pData(gse[[1]])[["gestational age:ch1"]], function(x) strsplit(x,";")[[1]][1])))
arousal <- as.numeric(Biobase::pData(gse[[1]])[["arousal:ch1"]])
batch <- as.numeric(as.factor((Biobase::pData(gse[[1]])[["batch:ch1"]])))
batch[is.na(batch)] <- 3
batch <- indicator_vars(batch)
# since gender information is not available on the GEO record, get gender information that was inferred from the IDAR files based on X chromosome methylation
gender <- read.table("https://raw.githubusercontent.com/cozygene/TCA/master/vignettes/gender.paquette.txt", row.names = 1)
cov.paquette <- data.frame(gender,gestational_age,arousal,batch)
colnames(cov.paquette) <- c("gender","gestational_age","arousal","batch1","batch2")
rownames(cov.paquette) <- colnames(X.paquette.processed.22)
# remove na values
keep <- setdiff(1:nrow(cov.paquette) ,which(rowSums(is.na(cov.paquette)) > 0))
cov.paquette <- cbind(cov.paquette[keep,], low_var_pcs.paquette[keep,])
X.paquette.processed.22 <- X.paquette.processed.22[,keep]
# Load cell-type proportions usign a reference-based approach; these were estimated using the raw IDAT files of the data (available on GEO) with the package ENmix.
W.paquette <- read.table("https://raw.githubusercontent.com/cozygene/TCA/master/vignettes/W.paquette.txt", row.names = 1)
W.paquette <- W.paquette[,2:ncol(W.paquette)]
W.paquette[W.paquette < 0] <- 0
# remove lowly abundant cell types
W.paquette <- W.paquette[keep, colMeans(W.paquette)>=0.01]
rownames(W.paquette) <- rownames(cov.paquette)
# normalize cell type proportions (output from the ENmix package didn't require the proportions of an individual to sum up to 1)
W.paquette <- W.paquette/t(repmat(rowSums(W.paquette),ncol(W.paquette),1))
# get the set of cord blood reference CpGs that were used for estimating cell-type proportions; we will use these CpGs for re-estimating W
library("FlowSorted.CordBlood.450k")
ref.cord <- get("FlowSorted.CordBlood.450k")
ref.cord <- preprocessRaw(ref.cord)
ref.cord.cpgs <- rownames(FlowSorted.CordBlood.450k.ModelPars)
ref.cord.cpgs.intersect <- intersect(ref.cord.cpgs, rownames(X.paquette))
paquette <- list(X = X.paquette.processed.22, cov = cov.paquette, W = W.paquette, X.ref_cpgs = X.paquette[ref.cord.cpgs.intersect,keep])
save(paquette, file = file_name3)
rm(paquette, gse, X.paquette, X.paquette.processed, X.paquette.processed.22)
file.remove(paste(data_path,.Platform$file.sep,"GSE75248_series_matrix.txt.gz",sep=""))
## note - the files W.paquette.txt and gender.paquette.txt were generated using the original IDAT files of the Paquette data.
# Gender information was inferred from methylation levels frmo the X chromosome as it is not available on Tthe GEO record.
# Code attched below; in order to run it, download the IDAT files from GEO accession number GSE75248 and set 'path' in 'readidat' below to the location of the files.
# require("ENmix")
# rgdat <- readidat(path = ".", manifestfile=NULL, recursive = FALSE, verbose = TRUE)
# meth <- getmeth(rgdat); rm(rgdat)
# W <- estimateCellProp(meth, refdata="FlowSorted.CordBlood.450k", cellTypes=NULL, nonnegative = TRUE, nProbes=50, normalize=TRUE, refplot=FALSE)
# write.table(W, file = "W.paquette.txt", quote = FALSE, sep = " ")
# beta <- getB(meth)
# map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
# chrs <- as.character(map[rownames(beta),1])
# sites.X <- which(chrs == "X")
# sex_cpgs.pca.paquette <- prcomp(t(beta[sites.X,]), center=TRUE, scale=TRUE, rank = 2)
# write.table(sex_cpgs.pca.paquette$x[,1:2], file = "sex_cpgs.pca.paquette.txt", quote = FALSE, sep = " ")
# plot(sex_cpgs.pca.paquette$x[,1],sex_cpgs.pca.paquette$x[,2]) # shows a separation between males and females; note that the males/females ratio here perfectly match the one in the Paquette et al. paper.
# df <-data.frame(as.numeric(x[,1] > 0))
# colnames(df) <-"gender"; rownames(df) <- rownames(x)
# write.table(df,file = "gender.paquette.txt", quote = FALSE, sep = " ")
}
return(file_names)
}
indicator_vars <- function(x){
x.indicators <- matrix(0,length(x),length(unique(x))-1)
counter <- 1
u <- unique(x)
for (i in setdiff(u,u[length(u)])){
x.indicators[,counter] <- as.numeric(x == i)
counter <- counter + 1
}
return(x.indicators)
}
low_var_pcs <- function(X, rank = 10, p = 1000){
site.variances <- matrixStats::rowVars(X)
names(site.variances) <- rownames(X)
low.var.sites <- names(head(sort(site.variances), p))
low.var.pca <- prcomp(t(X[low.var.sites,]), center=TRUE, scale=TRUE, rank = rank)
return(low.var.pca$x)
}
filter_data <- function(X, var_th = 0.0001){
nonspecific_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/nonspecific_probes.txt")[,1]
XY_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/HumanMethylationSites_X_Y.txt")[,1]
polymorphic_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/polymorphic_cpgs.txt")[,1]
# remove sites with very low variance that are unlikely to exhibit biological signal
site.variances <- matrixStats::rowVars(X)
names(site.variances) <- rownames(X)
low_var_sites <- names(which(site.variances < var_th))
exclude <- union(nonspecific_probes,union(XY_probes,union(low_var_sites,polymorphic_probes)))
return(X[setdiff(rownames(X),exclude),])
}
plot_qq <- function(pvals, labels, ggarrange.nrow = 1, ggarrange.ncol = 1, alpha = 0.05, experiment_wide_line = TRUE){
significance_th <- list(alpha/length(pvals[[1]]))
if(length(pvals)-1) significance_th[[2]] <- alpha/(length(pvals)*length(pvals[[1]]))
qqplots <- lapply(1:length(pvals), function(p){
df <- data.frame(pvals.obs = -log10(sort(pvals[[p]])), pvals.exp = -log10(sort((1:length(pvals[[1]]))/length(pvals[[1]]))));
qqplot <- ggplot(df, aes(x = pvals.exp, y = pvals.obs)) +
stat_binhex(geom = "point", bins=1000, size=1) +
geom_abline() +
ggtitle(labels[p]) +
xlab(expression(Expected~-log[10](P))) + ylab(expression(Observed~-log[10](P))) +
theme_bw() +
guides(fill="none") +
geom_hline(yintercept=-log10(significance_th[[1]]), linetype="dashed", color = "red", size=1)
if (length(significance_th)-1 & experiment_wide_line) qqplot <- qqplot + geom_hline(yintercept=-log10(significance_th[[2]]), linetype="dashed", color = "red", size=0.5)
return(qqplot)
})
ggarrange(plotlist = qqplots, ncol = ggarrange.ncol, nrow = ggarrange.nrow)
}
plot_scatter <- function(dfs, ggarrange.ncol, ggarrange.nrow, xlab, ylab, titles){
plots <- vector("list", length = length(dfs))
for (i in 1:length(dfs)){
df <- data.frame(y = dfs[[i]]$y, x = dfs[[i]]$x)
plots[[i]] <- ggplot(df, aes(x = x, y = y)) +
geom_point(alpha = 0.6) +
geom_smooth(method=lm) +
stat_cor(method = "pearson", colour = "blue") +
xlab(xlab) +
ylab(ylab) +
ggtitle(titles[i])
}
ggarrange(plotlist = plots, ncol = ggarrange.ncol, nrow = ggarrange.nrow)
}
if(FALSE){
# Set a path for storing data
data_path <- "./"
# Load the data
filenames <- prep_data(data_path)
for (filename in filenames) load(filename)
## Experiment #1: detecting CD4 differential methylation with age; working under the assumption X|Y (i.e. age affects methylation levels)
# Apply the TCA model to the Hannum whole-blood data
tca.mdl.hannum <- tca(X = hannum$X,
W = hannum$W,
C1 = hannum$cov[,c("gender","age")],
C2 = hannum$cov[,3:ncol(hannum$cov)])
# Extract p-values of a joint test
tca.mdl.hannum.pvals.joint <- tca.mdl.hannum$gammas_hat_pvals.joint[,"age"]
# Extract p-values for each cell type, under a marginal conditional test
tca.mdl.hannum.pvals.marg_cond <- tca.mdl.hannum$gammas_hat_pvals[,paste(colnames(hannum$W),".age",sep="")]
# qq-plots - for the p-values of the joint test, and for the p-values in CD4, under a marginal conditional test
plot_qq(list(tca.mdl.hannum.pvals.joint, tca.mdl.hannum.pvals.marg_cond[,"CD4T.age"]),
labels = c("Joint test with age", "CD4 marginal conditional test with age"),
ggarrange.nrow = 1,
ggarrange.ncol = 2,
experiment_wide_line = FALSE)
# Run ReFACTor to capture more of the variation of cell-type composition
refactor.mdl.hannum <- refactor(X = hannum$X,
k = 6,
C = hannum$cov[,3:ncol(hannum$cov)])
# Rerun TCA, this time include the ReFACTor components as additional covariates
tca.mdl.hannum.2 <- tca(X = hannum$X,
W = hannum$W,
C1 = hannum$cov[,c("gender","age")],
C2 = cbind(hannum$cov[,3:ncol(hannum$cov)],refactor.mdl.hannum$scores))
# Extract the updated p-values of a joint test
tca.mdl.hannum.2.pvals.joint <- tca.mdl.hannum.2$gammas_hat_pvals.joint[,"age"]
# Extract the updated marginal conditional p-values
tca.mdl.hannum.2.pvals.marg_cond <- tca.mdl.hannum.2$gammas_hat_pvals[,paste(colnames(hannum$W),".age",sep="")]
# qq-plots - for the p-values of the joint test, and for the p-values in CD4, under a marginal conditional test
plot_qq(list(tca.mdl.hannum.2.pvals.joint, tca.mdl.hannum.2.pvals.marg_cond[,"CD4T.age"]),
labels = c("Joint test with age", "CD4 marginal conditional test with age"),
ggarrange.nrow = 1,
ggarrange.ncol = 2,
experiment_wide_line = FALSE)
# Extract the hits based on the joint test
hits.joint <- names(which(tca.mdl.hannum.2.pvals.joint < 0.05/nrow(hannum$X)))
# Extract the hits from hits.joint where CD4 cells demonstrate the lowest p-value across all cell types
cd4.hits <- names(which(tca.mdl.hannum.2.pvals.marg_cond[hits.joint,"CD4T.age"] == rowMins(tca.mdl.hannum.2.pvals.marg_cond[hits.joint,])))
sprintf("Detected %d associations using a joint test, %d associations using a marginal conditional test, and %d associations in CD4 cells using a marginal conditional test.",
sum(tca.mdl.hannum.2.pvals.joint <= 0.05/nrow(hannum$X)),
sum(tca.mdl.hannum.2.pvals.marg_cond <= 0.05/(nrow(hannum$X)*nrow(hannum$W))),
sum(tca.mdl.hannum.2.pvals.marg_cond[,"CD4T.age"] <= 0.05/(nrow(hannum$X)*nrow(hannum$W))))
# Replicate the CD4 hits in the Liu purified CD4 data
cd4.hits.liu.pvals <- unlist(lapply(1:length(cd4.hits),
function(x) summary(lm(y~., data.frame(y=liu.cd4$X[cd4.hits[x],], liu.cd4$cov)))$coefficients["age","Pr(>|t|)"]))
# Plot adjusted methylation (i.e. adjusted for covariates) as a function of age in all four replicated CD4 associations
dfs <- vector("list", length = 4)
for (i in 1:4){
r <- scale(residuals(lm(y~., data.frame(y=liu.cd4$X[cd4.hits[i],], liu.cd4$cov[,2:ncol(liu.cd4$cov)]))))
dfs[[i]] <- data.frame(x = liu.cd4$cov[,"age"], y = r)
}
plot_scatter(dfs = dfs,
ggarrange.ncol = 2,
ggarrange.nrow = 2,
xlab = "Age",
ylab = "Adjusted methylation level",
titles = paste("CD4 methylation in ",cd4.hits,sep=""))
# Verify that p-values are calibrated in the Liu data
liu.cd4.regression.pvals <- unlist(lapply(1:nrow(liu.cd4$X),
function(x) summary(lm(y~.,data.frame(y = liu.cd4$X[x,],liu.cd4$cov)))$coefficients["age","Pr(>|t|)"]))
plot_qq(list(liu.cd4.regression.pvals), labels = "Linear regression (sorted CD4)")
## Experiment #2: detecting differential methylation with infant arousal in granulocytes; working under the assumption Y|X (i.e. methylation levels affect or mediate components affecting infant arousal).
# Apply the TCA model to the Paquette data
tca.mdl.paquette <- tca(X = paquette$X,
W = paquette$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE)
# Run tcareg with a joint test and extract p-values; generate a qq-plot
C3_names <- c("gender","gestational_age","batch1","batch2")
tcareg.mdl.paquette.joint <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "joint")
plot_qq(list(tcareg.mdl.paquette.joint$pvals), labels = "Joint test with infant arousal")
# Since there is an inflation in the qq-plot we try to re-estimate the cell-type proportions
# First, run TCA using only the reference sites for getting a new estimate of W
tca.mdl.paquette.refit_W <- tca(X = paquette$X.ref_cpgs,
W = paquette$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE,
refit_W = TRUE,
refit_W.features = rownames(paquette$X.ref_cpgs))
# Use the updated estimate of W in a new execution of TCA on the data
tca.mdl.paquette.2 <- tca(X = paquette$X,
W = tca.mdl.paquette.refit_W$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE)
# Run tcareg with a joint test and extract p-values; generate a qq-plot
tcareg.mdl.paquette.joint.2 <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette.2,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "joint")
plot_qq(list(tcareg.mdl.paquette.joint.2$pvals), labels = "Joint test with infant arousal")
# Run tcareg with a marginal conditional test; generate qq plots
tcareg.mdl.paquette.2.marg_cond <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette.2,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "marginal_conditional")
plot_qq(split(tcareg.mdl.paquette.2.marg_cond$pvals,rep(1:ncol(paquette$W), each = nrow(paquette$X))),
labels = paste(colnames(paquette$W)," marginal conditional with arousal", sep=""),
ggarrange.nrow = 2,
ggarrange.ncol = 2)
# Exctract the hit we got from the joint test
hit.joint <- rownames(paquette$X)[which(tcareg.mdl.paquette.joint.2$pvals < 0.05/nrow(paquette$X))]
# The p-values of the marginal conditional test in the hit we found with the joint test suggest that the association is in granulocytes
tcareg.mdl.paquette.2.marg_cond$pvals[hit.joint,]
# Extract from the tca model the part that is relevant for the hit found
tcasub.mdl.paquette.2 <- tcasub(tca.mdl = tca.mdl.paquette.2,
features = hit.joint)
# Calculate cell-type-specific methylation for the samples in the detected CpG
tensor.mdl.hit <- tensor(tca.mdl = tcasub.mdl.paquette.2,
X = paquette$X[hit.joint,,drop=F])
# Plot the estimated granulocyte-specific methylation with arasual; adjust methylation levels to covariates
# Compare with using the bulk data
r.gran <- scale(residuals(lm(y~., data.frame(y=tensor.mdl.hit[[2]][1,], paquette$cov[,setdiff(1:ncol(paquette$cov),3)]))))
df.gran <- data.frame(y = r.gran, x = paquette$cov[,"arousal"])
# The bulk data should be further adjusted for cell-type composition (i.e. tca.mdl.paquette.2$W)
r.bulk <- scale(residuals(lm(y~., data.frame(y=paquette$X[hit.joint,], tca.mdl.paquette.2$W, paquette$cov[,setdiff(1:ncol(paquette$cov),3)]))))
df.bulk <- data.frame(y = r.bulk, x = paquette$cov[,"arousal"])
plot_scatter(dfs = list(df.bulk, df.gran),
ggarrange.ncol = 2,
ggarrange.nrow = 1,
xlab = "Infant arousal",
ylab = "Adjusted methylation level",
titles = paste(c("Whole-blood methylation in ", "Granulocyte methylation in "),hit.joint,sep=""))
}
cozygene/TCA documentation built on Feb. 18, 2021, 1:17 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library("TCA")
library("ggplot2")
library("ggpubr")
library("pracma")
library("matrixStats")
prep_data <- function(data_path){
file_name1 <- paste(data_path,.Platform$file.sep,"hannum.chr22.RData",sep="")
file_name2 <- paste(data_path,.Platform$file.sep,"liu.cd4.chr22.RData",sep="")
file_name3 <- paste(data_path,.Platform$file.sep,"paquette.chr22.RData",sep="")
file_names <- c(file_name1, file_name2, file_name3)
if(any(!file.exists(file_name1, file_name2, file_name3))){
library("GEOquery")
library("data.table")
library("EpiDISH")
}
if (!file.exists(file_name1)){
# Download the Hannum et al. data
gse <- GEOquery::getGEO("GSE40279", destdir = data_path, GSEMatrix = TRUE)
# Extract methylation data
X.hannum <- Biobase::exprs(gse[[1]])
# covariates; note that there's also 'ethnicity' covariate in the data, however, it is perfectly captured by the plateinformation
plate.hannum <- as.numeric(as.factor((Biobase::pData(gse[[1]])[["plate:ch1"]])))
cov.hannum <- cbind(as.numeric((Biobase::pData(gse[[1]])[["age (y):ch1"]])), as.numeric(as.factor((Biobase::pData(gse[[1]])[["gender:ch1"]]))), indicator_vars(plate.hannum))
rownames(cov.hannum) <- colnames(X.hannum)
colnames(cov.hannum) <- c("age", "gender", "plate1", "plate2", "plate3", "plate4", "plate5", "plate6", "plate7", "plate8")
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.hannum <- low_var_pcs(X.hannum, rank = 10, p = 1000)
cov.hannum <- cbind(cov.hannum, low_var_pcs.hannum)
# estimate cell-type proportions usign a reference-based approach
ref <- as.matrix(EpiDISH::centDHSbloodDMC.m[,c("Neutro","Eosino","CD4T","CD8T","Mono","B","NK")])
W.hannum <- EpiDISH::epidish(X.hannum, ref)$estF
# merge Neutro and Eosino
W.hannum.gran <- W.hannum[,"Neutro",drop=F] + W.hannum[,"Eosino",drop=F]
colnames(W.hannum.gran) <- "Gran"
W.hannum <- cbind(W.hannum.gran,W.hannum[,c("CD4T","CD8T","Mono","B","NK")])
# remove polymorphic sites, cross-reactive sites, and non-autosomal sites according to Chen et al.; in addition, remove low variance sites, as those are unlikely to demonstrate biological signal.
X.hannum.processed <- filter_data(X.hannum)
# keep only sites on chromosome 22 (so that tutorial can run quickly)
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- map[rownames(X.hannum.processed),1]
chr22sites <- which(chrs == 22)
X.hannum.processed.22 <- X.hannum.processed[chr22sites,]
hannum <- list(X = X.hannum.processed.22, cov = cov.hannum, W = W.hannum)
save(hannum, file = file_name1, compress = "bzip2")
rm(hannum, gse, X.hannum, X.hannum.processed)
file.remove(paste(data_path,.Platform$file.sep,"GPL13534.soft",sep=""))
file.remove(paste(data_path,.Platform$file.sep,"GSE40279_series_matrix.txt.gz",sep=""))
}
if (!file.exists(file_name2)){
# Download the Liu et al. CD4 data
# covariates
gse <- GEOquery::getGEO("GSE56581", destdir = data_path, GSEMatrix = TRUE)
liu.cd4.age <- as.matrix(as.numeric(Biobase::pData(gse[[1]])[["age (yrs):ch1"]]))
colnames(liu.cd4.age) <- "age"
# methylation data
liu.cd4.methfile <- paste(data_path,.Platform$file.sep,"GSE56581_methylome_normalized.txt",sep="")
download.file("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE56nnn/GSE56581/suppl/GSE56581_methylome_normalized.txt.gz", paste(liu.cd4.methfile,".gz",sep=""))
gunzip(paste(liu.cd4.methfile,".gz",sep=""))
X.liu.cd4 <- fread(liu.cd4.methfile)
X.liu.cd4.rownames <- X.liu.cd4$ID_REF
X.liu.cd4 <- as.matrix(X.liu.cd4[,2:ncol(X.liu.cd4)])
rownames(X.liu.cd4) <- X.liu.cd4.rownames
X.liu.cd4 <- X.liu.cd4[,seq(1,ncol(X.liu.cd4),2)]
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.liu.cd4 <- low_var_pcs(X.liu.cd4, rank = 10)
cov.liu.cd4 <- cbind(liu.cd4.age, low_var_pcs.liu.cd4)
# keep only sites on chromosome 22 (so that tutorial can run quickly)
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- map[rownames(X.liu.cd4),1]
chr22sites <- which(chrs == 22)
X.liu.cd4.22 <- X.liu.cd4[chr22sites,]
liu.cd4 <- list(X = X.liu.cd4.22, cov = cov.liu.cd4)
save(liu.cd4, file = file_name2)
rm(liu.cd4, gse, X.liu.cd4.22, X.liu.cd4)
file.remove(paste(data_path,.Platform$file.sep,"GSE56581_methylome_normalized.txt",sep=""))
file.remove(paste(data_path,.Platform$file.sep,"GSE56581_series_matrix.txt.gz",sep=""))
}
if (!file.exists(file_name3)){
# Paquette et al. (Epigenetics 2016)
gse <- GEOquery::getGEO("GSE75248", destdir = data_path, GSEMatrix = TRUE)
# methylation data
X.paquette <- Biobase::exprs(gse[[1]])
# Calculate principal components from control probes that reflect technical variability; since we are not working with IDAT files we don't have "real" control probes - instead, we use low variance sites that are not expected to capture any true biological variability.
low_var_pcs.paquette <- low_var_pcs(X.paquette, rank = 10, p = 1000)
# remove polymorphic sites, cross-reactive sites, and non-autosomal sites according to Chen et al.; in addition, remove low variance sites, as those are unlikely to demonstrate biological signal.
X.paquette.processed <- filter_data(X.paquette)
# keep only sites on chromosome 22 (so that tutorial can run quickly); Paquette et al. reported some signal in chromosome 16
map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
chrs <- as.character(map[rownames(X.paquette.processed),1])
chr22sites <- which(chrs == "22")
X.paquette.processed.22 <- X.paquette.processed[chr22sites,]
# Extract covaraites
gestational_age <-as.numeric(unlist(lapply(Biobase::pData(gse[[1]])[["gestational age:ch1"]], function(x) strsplit(x,";")[[1]][1])))
arousal <- as.numeric(Biobase::pData(gse[[1]])[["arousal:ch1"]])
batch <- as.numeric(as.factor((Biobase::pData(gse[[1]])[["batch:ch1"]])))
batch[is.na(batch)] <- 3
batch <- indicator_vars(batch)
# since gender information is not available on the GEO record, get gender information that was inferred from the IDAR files based on X chromosome methylation
gender <- read.table("https://raw.githubusercontent.com/cozygene/TCA/master/vignettes/gender.paquette.txt", row.names = 1)
cov.paquette <- data.frame(gender,gestational_age,arousal,batch)
colnames(cov.paquette) <- c("gender","gestational_age","arousal","batch1","batch2")
rownames(cov.paquette) <- colnames(X.paquette.processed.22)
# remove na values
keep <- setdiff(1:nrow(cov.paquette) ,which(rowSums(is.na(cov.paquette)) > 0))
cov.paquette <- cbind(cov.paquette[keep,], low_var_pcs.paquette[keep,])
X.paquette.processed.22 <- X.paquette.processed.22[,keep]
# Load cell-type proportions usign a reference-based approach; these were estimated using the raw IDAT files of the data (available on GEO) with the package ENmix.
W.paquette <- read.table("https://raw.githubusercontent.com/cozygene/TCA/master/vignettes/W.paquette.txt", row.names = 1)
W.paquette <- W.paquette[,2:ncol(W.paquette)]
W.paquette[W.paquette < 0] <- 0
# remove lowly abundant cell types
W.paquette <- W.paquette[keep, colMeans(W.paquette)>=0.01]
rownames(W.paquette) <- rownames(cov.paquette)
# normalize cell type proportions (output from the ENmix package didn't require the proportions of an individual to sum up to 1)
W.paquette <- W.paquette/t(repmat(rowSums(W.paquette),ncol(W.paquette),1))
# get the set of cord blood reference CpGs that were used for estimating cell-type proportions; we will use these CpGs for re-estimating W
library("FlowSorted.CordBlood.450k")
ref.cord <- get("FlowSorted.CordBlood.450k")
ref.cord <- preprocessRaw(ref.cord)
ref.cord.cpgs <- rownames(FlowSorted.CordBlood.450k.ModelPars)
ref.cord.cpgs.intersect <- intersect(ref.cord.cpgs, rownames(X.paquette))
paquette <- list(X = X.paquette.processed.22, cov = cov.paquette, W = W.paquette, X.ref_cpgs = X.paquette[ref.cord.cpgs.intersect,keep])
save(paquette, file = file_name3)
rm(paquette, gse, X.paquette, X.paquette.processed, X.paquette.processed.22)
file.remove(paste(data_path,.Platform$file.sep,"GSE75248_series_matrix.txt.gz",sep=""))
## note - the files W.paquette.txt and gender.paquette.txt were generated using the original IDAT files of the Paquette data.
# Gender information was inferred from methylation levels frmo the X chromosome as it is not available on Tthe GEO record.
# Code attched below; in order to run it, download the IDAT files from GEO accession number GSE75248 and set 'path' in 'readidat' below to the location of the files.
# require("ENmix")
# rgdat <- readidat(path = ".", manifestfile=NULL, recursive = FALSE, verbose = TRUE)
# meth <- getmeth(rgdat); rm(rgdat)
# W <- estimateCellProp(meth, refdata="FlowSorted.CordBlood.450k", cellTypes=NULL, nonnegative = TRUE, nProbes=50, normalize=TRUE, refplot=FALSE)
# write.table(W, file = "W.paquette.txt", quote = FALSE, sep = " ")
# beta <- getB(meth)
# map <- read.table("https://github.com/cozygene/glint/blob/master/utils/assets/HumanMethylationSites?raw=true",sep=",",row.names = 1)
# chrs <- as.character(map[rownames(beta),1])
# sites.X <- which(chrs == "X")
# sex_cpgs.pca.paquette <- prcomp(t(beta[sites.X,]), center=TRUE, scale=TRUE, rank = 2)
# write.table(sex_cpgs.pca.paquette$x[,1:2], file = "sex_cpgs.pca.paquette.txt", quote = FALSE, sep = " ")
# plot(sex_cpgs.pca.paquette$x[,1],sex_cpgs.pca.paquette$x[,2]) # shows a separation between males and females; note that the males/females ratio here perfectly match the one in the Paquette et al. paper.
# df <-data.frame(as.numeric(x[,1] > 0))
# colnames(df) <-"gender"; rownames(df) <- rownames(x)
# write.table(df,file = "gender.paquette.txt", quote = FALSE, sep = " ")
}
return(file_names)
}
indicator_vars <- function(x){
x.indicators <- matrix(0,length(x),length(unique(x))-1)
counter <- 1
u <- unique(x)
for (i in setdiff(u,u[length(u)])){
x.indicators[,counter] <- as.numeric(x == i)
counter <- counter + 1
}
return(x.indicators)
}
low_var_pcs <- function(X, rank = 10, p = 1000){
site.variances <- matrixStats::rowVars(X)
names(site.variances) <- rownames(X)
low.var.sites <- names(head(sort(site.variances), p))
low.var.pca <- prcomp(t(X[low.var.sites,]), center=TRUE, scale=TRUE, rank = rank)
return(low.var.pca$x)
}
filter_data <- function(X, var_th = 0.0001){
nonspecific_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/nonspecific_probes.txt")[,1]
XY_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/HumanMethylationSites_X_Y.txt")[,1]
polymorphic_probes <- read.table("https://raw.githubusercontent.com/cozygene/glint/master/parsers/assets/polymorphic_cpgs.txt")[,1]
# remove sites with very low variance that are unlikely to exhibit biological signal
site.variances <- matrixStats::rowVars(X)
names(site.variances) <- rownames(X)
low_var_sites <- names(which(site.variances < var_th))
exclude <- union(nonspecific_probes,union(XY_probes,union(low_var_sites,polymorphic_probes)))
return(X[setdiff(rownames(X),exclude),])
}
plot_qq <- function(pvals, labels, ggarrange.nrow = 1, ggarrange.ncol = 1, alpha = 0.05, experiment_wide_line = TRUE){
significance_th <- list(alpha/length(pvals[[1]]))
if(length(pvals)-1) significance_th[[2]] <- alpha/(length(pvals)*length(pvals[[1]]))
qqplots <- lapply(1:length(pvals), function(p){
df <- data.frame(pvals.obs = -log10(sort(pvals[[p]])), pvals.exp = -log10(sort((1:length(pvals[[1]]))/length(pvals[[1]]))));
qqplot <- ggplot(df, aes(x = pvals.exp, y = pvals.obs)) +
stat_binhex(geom = "point", bins=1000, size=1) +
geom_abline() +
ggtitle(labels[p]) +
xlab(expression(Expected~-log[10](P))) + ylab(expression(Observed~-log[10](P))) +
theme_bw() +
guides(fill="none") +
geom_hline(yintercept=-log10(significance_th[[1]]), linetype="dashed", color = "red", size=1)
if (length(significance_th)-1 & experiment_wide_line) qqplot <- qqplot + geom_hline(yintercept=-log10(significance_th[[2]]), linetype="dashed", color = "red", size=0.5)
return(qqplot)
})
ggarrange(plotlist = qqplots, ncol = ggarrange.ncol, nrow = ggarrange.nrow)
}
plot_scatter <- function(dfs, ggarrange.ncol, ggarrange.nrow, xlab, ylab, titles){
plots <- vector("list", length = length(dfs))
for (i in 1:length(dfs)){
df <- data.frame(y = dfs[[i]]$y, x = dfs[[i]]$x)
plots[[i]] <- ggplot(df, aes(x = x, y = y)) +
geom_point(alpha = 0.6) +
geom_smooth(method=lm) +
stat_cor(method = "pearson", colour = "blue") +
xlab(xlab) +
ylab(ylab) +
ggtitle(titles[i])
}
ggarrange(plotlist = plots, ncol = ggarrange.ncol, nrow = ggarrange.nrow)
}
if(FALSE){
# Set a path for storing data
data_path <- "./"
# Load the data
filenames <- prep_data(data_path)
for (filename in filenames) load(filename)
## Experiment #1: detecting CD4 differential methylation with age; working under the assumption X|Y (i.e. age affects methylation levels)
# Apply the TCA model to the Hannum whole-blood data
tca.mdl.hannum <- tca(X = hannum$X,
W = hannum$W,
C1 = hannum$cov[,c("gender","age")],
C2 = hannum$cov[,3:ncol(hannum$cov)])
# Extract p-values of a joint test
tca.mdl.hannum.pvals.joint <- tca.mdl.hannum$gammas_hat_pvals.joint[,"age"]
# Extract p-values for each cell type, under a marginal conditional test
tca.mdl.hannum.pvals.marg_cond <- tca.mdl.hannum$gammas_hat_pvals[,paste(colnames(hannum$W),".age",sep="")]
# qq-plots - for the p-values of the joint test, and for the p-values in CD4, under a marginal conditional test
plot_qq(list(tca.mdl.hannum.pvals.joint, tca.mdl.hannum.pvals.marg_cond[,"CD4T.age"]),
labels = c("Joint test with age", "CD4 marginal conditional test with age"),
ggarrange.nrow = 1,
ggarrange.ncol = 2,
experiment_wide_line = FALSE)
# Run ReFACTor to capture more of the variation of cell-type composition
refactor.mdl.hannum <- refactor(X = hannum$X,
k = 6,
C = hannum$cov[,3:ncol(hannum$cov)])
# Rerun TCA, this time include the ReFACTor components as additional covariates
tca.mdl.hannum.2 <- tca(X = hannum$X,
W = hannum$W,
C1 = hannum$cov[,c("gender","age")],
C2 = cbind(hannum$cov[,3:ncol(hannum$cov)],refactor.mdl.hannum$scores))
# Extract the updated p-values of a joint test
tca.mdl.hannum.2.pvals.joint <- tca.mdl.hannum.2$gammas_hat_pvals.joint[,"age"]
# Extract the updated marginal conditional p-values
tca.mdl.hannum.2.pvals.marg_cond <- tca.mdl.hannum.2$gammas_hat_pvals[,paste(colnames(hannum$W),".age",sep="")]
# qq-plots - for the p-values of the joint test, and for the p-values in CD4, under a marginal conditional test
plot_qq(list(tca.mdl.hannum.2.pvals.joint, tca.mdl.hannum.2.pvals.marg_cond[,"CD4T.age"]),
labels = c("Joint test with age", "CD4 marginal conditional test with age"),
ggarrange.nrow = 1,
ggarrange.ncol = 2,
experiment_wide_line = FALSE)
# Extract the hits based on the joint test
hits.joint <- names(which(tca.mdl.hannum.2.pvals.joint < 0.05/nrow(hannum$X)))
# Extract the hits from hits.joint where CD4 cells demonstrate the lowest p-value across all cell types
cd4.hits <- names(which(tca.mdl.hannum.2.pvals.marg_cond[hits.joint,"CD4T.age"] == rowMins(tca.mdl.hannum.2.pvals.marg_cond[hits.joint,])))
sprintf("Detected %d associations using a joint test, %d associations using a marginal conditional test, and %d associations in CD4 cells using a marginal conditional test.",
sum(tca.mdl.hannum.2.pvals.joint <= 0.05/nrow(hannum$X)),
sum(tca.mdl.hannum.2.pvals.marg_cond <= 0.05/(nrow(hannum$X)*nrow(hannum$W))),
sum(tca.mdl.hannum.2.pvals.marg_cond[,"CD4T.age"] <= 0.05/(nrow(hannum$X)*nrow(hannum$W))))
# Replicate the CD4 hits in the Liu purified CD4 data
cd4.hits.liu.pvals <- unlist(lapply(1:length(cd4.hits),
function(x) summary(lm(y~., data.frame(y=liu.cd4$X[cd4.hits[x],], liu.cd4$cov)))$coefficients["age","Pr(>|t|)"]))
# Plot adjusted methylation (i.e. adjusted for covariates) as a function of age in all four replicated CD4 associations
dfs <- vector("list", length = 4)
for (i in 1:4){
r <- scale(residuals(lm(y~., data.frame(y=liu.cd4$X[cd4.hits[i],], liu.cd4$cov[,2:ncol(liu.cd4$cov)]))))
dfs[[i]] <- data.frame(x = liu.cd4$cov[,"age"], y = r)
}
plot_scatter(dfs = dfs,
ggarrange.ncol = 2,
ggarrange.nrow = 2,
xlab = "Age",
ylab = "Adjusted methylation level",
titles = paste("CD4 methylation in ",cd4.hits,sep=""))
# Verify that p-values are calibrated in the Liu data
liu.cd4.regression.pvals <- unlist(lapply(1:nrow(liu.cd4$X),
function(x) summary(lm(y~.,data.frame(y = liu.cd4$X[x,],liu.cd4$cov)))$coefficients["age","Pr(>|t|)"]))
plot_qq(list(liu.cd4.regression.pvals), labels = "Linear regression (sorted CD4)")
## Experiment #2: detecting differential methylation with infant arousal in granulocytes; working under the assumption Y|X (i.e. methylation levels affect or mediate components affecting infant arousal).
# Apply the TCA model to the Paquette data
tca.mdl.paquette <- tca(X = paquette$X,
W = paquette$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE)
# Run tcareg with a joint test and extract p-values; generate a qq-plot
C3_names <- c("gender","gestational_age","batch1","batch2")
tcareg.mdl.paquette.joint <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "joint")
plot_qq(list(tcareg.mdl.paquette.joint$pvals), labels = "Joint test with infant arousal")
# Since there is an inflation in the qq-plot we try to re-estimate the cell-type proportions
# First, run TCA using only the reference sites for getting a new estimate of W
tca.mdl.paquette.refit_W <- tca(X = paquette$X.ref_cpgs,
W = paquette$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE,
refit_W = TRUE,
refit_W.features = rownames(paquette$X.ref_cpgs))
# Use the updated estimate of W in a new execution of TCA on the data
tca.mdl.paquette.2 <- tca(X = paquette$X,
W = tca.mdl.paquette.refit_W$W,
C1 = paquette$cov[,c("gender","gestational_age")],
C2 = paquette$cov[,4:ncol(paquette$cov)],
constrain_mu = TRUE)
# Run tcareg with a joint test and extract p-values; generate a qq-plot
tcareg.mdl.paquette.joint.2 <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette.2,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "joint")
plot_qq(list(tcareg.mdl.paquette.joint.2$pvals), labels = "Joint test with infant arousal")
# Run tcareg with a marginal conditional test; generate qq plots
tcareg.mdl.paquette.2.marg_cond <- tcareg(X = paquette$X,
tca.mdl = tca.mdl.paquette.2,
y = paquette$cov[,"arousal",drop=F],
C3 = paquette$cov[,C3_names],
test = "marginal_conditional")
plot_qq(split(tcareg.mdl.paquette.2.marg_cond$pvals,rep(1:ncol(paquette$W), each = nrow(paquette$X))),
labels = paste(colnames(paquette$W)," marginal conditional with arousal", sep=""),
ggarrange.nrow = 2,
ggarrange.ncol = 2)
# Exctract the hit we got from the joint test
hit.joint <- rownames(paquette$X)[which(tcareg.mdl.paquette.joint.2$pvals < 0.05/nrow(paquette$X))]
# The p-values of the marginal conditional test in the hit we found with the joint test suggest that the association is in granulocytes
tcareg.mdl.paquette.2.marg_cond$pvals[hit.joint,]
# Extract from the tca model the part that is relevant for the hit found
tcasub.mdl.paquette.2 <- tcasub(tca.mdl = tca.mdl.paquette.2,
features = hit.joint)
# Calculate cell-type-specific methylation for the samples in the detected CpG
tensor.mdl.hit <- tensor(tca.mdl = tcasub.mdl.paquette.2,
X = paquette$X[hit.joint,,drop=F])
# Plot the estimated granulocyte-specific methylation with arasual; adjust methylation levels to covariates
# Compare with using the bulk data
r.gran <- scale(residuals(lm(y~., data.frame(y=tensor.mdl.hit[[2]][1,], paquette$cov[,setdiff(1:ncol(paquette$cov),3)]))))
df.gran <- data.frame(y = r.gran, x = paquette$cov[,"arousal"])
# The bulk data should be further adjusted for cell-type composition (i.e. tca.mdl.paquette.2$W)
r.bulk <- scale(residuals(lm(y~., data.frame(y=paquette$X[hit.joint,], tca.mdl.paquette.2$W, paquette$cov[,setdiff(1:ncol(paquette$cov),3)]))))
df.bulk <- data.frame(y = r.bulk, x = paquette$cov[,"arousal"])
plot_scatter(dfs = list(df.bulk, df.gran),
ggarrange.ncol = 2,
ggarrange.nrow = 1,
xlab = "Infant arousal",
ylab = "Adjusted methylation level",
titles = paste(c("Whole-blood methylation in ", "Granulocyte methylation in "),hit.joint,sep=""))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.