This is the WinBUGS (Spiegelhalter et al., 1999) code for the Bayesian spatially structured mixture (three components) model considered in the article of Broet and Richardson (Bioinformatics, 22(8): 911-918) but presented (for simplicity) for one chromosome. The code shown here is for this simplified version of the model but it can easily be modified to include clones on more than one chromosome, indexing the precisions in the random fields by chromosome. We used version 1.4 of WinBUGS to run this model.

model
{
for( i in 1 : N ) {Y[i] ~dnorm(mu[i], tauR[i]) # genomic sequence measure-
ments
mu[i] <- lambda[T[i]]
tauR[i] <- tau[T[i]]
p[i,1] <-exp(alpha1[i])/(exp(alpha1[i])+exp(alpha2[i])+exp(alpha3[i]))
p[i,2] <-exp(alpha2[i])/(exp(alpha1[i])+exp(alpha2[i])+exp(alpha3[i]))
p[i,3] <-exp(alpha3[i])/(exp(alpha1[i])+exp(alpha2[i])+exp(alpha3[i]))
T[i] ~dcat(p[i,1:3]) # latent variable (deletion/gain/modal)
T1[i] <-equals(T[i],1) ;T2[i] <-equals(T[i],2) ;T3[i] <-equals(T[i],3) # latent
variables}
# Random Fields specification
alpha1[1:N] ~car.normal(adj[], weights[], num[], tauC[1])
alpha2[1:N] ~car.normal(adj[], weights[], num[], tauC[2])
alpha3[1:N] ~car.normal(adj[], weights[], num[], tauC[3])
# weights specification
for(i in 1:1){weights[i] <- 1;adj[i] <- i+1 ;num[i]<- 1}
for(i in 2:(N-1)) {weights[2+(i-2)*2] <- 1;adj[2+(i-2)*2] <- i-1
weights[3+(i-2)*2] <- 1; adj[3+(i-2)*2] <- i+1;num[i] <- 2}
for(i in N:N){weights[(N-2)*2 + 2] <- 1;adj[(N-2)*2 + 2] <- i-1;num[i] <- 1}
# prior (precision/variance spatial .field)
sig2C[1] ~dgamma(0.01, 0.01)I(0.0001, )
sig2C[2] ~dgamma(0.01, 0.01)I(0.0001, )
sig2C[3] ~dgamma(0.01, 0.01)I(0.0001, )
tauC1[1] <- 1/sig2C[1]; tauC1[2] <- 1/sig2C[2]; tauC1[3] <- 1/sig2C[3]
# prior (range [a-b])
lambda[1] <- 0 # modal copy
lambda[3] ~dnorm(0, 1.0E-6)I(a,0.0 ) # loss copy
lambda[2] ~dnorm(0, 1.0E-6)I(0.0,b) # gain copy
# precision/variance (density)
tau[1] ~dgamma(0.1, 0.1)
tau[2]~dgamma(0.1, 0.1)
tau[3]~dgamma(0.1, 0.1)
sigma[1] <- 1/sqrt(tau[1]) ;sigma[2] <- 1/sqrt(tau[2]) ;sigma[3] <- 1/sqrt(tau[3])
}