setup: DNA from ten individuals with known genotypes (e.g. each person
is either AA, AC, or CC at a particular locus) at 200 genetic loci are
pooled at equivalent concentrations. Now we have sequencing data from
this pool, such that at each of the 200 loci we have counts for each
observed nucleotide (e.g. 40 A and 25 C; total counts for each loci
vary randomly). All 200 loci vary in at least one individual (i.e.
there is no loci at which all 10 individuals are identically AA or all
identically CC). The inference problem is to estimate the actual
fraction of each individual contributing to the overall counts, using
information across all 200 loci to disambiguate the pooling.
approach: model the counts of one allele seen at each loci as
binomial-distributed with unknown mean mu; mu is the inner product
between the vector "fracs" of unknown true fractions of each
individual and the vector of known genotypes (expressed as a ratio -
0, 0.5 or 1). The vector of unknown true fractions is dirichlet with
10 parameters equal to 1 (or 10, if I wish to make the prior more
informative of the expected frequencies).
BUGS code:
model {
fracs[1:10] ~ ddirch(alpha[1:10])
# for (k in 1:10) {
# fracs[k] <- delta[k]/sum(delta[1:10]);
# delta[k] ~ dgamma(alpha[k], 1)
# }
for (i in 1:N) {
mu[i] <- inprod(fracs[1:10], ratio[i,1:10]);
snpA[i] ~ dbin(mu[i], total[i]);
}
}
The commented-out bit there is what I have to do to get this code to
even work without an error, since otherwise I get "cannot find an
updater for fracs". With the delta/dgamma version, everything works,
but autocorrelation is very bad.
Suggestions?
Thanks,
-Aaron
P.S. here's a reduced set of data you can try this with:
list(alpha=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), N=40, ratio=
structure(.Data= c(0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0.5, 0.5, 0, 0.5,
0.5, 0.5, 0.5, 1, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0.5, 0.5, 1, 0.5, 0,
0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0, 1, 1, 0, 0.5, 0, 0, 0.5, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0.5, 0.5, 1, 0.5, 0.5, 1, 0.5,
0.5, 0.5, 0.5, 0.5, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0.5, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 1, 1, 0.5, 1, 0.5, 1, 0.5, 1, 0.5,
0, 1, 1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 1, 0, 0.5, 0, 0.5,
0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0.5, 1, 1, 0.5, 0.5, 0.5, 0.5, 1,
0.5, 1, 1, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 1, 0, 0.5, 1, 1, 1, 0.5,
1, 0.5, 1, 0.5, 1, 0.5, 0, 0.5, 0.5, 0.5, 0, 0.5, 0, 0.5, 0, 0.5, 1,
0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0.5, 1, 0.5, 1, 0.5, 0.5, 0.5, 0, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0.5, 0.5, 1, 1, 0, 1, 1,
0, 1, 1, 0.5, 0.5, 1, 1, 0, 1, 0, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 0, 1, 0,
1, 1, 0.5, 0.5, 1, 1, 0, 1, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 1, 0.5, 1,
0.5, 1, 0.5, 1, 0.5, 0.5, 1, 0, 0.5, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 1,
0.5, 0, 0.5, 0, 0.5, 0, 0.5, 0.5, 0, 1, 0.5, 1, 0.5, 1, 1, 1, 1, 0.5,
1, 0, 0.5, 0, 0.5, 0, 0.5, 0, 0.5, 0.5, 0.5, 1, 0.5, 0, 0.5, 0, 0, 0,
0, 0.5, 0, 1, 0.5, 1, 0.5, 1, 0.5, 1, 0.5, 0.5, 1, 0, 0.5, 1, 0.5, 1,
0.5, 1, 0.5, 0.5, 1, 0, 0.5, 0, 0.5, 0, 0.5, 0, 0.5, 0.5, 0, 1, 0.5,
1, 0.5, 1, 1, 1, 1, 0.5, 1, 0, 0.5), .Dim=c(40, 10)), snpA=c(7, 38,
23, 72, 12, 68, 51, 20, 93, 0, 1, 40, 15, 0, 10, 22, 36, 48, 45, 29,
28, 94, 26, 0, 64, 24, 30, 18, 48, 13, 20, 31, 23, 23, 31, 10, 39, 44,
24, 41), total=c(85, 93, 69, 106, 103, 68, 91, 66, 93, 85, 67, 60, 20,
30, 40, 56, 46, 78, 66, 69, 30, 94, 50, 118, 65, 44, 39, 75, 61, 40,
41, 94, 79, 30, 67, 41, 65, 68, 77, 55))
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|