Description Usage Arguments Value See Also Examples
Marginalises over the w parameter in a fitted epiallele model
1 | bed_marginal(Y, X, W)
|
Y |
an N by d numeric matrix of observed DNAm reads |
X |
an Q by d numeric matrix of inferred epialleles |
W |
an N by 1 numeric vector assigning reads to epiallels |
An N by d numerc matrix where where element [i,mu] is the probability that read i corresponds to epiallele mu.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | set.seed(1314)
# Generate 50 reads that are fully methylated
Y1 <- matrix(rep(c(1,1,1,1,1,1,1,1),50), nrow=50, ncol=8)
# Generate 30 reads that are half methylated
Y2 <- matrix(rep(c(1,1,1,1,0,0,0,0),30),nrow=30, ncol=8, byrow=TRUE)
# Generate 20 reads that are fully unmethylated
Y3 <- matrix(rep(c(0,0,0,0,0,0,0,0),20),nrow=20,ncol=8)
# Combine them
Y <- rbind(Y1,Y2,Y3)
# Add some noise by randomly regenerating 20 percent of values
Y[sample(1:800,80,replace=FALSE)] <- round(runif(80))
# Make a few data missing
Y[sample(1:800,5,replace=FALSE)] <- NA
# Infer which epialleles are present
fit <- bed_fit(Y)
# Plot AIC score
plot(fit$AIC[1:10], type='b', xlab='Q = number of epialleles',
ylab='AIC score')
# The fit$rho field contains the optimal number of epialleles
cat(' Optimal number of epialleles =', fit$rho,'\n')
#Marginalise over W posterior
w.margin <- bed_marginal(Y, fit$model[[fit$rho]]$X, fit$model[[fit$rho]]$W)
# Summing over the columns gives the total proportion of reads that are
# attributed to each epiallele
p <- apply(w.margin, 2, FUN='sum')/100
cat(' Inferred Epiallele [1] =',paste(fit$model[[3]]$X[1,],collapse=""),
', Proportion =', p[1],'\n',
'Inferred Epiallele [2] =',paste(fit$model[[3]]$X[2,],collapse=""),
', Proportion =', p[2],'\n',
'Inferred Epiallele [3] =',paste(fit$model[[3]]$X[3,],collapse="")
,', Proportion =', p[3])
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.