bed_marginal: Marginal epiallele distribution

In james-e-barrett/bed: Bayesian Epiallele Detection

Description

Marginalises over the w parameter in a fitted epiallele model

Usage

bed_marginal(Y, X, W)

Arguments

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

Value

An N by d numerc matrix where where element [i,mu] is the probability that read i corresponds to epiallele mu.

See Also

bed_fit

Examples

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])

james-e-barrett/bed documentation built on May 18, 2019, 11:19 a.m.