library(rjags)

# Load variant read counts from amplicon deep sequencing experiment and the exome
# sequencing experiment
deep.seq <- read.delim('readcounts.txt')

# Calculate probabilities of the null hypothesis with JAGS for primary
load.module('mix')
run <- jags.model('same_overdispersion_model.bugs',
         data=with(deep.seq, list(x=normal_variant,
                                  y=primary_variant,
                                  nx=normal_total,
                                  ny=primary_total,
                                  k=nrow(deep.seq))))
samples <- jags.samples(run, c('tau', 'is.null'), 1000)
prim.p.is.null <- apply(samples$is.null, 1, mean)
names(prim.p.is.null) <- deep.seq$gene

# Calculate probabilities of the null hypothesis with JAGS for met
run2<- jags.model('different_overdispersion_model.bugs',
         data=with(deep.seq, list(x=normal_variant,
                                  y=met_variant,
                                  nx=normal_total,
                                  ny=met_total,
                                  k=nrow(deep.seq))))
samples2<- jags.samples(run2, c('tau', 'tau2', 'is.null'), 1000)
met.p.is.null <- apply(samples2$is.null, 1, mean)
names(met.p.is.null) <- deep.seq$gene

# Write null probabilities to file
write.table(data.frame(gene=deep.seq$gene, prim_null_prob=prim.p.is.null, met_null_prob=met.p.is.null), "null_probabilities.txt", row.names=FALSE, quote=FALSE)
