Bayesian Manifold-Constrained-Prior Model for an Experiment to Locate Xce

TL;DR

This paper introduces a Bayesian manifold-constrained-prior model to analyze an experiment aimed at locating the Xce gene, accounting for biases in X-inactivation and using a novel physical model based on copy number variation.

Contribution

It presents a new Bayesian analysis framework with manifold-constrained priors and a physical model for predicting X-inactivation bias in genetic crosses.

Findings

The model effectively accounts for measurement biases and differences in gene expression precision.

Reparameterized slice-sampling handles complex constrained priors successfully.

The physical model links Xce alleles to copy number variation, explaining their differences.

Abstract

We propose an analysis for a novel experiment intended to locate the genetic locus Xce (X-chromosome controlling element), which biases the stochastic process of X-inactivation in the mouse. X-inactivation bias is a phenomenon where cells in the embryo randomly choose one parental chromosome to inactivate, but show an average bias towards one parental strain. Measurement of allele-specific gene-expression through pyrosequencing was conducted on mouse crosses of an uncharacterized parent with known carriers. Our Bayesian analysis is suitable for this adaptive experimental design, accounting for the biases and differences in precision among genes. Model identifiability is facilitated by priors constrained to a manifold. We show that reparameterized slice-sampling can suitably tackle a general class of constrained priors. We demonstrate a physical model, based upon a "weighted-coin" hypothesis, that predicts X-inactivation ratios in untested crosses. This model suggests that Xce alleles differ due to a process known as copy number variation, where stronger Xce alleles are shorter sequences.

Tables

Table 1. Table 1: Information on common inbred mouse strains used in this experiment.

Strain	Abbrev.	details	Sub-species	derived	Xce allele
A/J	AJ	Albino (hearing loss)	domesticus	lab	$a$
129S1/SvlmJ	129	Brown (testicular cancer)	dom.	lab	$a$
C57BL/6J	B6	“Black-six”	dom.	lab	$b$
ALS/LtJ	ALS	Lou Gehrig’s disease	dom.	lab	$b$
CAST/EiJ	CAST		castaneus	wild	$c$
WSB/EiJ	WSB	“Watkin Star”, aka “White Spot”	domesticus	wild	??
WLA/Pas	WLA	Malaria Resistant	dom.	wild	??
LEWES/EiJ	LEWES	Caught in Lewes, DE	dom.	wild	??
PWK/PhJ	PWK		musculus	wild	??
SJL/J	SJL	Albino (sight loss)	dom.	lab	??
TIRANO/EiJ	TIRANO	Poschiavinus Valley, Tirano, Italy	dom.	wild	??
ZALENDE/Ei	ZALENDE	Zalende, Switzerland	dom.	wild	??
DDK/Pas	DDK	Fertility loss, Japanese	dom.	lab	$a$
PERA/EiJ	PERA	Rimac Valley, Peru	dom.	wild	??
PANCEVO/EiJ	PANCEVO	Pancevo, Serbia	spicilegus	wild	??

Table 2. Table 2: Example section of pyrosequencing maternal expression ratios dataset.

pup	cross	dam	sire	$\begin{array}{c}\text{Brain}\\ \text{Ddx26b}\end{array}$	$\begin{array}{c}\text{Brain}\\ \text{Rragb}\end{array}$	$\begin{array}{c}\text{Kidney}\\ \text{Ddx26b}\end{array}$	$\begin{array}{c}\text{Kidney}\\ \text{Rragb}\end{array}$	…
3-7	1Wl	129S1-3	WLA-3	0.234	NA	0.183	NA	…
5-15	1Wl	129S1-5	WLA-5	0.316	NA	NA	NA
…	…
2-1	1Wl	129S1-2	WLA-7	0.606	NA	0.483	NA
7-4	1Wl	129S1-7	WLA-7	0.648	NA	NA	NA
16-3	AlAj	ALS-16	AJ-873	0.456	.501	0.447	.51
16-5	AlAj	ALS-16	AJ-873	0.574	.432	0.532	.59	…
…

Table 3. Table 5: Simulating 1000 datasets from a fixed μ → → → 𝜇 \vec{\overset{{}_{\shortrightarrow}}{\mu}} vector in model

	$\text{True }\mu$	0.25	0.45	0.5	0.65	0.75
Bayes	$\text{Mean (Bias)}\hat{\mu }-\mu$	-0.003	-0.001	-0.001	0.001	0
Method	$\text{Mean }\sqrt{{\left(\hat{\mu }-\mu \right)}^2}$	0.015	0.02	0.022	0.025	0.027
	$95$% HPD Width	0.056	0.063	0.063	0.061	0.056
	$95$% HPD Coverage	0.974	0.972	0.962	0.962	0.969
Sample	$\text{Mean (Bias)}\hat{\mu }-\mu$	0.022	0.003	-0.002	-0.016	-0.025
Mean	$\text{Mean }\sqrt{{\left(\hat{\mu }-\mu \right)}^2}$	0.026	0.015	0.016	0.022	0.028
	$95$% CI Width	0.054	0.061	0.06	0.056	0.051
	$95$% CI Coverage	0.619	0.957	0.933	0.786	0.487

Table 4. Table 6: Simulation performance in estimating population simulated mean.

	$\text{True }{\overline{Y}}_g^{\text{pop}}$	0.273	0.453	0.498	0.633	0.725
Bayes	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	0.004	0	0	-0.002	-0.004
	$\text{Mean }\sqrt{{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}){}^2}$	0.013	0.012	0.013	0.013	0.014
	$95$% HPD Width	0.049	0.055	0.055	0.053	0.05
	$95$% HPD Coverage	0.974	0.972	0.957	0.962	0.976
Sample	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	0	0	0	0.001	0
Mean	$\text{Mean }\sqrt{{\left({\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}\right)}^2}$	0	0	0	0	0
	$95$% CI Width	0.047	0.052	0.052	0.047	0.042
	$95$% CI Coverage	0.914	0.936	0.894	0.9	0.899
Bootstrap	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	-0.001	0	0	0.001	0
	$\text{Mean }\sqrt{{\left({\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}\right)}^2}$	0.014	0.015	0.016	0.015	0.012
	$95$% BI Width	0.054	0.061	0.06	0.056	0.051
	$95$% BI Coverage	0.939	0.955	0.936	0.937	0.951

Table 5. Table 7: Bayesian estimator performance measuring other model parameters

	$\overrightarrow{\alpha }$	$\overrightarrow{\mu }$	$𝐏$	$\overrightarrow{S}$	$\overrightarrow{R}$	$𝜼$
$\sqrt{}\text{Average Squared Bias }\widehat{𝜽}$	$6.994$	$0.001$	$0.002$	$31.091$	$0.001$	$0.026$
$\sqrt{}\text{MSE for }\widehat{𝜽}$	$13.383$	$0.022$	$0.046$	$52.129$	$0.019$	$0.591$
Mean HPD Width	$53.007$	$0.06$	$0.157$	$189.364$	$0.039$	$2.482$
Mean HPD Coverage	$0.882$	$0.968$	$0.931$	$0.969$	$0.927$	$0.964$

Table 6. Table 8: Simulation estimating Y ¯ ^ g pop superscript subscript ^ ¯ 𝑌 𝑔 pop \hat{\bar{Y}}_{g}^{\mbox{\tiny{pop}}} from alternate linear mixed-effects

	$\text{True }{\hat{\overline{Y}}}_g^{\text{pop}}$	0.256	0.451	0.499	0.645	0.744
Bayes	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	0.001	0	0	-0.001	-0.001
	$\text{Mean }\sqrt{{\left({\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}\right)}^2}$	0.005	0.006	0.006	0.006	0.005
	$95$% HPD Width	0.025	0.031	0.031	0.029	0.025
	$95$% HPD Coverage	0.99	0.985	0.982	0.986	0.986
Sample	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	0	0	-0.001	0	0
Mean	$\text{Mean }\sqrt{{\left({\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}\right)}^2}$	0.008	0.009	0.01	0.008	0.007
	$95$% CI Width	0.031	0.038	0.038	0.034	0.028
	$95$% CI Coverage	0.947	0.956	0.946	0.951	0.948
Bootstrap	$\text{Mean }{\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}$	0	0	0	0	0
	$\text{Mean }\sqrt{{\left({\hat{\overline{Y}}}_g-{\overline{Y}}_g^{\text{pop}}\right)}^2}$	0.009	0.011	0.011	0.009	0.008
	$95$% BI Width	0.033	0.041	0.042	0.037	0.03
	$95$% BI Coverage	0.918	0.935	0.939	0.945	0.945

Table 7. Table 9: Performance of a ^ k subscript ^ 𝑎 𝑘 \hat{a}_{k} estimates in Bayesian WBC model

	True $a_k$	-0.625	-0.125	0	0.125	0.25	0.375
Bayes	$\text{Mean }{\hat{a}}_k-a_k$	-0.003	0.006	0	-0.001	0	-0.043
	$\text{Mean }\sqrt{{\left({\hat{a}}_k-a_k\right)}^2}$	0.108	0.077	0.075	0.077	0.08	0.09
	$95$% HPD Width	0.624	0.608	0.611	0.609	0.602	0.616
	$95$% HPD Coverage	0.998	0.998	0.999	1	0.999	0.997

Table 8. Table 10: Repeat of Table  9 where m = 0 → 𝑚 0 \overset{{}_{\shortrightarrow}}{m}=0 is assumed

	True $a_k$	-0.625	-0.125	0	0.125	0.25	0.375
Bayes	$\text{Mean }{\hat{a}}_k-a_k$	0	0	0	0	0	0
	$\text{Mean }\sqrt{{\left({\hat{a}}_k-a_k\right)}^2}$	0.044	0.042	0.043	0.041	0.042	0.04
	$95$% HPD Width	0.176	0.172	0.171	0.17	0.173	0.172
	$95$% HPD Coverage	0.965	0.964	0.956	0.961	0.96	0.964

Table 9. Table 11: Simulation where m → 𝑚 \overset{{}_{\shortrightarrow}}{m} is estimated featuring 16 randomly sampled crosses

	True $a_k$	-0.625	-0.125	0	0.125	0.25	0.375
Bayes	$\text{Mean }{\hat{a}}_k-a_k$	-0.002	0.001	0.002	-0.001	0	0.001
	$\text{Mean }\sqrt{{\left({\hat{a}}_k-a_k\right)}^2}$	0.029	0.029	0.029	0.028	0.028	0.028
	$95$% HPD Width	0.124	0.119	0.119	0.118	0.119	0.121
	$95$% HPD Coverage	0.967	0.965	0.954	0.959	0.971	0.968

Table 10. Table 12: Estimate μ ^ g subscript ^ 𝜇 𝑔 \hat{\mu}_{g} , tail probability ( p g subscript 𝑝 𝑔 p_{g} ) against the hypothesis ℋ 0 : μ g = .5 : subscript ℋ 0 subscript 𝜇 𝑔 .5 \mathcal{H}_{0}:\mu_{g}=.5 .

	${\text{<em>Xce</em>}}^a$		${\text{<em>Xce</em>}}^b$		${\text{<em>Xce</em>}}^c$
	${\mu }_g$	$p_g$	${\mu }_g$	$p_g$	${\mu }_g$	$p_g$
ALS	$\begin{array}{c}0.59\\ (0.54,0.63)\end{array}$	$0.003$	$\begin{array}{c}0.37\\ (0.33,0.41)\end{array}$	$\text{**}0.003$	$\begin{array}{c}0.34\\ (0.25,0.43)\end{array}$	$0.002$
SJL	$\begin{array}{c}0.57\\ (0.54,0.61)\end{array}$	$0.003$	$\begin{array}{c}0.48\\ (0.45,0.51)\end{array}$	$\text{*}0.142$	$\begin{array}{c}0.31\\ (0.27,0.35)\end{array}$
LEWES	$\begin{array}{c}0.6\\ (0.56,0.64)\end{array}$	$0.003$	$\begin{array}{c}0.49\\ (0.44,0.53)\end{array}$	$\text{*}0.683$
WLA	$\begin{array}{c}0.55\\ (0.51,0.59)\end{array}$	$0.009$	$\begin{array}{c}0.48\\ (0.44,0.53)\end{array}$	$\text{*}0.45$
WSB	$\begin{array}{c}0.6\\ (0.52,0.67)\end{array}$	$0.022$	$\begin{array}{c}0.47\\ (0.37,0.57)\end{array}$	$\text{*}0.545$	$\begin{array}{c}0.25\\ (0.17,0.33)\end{array}$	$0.006$
PWK	$\begin{array}{c}0.67\\ (0.57,0.76)\end{array}$	$0.008$	$\begin{array}{c}0.61\\ (0.54,0.67)\end{array}$	$0.007$	$\begin{array}{c}0.21\\ (0.15,0.28)\end{array}$	$0.005$
TIRANO	$\begin{array}{c}0.61\\ (0.56,0.66)\end{array}$		$\begin{array}{c}0.51\\ (0.47,0.55)\end{array}$	$\text{*}0.511$
ZALENDE	$\begin{array}{c}0.6\\ (0.56,0.63)\end{array}$		$\begin{array}{c}0.56\\ (0.5,0.61)\end{array}$	$\text{*}0.076$
PERA			$\begin{array}{c}0.54\\ (0.49,0.59)\end{array}$	$\text{*}0.1$
PANCEVO	$\begin{array}{c}0.46\\ (0.41,0.5)\end{array}$	$0.082$	$\begin{array}{c}0.41\\ (0.37,0.44)\end{array}$	$0.003$

Table 11. Table 13: Logit-scale estimated allele effects

Allele	$\begin{array}{c}\text{Additive}\\ \stackrel{}{a}\end{array}$
${\text{<em>Xce</em>}}^a$	-0.36(-0.52, -0.18)
${\text{<em>Xce</em>}}^b$	-0.01(-0.19, 0.19)
${\text{<em>Xce</em>}}^c$	1.08(0.47, 1.67)
${\text{<em>Xce</em>}}^p$	-0.23(-0.42, -0.03)
${\text{<em>Xce</em>}}^v$	-0.67(-0.88, -0.45)

Table 12. Table 14: S   ← ⁡ ←   absent 𝑆 \overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S} precisions from pyrosequencing (brain/kidney/liver) measurements

	Fgd1	Ddx26b	Rragb1	Mid2
brain	$440.6$	$103.4$	$123$	$385.2$
	($214.3$, $692.2$)	($77.1$, $130.6$)	($70.2$, $182.8$)	($138.5$, $683.7$)
kidney	$251.8$	$49.4$	$120.7$	$99.3$
	($169.6$, $349.3$)	($32.9$, $66.3$)	($74.2$, $169.2$)	($0.04$, $296.2$)
liver	$81.2$	$31.4$	$30.8$	$105.9$
	($61.6$, $101$)	($23.1$, $39.9$)	($20$, $42.4$)	($30.2$, $192.4$)
	Zfp182	Xist	Hprt1	Acs14
brain	$289.2$	$25.5$	$29.5$	$26.7$
	($62.2$, $581.5$)	($1.6\times {10}^{-14}$, $135$)	($1.2\times {10}^{-15}$, $151.2$)	($1.3\times {10}^{-26}$, $139.9$)
kidney	$109.6$	$30.4$	$25.5$	$26.5$
	($0.1$, $330$)	($1.2\times {10}^{-11}$, $152.2$)	($4.3\times {10}^{-18}$, $131.6$)	($1.4\times {10}^{-15}$, $137.2$)
liver	$349$	$28.5$	$27.3$	$26.8$
	($84.9$, $650.4$)	($4.3\times {10}^{-13}$, $143.3$)	($1.2\times {10}^{-18}$, $142.6$)	($4.2\times {10}^{-16}$, $138.9$)

Table 13. Table 15: Precision in genes from RNA-seq Measurements (brain only)

Gripap1	Zdhhc9	Ids	Pls3	Arhgef9
$232.4$	$114.7$	$534.2$	$175.4$	$42.1$
($102.5$, $376.2$)	($45.3$, $191.9$)	($236.1$, $873.8$)	($73.1$, $303.4$)	($11.9$, $76$)
Zmym3	Sh3bgr1	Gprasp1	Gnl3l	Tspyl5
$256.7$	$120.1$	$8.7$	$126.9$	$220.1$
($108.2$, $434.4$)	($46.9$, $206.3$)	($1.1$, $17$)	($49.1$, $213.1$)	($85.8$, $389$)

Equations

$$\begin{array}{rl}\text{mbox}Observed:Y_{ij}\sim \text{mbox}Beta(\text{mbox}{P}_i{R}_j& S_j{e}^{{\eta }_{g,j}}+1,\text{mbox}(1-P_i)(1-R_j)S_j{e}^{{\eta }_{g,j}}+1)\text{mbox},\\ \text{mbox}WhereLatent:\underset{\text{mbox}WholeBody}{\underbrace{P_i}}\text{mbox}\sim \text{mbox}\text{mbox}Beta(\text{mbox}& {\mu }_g{\alpha }_g+1,(1-{\mu }_g){\alpha }_g+1)\text{mbox}.\end{array}$$

$${\mu }_1,{\mu }_2,\dots ,{\mu }_G\sim \text{mbox}i.i.d.\text{mbox}Beta({\mu }_{\text{mbox}All}{\alpha }_{\text{mbox}All}+1,(1-{\mu }_{\text{mbox}All}){\alpha }_{\text{mbox}All}+1)\text{mbox}.$$

$$R_j\sim \text{mbox}Marginally\text{mbox}Beta(u_R,u_R)\text{mbox}conditionalon\sum _j\log \frac{R_j}{1-R_j}=0\text{mbox}.$$

$$\sqrt{n_g}({\hat{\mu }}_g-{\mu }_g)⟹N\left(0,{\mu }_g(1-{\mu }_g)({\sigma }_u^2+{\sigma }_e^2/J)\right)+\mathcal{O}\left(1/S_j^2\right)$$

$$\begin{array}{rl}\mathcal{P}(\stackrel{{}_{\text{shortrightarrow}}}{\mu },\stackrel{{}_{\text{shortrightarrow}}}{\alpha },{\mu }_{\text{mbox}All},\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{S},& \stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{R},\mathbf{\eta },{\xi }_S,{\chi }_S,u_R,{\tau }_{\eta }^2,\stackrel{{}_{\text{mbox}\text{shortrightarrow}}}{P}|\mathbf{Y})\propto \text{mbox}Prior(\stackrel{{}_{\text{shortrightarrow}}}{\mu },\stackrel{{}_{\text{shortrightarrow}}}{\alpha },\dots ,u_R,{\tau }_{\eta }^2)\times \\ & \prod _{i=1}^N\text{mbox}Be\text{mbox}ta(P_i;{\mu }_g{\alpha }_g+1,(1-{\mu }_g){\alpha }_g+1)\times \\ & \text{mbox}\prod _{j=1}^J\text{mbox}Beta(Y_{ij};P_i{R}_j{S}_j{e}^{{\eta }_{g,j}}+1,(1-P_i)(1-R_j)S_j{e}^{{\eta }_{g,j}}+1)\text{mbox}.\end{array}$$

$$U_{1}\sim\mbox{Uniform}\left(0,f_{\alpha_{g}}(\alpha_{g}^{(t)}|\overset{{}_{\shortrightarrow}}{\alpha}^{(t)}_{\mbox{\footnotesize{$\backslash g$}}},\overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)})\right)\mbox{.}$$

$${\stackrel{{}_{\text{shortrightarrow}}}{\eta }}_{,j}^{(t)}(g,\Delta )=\left({\eta }_{1,j}^{(t)}-\frac{\Delta }{G-1},{\eta }_{2,j}^{(t)}-\frac{\Delta }{G-1},\dots ,{\eta }_{g,j}^{(t)}+\Delta ,\dots ,{\eta }_{G,j}^{(t)}-\frac{\Delta }{G-1}\right)$$

$$q_{\mathbf{\omega }}(\mathbf{\omega })=\frac{f_{\mathbf{\theta }}({\mathbf{A}}_1\mathbf{\omega }+{\vec{b}}_{1:p^{\prime }},{\mathbf{A}}_2\mathbf{\omega }+{\vec{b}}_{(p^{\prime }+1):p})|{\mathbf{A}}_1|}{\int f_{\mathbf{\theta }}({\mathbf{A}}_1{\mathbf{\omega }}^{\prime }+{\vec{b}}_{1:p^{\prime }},{\mathbf{A}}_2{\mathbf{\omega }}^{\prime }+{\vec{b}}_{(p^{\prime }+1):p})|{\mathbf{A}}_1|d{\omega }_1^{\prime }\dots d{\omega }_p^{\prime }}.$$

$$\left|\begin{array}[]{cccc}\frac{\partial}{\partial L_{1}}\frac{e^{L_{1}}}{1+e^{L_{1}}}&\frac{\partial}{\partial L_{2}}\frac{e^{L_{1}}}{1+e^{L_{1}}}&\ldots&\frac{\partial}{\partial L_{J}}\frac{e^{L_{1}}}{1+e^{L_{1}}}\\ \frac{\partial}{\partial L_{1}}\frac{e^{L_{2}}}{1+e^{L_{2}}}&\frac{\partial}{\partial L_{2}}\frac{e^{L_{2}}}{1+e^{L_{2}}}&\ldots&\ldots\\ \ldots&\ldots&\ldots&\ldots\\ \frac{\partial}{\partial L_{J}}\frac{e^{L_{1}}}{1+e^{L_{1}}}&\ldots&\ldots&\frac{\partial}{\partial L_{J}}\frac{e^{L_{J}}}{1+e^{L_{J}}}\end{array}\right|=\prod_{j=1}^{J}R_{j}(1-R_{j})\mbox{ .}$$

$$f_{\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{L}}(\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{L}|\mathbf{Y},\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{S},\stackrel{{}_{\text{mbox}\text{shortrightarrow}}}{P},\mathbf{\eta },u_R)=f_{\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{R}}(\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{R}(\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{L})|\mathbf{Y},\stackrel{{}_{\text{mbox}\text{shortleftarrow}}}{S},\stackrel{{}_{\text{mbox}\text{shortrightarrow}}}{P},\mathbf{\eta },u_R)\prod _j{R}_j(L_j)\left(1-R_j(L_j)\right).$$

$$\text{mbox}logit({\mu }_g)=(a_{d(g)}+m_{d(g)})-(a_{s(g)}-m_{s(g)})\text{mbox}.$$

$$\text{mbox}logit(Y_{ij})={\mu }_{g(i)}+{\delta }_i+\text{mbox}{Bias}_j+N(0,{\sigma }_j^2)\text{mbox}.$$

$$\begin{split}\mathcal{P}(\bm{\Theta}|\mathbf{Y})\propto&\prod_{i=1}^{n}f_{\text{trans}}^{i}(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i}|Y_{\mbox{\tiny{$<i$}}},\mathbf{I}_{\mbox{\tiny{$<i$}}})\times\\ &\mbox{ }\mbox{ }\mbox{ }\mbox{ }\mbox{ }\mbox{\small{Beta}}(P_{i};\sum_{g}I_{i,g}\mu_{g}\alpha_{g}+1,\sum_{g}I_{i,g}(1-\mu_{g})\alpha_{g}+1)\times\\ &\mbox{ }\mbox{ }\mbox{ }\mbox{ }\mbox{ }\prod_{j}\mbox{\small{Beta}}(Y_{ij};P_{i}R_{j}S_{j}e^{\eta_{g,j}}+1,(1-P_{i})(1-R_{j})S_{j}e^{\eta_{g,j}}+1)\\ &\times\mbox{Priors}(\mathbf{S},\mathbf{R},\bm{\mu},\bm{\alpha},\bm{\eta})\mbox{ . }\end{split}$$

Full text

1 Introduction

X-inactivation bias is a phenomenon observed in the female mouse, first proposed in Cattanach and Isaacson (1967) to be caused by a gene Xce (X-controlling element), which has unknown location and mechanism. In mammals, cells of the early female (XX) embryo inactivate one of their parental chromosomes through a random choice, and pass this decision onto their daughter cells (Lyon, 1962; Gendrel and Heard, 2011), as shown in Figure 1. Thus, in mammals, tissues have mosaic regions of parental chromosome preference; this can be seen, for example, in the fur of a calico cat, whose two fur colors reflect its two parents. Cattanach and Isaacson (1965, 1967) first discovered an X-inactivation bias in mice by shaving mice with parents of white and black hair follicle color, counting the ratio of colors in the hybrid animal, and observing an average bias toward one color. In our related report, Calaway et al. (2013), we proposed to investigate for the location of Xce, refining an interval proposed in Chadwick et al. (2006), by ascertaining the allele carried by 10 unclassified strains. For experimental reasons relating to cost and accuracy, that study measured allele-specific expression in infant hybrids using the technology of pyrosequencing.

Here we provide a model for these data, specific to the needs of the experiment, including the adaptive nature of the experimental design. We propose a Bayesian MCMC-based hierarchical model, based upon the beta-distribution. The model relates allele-specific proportions to a summary proportion reflecting the initial proportion of progenitor cells that conducted the initial choice around 5 days in the embryo development. Our use of the beta distribution suggests a Pòlya-urn model, which can estimate the number of cells in each organ at X-inactivation. We use posterior tail values to argue for the identity of the allele carried by an unknown strain, based on its average bias against known allele-carriers.

To make our model identifiable, we propose priors for sets of our parameters that are constrained to a manifold condition. Through reparameterizing the slice-sample distribution, we can generate acceptable draws from the correct, constrained posterior. Constrained priors become especially necessary when we attempt to design a physical model for the X-inactivation bias observed in the crosses. To predict X-inactivation bias in potential crosses between alleles, we fit a “weight-biased coin” physical model. This physical model supports the hypothesis that Xce is a genetic locus based upon copy-number variation (Stankiewicz and Lupski, 2010; Sebat et al., 2004). Our simulations and analysis here strengthen the case for the findings of the  Calaway et al. (2013) Xce experiment and demonstrate value for Bayesian hierarchical modeling in future allele-specific-expression experiments.

2 Allele-specific Expression Experiment for Locating Xce

Prior to this experiment, four alleles ($a,b,c,e$) of Xce were known in mice, found, for example, in standard strains of mice AJ: $a$, C657B6/J (B6): $b$, CAST/EiJ(CAST): $c$, and in SPRET/EiJ, of the separate mouse species Spretus:$e$ (Cattanach and Isaacson, 1967; Cattanach and Rasberry, 1994, 1991). We have related these alleles with respect to the order of their “strength”. An AJ $\times$B6 hybrid which is $a\times b$ will favor expressing the B6 X-chromosome on average in approximately 65% of their cells; an $a\times c$ cross produces approximately a 30-70 ratio; for $a\times e$ it is 20-80. Even when the population average for a given hybrid is 50-50, there is considerable variation in individual “skew”, i.e. the proportion an individual might over-express one parental X-chromosome. Skewed ratios of 75-25, 27-72, etc. are observed in individuals, and the allele specific expression individual genes might present additional imprecision or biases.

Chadwick et al. (2006) proposed a location of Xce to somewhere within a $\sim$1.9 Mb (million base-pairs) region on the X-chromosome. This was based upon identification of the allele carried by strains of mice, originally conducted with techniques like shaving and counting hair follicles of hybrids with differently colored parents. But counting hair follicles imprecisely measures X-inactivation only in mice old enough to be shaven and only in skin tissues. Given the costs involved in studying embryos and laboratory expertise, Calaway et al. (2013) sought an efficient and reliable measure of aggregate whole-body expression taken from newborn tissue. Due to variation of X-inactivation skew among clones of identical genome, multiple individuals must be measured to establish that a given “first-generation hybrid” ($\text{F}_{1}$) cross shows a bias, favoring measurement techniques of low cost.

Pyrosequencing, a “sequencing by synthesis method” introduced in Ronaghi et al. (1998), was viewed as a cost-effective technology to measure gene expression for strains whose genome sequences were known but whose Xce alleles were not. Calaway et al. (2013) chose commonly expressed X-located genes: Ddx26b, Fgd1, Rragb, Mid2, Zfp185, Xist, Hprt1, and Acsl4, which featured known SNPs that differentiate two strains. Pyrosequencing reads target-sequences of length 20-30 base pairs including either allele of a SNP: Exponential range PCR amplifies up to 1 micro-liter ($\mu l$) of the target, or approximately $3.7\times 10^{10}$ copies, which, after purification, are counted through synthesis. Pyrosequencers report measurements in proprietary fluorescence units, where, conservatively, the final ratio of allele-specific gene expression is based upon more than $1$ million targeted reads.

Table 1 gives a short reference background to the mouse types used in this experiment. There are three mouse sub-species: mus musculus domesticus from Europe, mus musculus castaneus from Africa, India and the Arabian subcontinent, and mus musculus musculus from Asia, that rarely interbreed in the wild (PANCEVO is a different species: mus spicilegus, “Steppe Mouse”). Laboratory-derived strains were initially generated from domesticated toy mice. Wild-derived strains come from wild-captured samples, inbred to be representatives of a wild subpopulation. We wish to classify Xce allele membership for strains marked with unknown allele “??”, and thus narrow down a candidate region for Xce.

2.1 Implied Sequential Design of the Experiment

The Calaway et al. (2013) first established which Xce allele was carried by a few unclassified strains. Based upon the observations, a provisional, semi-statistical classification was assumed, and this motivated the choice to classify additional strains. First, strains SJL/J (SJL) and ALS/LtJ (ALS) were analyzed because their sequences split the known region between $\mbox{\emph{Xce}}^{a}$ and $\mbox{\emph{Xce}}^{b}$. Then strains LEWES/EiJ (LEWES) and WLA/Pas (WLA) followed. Additional crosses with WSB/EiJ (WSB), PWK/PhJ (PWK) were generated and measured, including some RNA-seq measurements. Available frozen $\text{F}_{2}$ samples of crosses including DDK/Pas (DDK), ZALENDE/Ei (ZALENDE), PERA/EiJ (PERA) (Nishioka, 1995) were acquired. Some successful breedings with PANCEVO/EiJ (PANCEVO) (Kim et al., 2005) were conducted to see if this second species of mouse had a different Xce profile. Only a subset of possible crosses were performed, and not all reciprocal crosses were generated. Count of replicates per unique crosses varied from 4 to 50.

For convenience, summarizing details for the inbred parental strains are given in Table 1. A statistical analysis must respect that the experiment is adaptive and sequential in nature, that many crosses which might have aided in the statistical estimate will be missing, and that we should provide a statistical criterion not only to assert that a given cross SJL $\times\texttt{B6}$ has an average “$50.1\%\pm 2.1\%$”, but that we should also be able to express how certain we are that SJL is an $a$ allele carrier. By choosing closely related potential $\mbox{\emph{Xce}}^{a}$ and $\mbox{\emph{Xce}}^{b}$ allele carriers, predominantly from mus domesticus,  Calaway et al. (2013) proceeded with an investigation that, aided by this statistical analysis, refined the Xce location from the 1.9 Mega-base (Mb) Chadwick et al. (2006) candidate interval to a 176 kilo-base (kb) region approximately 500 kb from the gene Xist whose expression causes X-inactivation. Calaway et al. (2013) established that allele $\mbox{\emph{Xce}}^{a}$ is rare in the wild and that most individuals within a sub-species share the same Xce allele. They found that Xce-allele strengths can potentially increase or decrease during speciation, which suggests Xce is mutated through copy number variation.

2.2 Format of Pyrosequencing Data

A partial Table 2 shows the challenges of the Calaway et al. (2013) data. Every mouse is a member of a single cross, which is coded to correspond to a unique hybrid combination of mother and father inbred strains. We encode every cross with a short code, such as “1Wl” for a 129 $\times$WLA cross, where 129 is the maternal strain. Note that a WLA $\times$129 hybrid, with WLA as maternal strain, would be coded “Wl1” and coded as a different reciprocal cross. The other columns are measurements of maternal proportion expressed for given tissue-gene combinations. Allele-specific proportion is measured in multiple tissue-genes, with considerable missingness: “NA” for genes that were either dropped at random for an individual mouse during experiment, or because they are not experimentally viable for the cross when no SNP exists on that gene to differentiate maternal and paternal copies. Since SNPs are single base-pair substitutions, between any three parental strains there is no single SNP that can be used to measure allele-specific expression for all three pairings of two strains.

If we naively take sample means of all observed gene expression in individuals stemming from one type of cross, we can summarize those averages in a table like Table 3, where we see that of candidates with historically established alleles $a$, $b$, $c$, there are similar behaviors encountered in crosses with previously unknown PWK, SJL, WSB strains. Reciprocal crosses, such as the $b\times a$ cross, have often not been measured. We inevitably hope for a statistical method that can establish the allele carried by unestablished strains.

Allele-specific expression could have been measured through Illumina RNA-seq (Morin et al., 2008) or microarrays (Chang, 1983; Hall et al., 2007). RNA-seq whole-transcriptome sequencing enables simultaneous measurement of allele-specific expression of all 473 X-located genes. But typical RNA-seq will only have on order 500 read counts of any SNP region, offering less precision than pyrosequencing. A Bayesian method for estimating precision and bias with RNA-seq count data can be found in Graze et al. (2012). Both RNA-seq and pyrosequencing are unbiased at preserving true, initial allelic stoichiometry. But the non-linear relationship between expression and luminescence in microarrays can lead to biased estimation of allele-specific proportion.

Late in the Calaway et al. (2013) experiment, brain-tissue RNA-seq measurements were made available for 34 individuals from $\text{F}_{1}$ crosses between wild-derived sub-species representatives WSB, PWK, and CAST, about $6$ individuals of each cross, including reciprocal crosses. It was hoped that our statistical procedure designed for pyrosequencing output could accommodate RNA-seq data. As we show, RNA-seq measurements serve the same role in our model as pyrosequencing, although RNA-seq measurements often have less precision when baseline expression of a gene is low.

3 Statistical Model for Gene-expression Ratios

Denote the proportion of maternal gene-expression measurements $Y_{ij}\in[0.001,.999]$, where $i\in\{1,\ldots,n\}$ is the index of the mouse specimen used in the experiment, and $j\in\left\{1,\ldots,J\right\}$ is an index of tissue-gene combination (examples: “kidney-Fgd1”, “brain-Ddx26b”). We are given $Y_{ij}$ as a fractional proportion of maternal expression over total expression. Certain genes, dependent on tissues, could be biased or imprecise measures of true overall X-inactivation. For instance, the gene Fgd1, which differentiates AJ and CAST, might be biased in favor of AJ expression within brain tissue but not in the liver. Measuring multiple tissues in multiple genes, we hope to estimate a quantity, $P_{i}$, representing a best whole-body proportion of active maternal X-chromosome.

Even if mouse $i$ has a whole-body proportion of $P_{i}$, replicate mice from the same cross $i^{\prime},i^{\prime\prime},\ldots$ will have different proportions $P_{i^{\prime}},P_{i^{\prime\prime}},\ldots$. The cross-specific mean $\mu_{g}$ denote the population mean for mice in cross/group $g\in\{1,\ldots,G\}$. We seek then a hierarchical method that estimates $P_{i}$ from observed $Y_{ij}$, and relates $P_{i}$ for individuals $i$ in group $g$ to their mean $\mu_{g}$. Only a subset of genes $j$ can be observed for certain crosses $g$. Some data is missing-at-random, but much $Y_{ij}$ is systematically unmeasurable, due to shared SNPs for that gene. Our inference model follows.

3.1 A Hierarchical Beta Model

Our scientific objectives are served by a Bayesian hierarchical model. Observed data are proportions $Y_{ij}$, so a beta-distribution regression model (Ferrari and Cribari-Neto, 2004) can model unimodal proportions. For our problem, the mean and mode for the beta densities should be constrained within $(0,1)$, excluding [math] and $1$. Thus, we choose parameterized beta distributions of the form $\mbox{Beta}(\mathbb{A}+1,\mathbb{B}+1)$ where $\mathbb{A},\mathbb{B}\geq 0$. We first give the parameterized mathematical model, and then describe the purpose of the parameters in the model. Let:

[TABLE]

In the above equation, we have latent, unobserved ratios $P_{i}$ for $i\in\{1,\ldots,n\}$ for each individual and parameter vectors $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}\in\mathbb{R}^{J}$, $\overset{{}_{\shortrightarrow}}{\alpha},\overset{{}_{\shortrightarrow}}{\mu}\in\mathbb{R}^{G}$ and a perturbation matrix $\bm{\eta}\in\mathbb{R}^{G\times J}$. In our notation, $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ and $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ are row-vectors each of length $J$, the number of columns of Table 2, while $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P},\overset{{}_{\shortrightarrow}}{\mu},\overset{{}_{\shortrightarrow}}{\alpha}$ are column-vectors for the $n$ individuals or $G$ crosses. The matrix $\bm{\eta}$ is of dimension $G\times J$.

$R_{j}\in(0,1)$ represents bias for tissue-gene combination $j\in\{1,\ldots,J\}$. The value $R_{j}=.5$ suggests that this tissue-gene is relatively unbiased as a predictor of overall X-inactivation. If a tissue-gene has $R_{j}=.6$, then we expect this to push $Y_{ij}$, on average for all crosses, in the direction of maternal expression.

In contrast, $S_{j}\in(0,\infty)$, represent precision of measurements for a given tissue-gene combination $j$. Larger values of $S_{j}$ imply that tissue-gene $j$ is a more precise measurement of $P_{i}$. To add further difference in precision for certain crosses, parameters $\eta_{g,j}\in\{-\infty,\infty\}$ model deviation in precision for crosses $g$ oftissue-gene measurements $j$ as compared to other crosses. The precision for cross $g$ on tissue-gene $j$ will be $S_{j}\times e^{\eta_{g,j}}$. Because the Beta-distribution implies a Pòlya-urn model for cell proliferation, we can use measurement of $S_{j}$ to infer in some measure how many progenitor cells exist in each organ during the X-inactivation event. Let priors for $S_{j}$ be $\mbox{Gamma}(\chi_{S},\xi_{S})\propto e^{-x/\xi_{S}}x^{\chi_{S}-1}$ with prior for hyper-parameters $(\xi_{S},\chi_{S})\sim\mbox{Gamma}(.1,.1)$.

The $+1$ to both terms of $\mbox{Beta}(\mathbb{A}+1,\mathbb{B}+1)$ ensure stability of the posterior, keeping the beta density $\propto x^{\mathbb{A}+1-1}(1-x)^{\mathbb{B}+1-1}$ from having infinite boundary density if $P_{i},\mu_{g}\approx 1$ or if $P_{i},\mu_{g}\approx 0$. Since we restrict our population of genes and individuals to those expressing between $1\%$ and $99\%$ of the maternal X (individuals expressing more or less are presumed either males or XO females carrying only a single X-chromosome due to a XY separation error in the gamete), the $(+1,+1)$ information has negligible effect on mid-range values and reflects our sampling choices.

Hyperparameters $\mu_{g},\alpha_{g}$ represent the average and variance for $P_{i}$ within a cross group $g$. $\mu_{g}\in(0,1)$ reflects the average $P_{i}$ proportion of cross $g$, and is the most important parameter for establishing Xce membership. Nuisance parameter $\alpha_{g}\in(0,\infty)$ represents the variability of $P_{i}$ about $\mu_{g}$. For large $\alpha_{g}$ values, $P_{i}$ variation about $\mu_{g}$ will be less. We assign $\overset{{}_{\shortrightarrow}}{\mu}$ a hierarchical prior:

[TABLE]

For Equation 2, the global mean and precision parameters, $\mu_{\mbox{\tiny{All}}},\alpha_{\mbox{\tiny{All}}}$ are assigned weak hyper-priors, $\mbox{Unif}(0,1)$ and $\mbox{Gamma}(.1,.1)$ respectively. For vector $\overset{{}_{\shortrightarrow}}{\alpha}$, which represents the spread of observed $P_{i}$ around their group mean $\mu_{g}$, we assign an i.i.d. $\mbox{Gamma}(1,1)$ exponential prior.

3.2 Constrained Priors

Due to identifiability, canonical i.i.d. priors do not suit our purposes for parameters $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ and $\bm{\eta}$. If all $J$ bias parameters $R_{j}$ deviated from $.5$ together, such as if \bar{R}$$= $\frac{1}{J}\sum_{j}R_{j}=.6$, the model would shift $\hat{P}_{i}$ to the $\mbox{sign}(.5-\bar{R})$ direction of all $Y_{ij}$. But a global movement of $\hat{P}_{i}$ would shift all $\hat{\mu}_{g}$. Relative differences $\hat{\mu}_{g_{1}}-\hat{\mu}_{g_{2}}$ might be preserved, but we would not obtain a precise estimate of $\hat{\mu}_{g}$, if confounded with $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$.

We explicitly assume (because we have no way of discovering otherwise) that tissue-gene combinations are “on-average” unbiased: $\sum_{j}\mbox{logit}(R_{j})=0$, where $\mbox{logit}(R_{j})=\log_{e}\frac{R_{j}}{1-R_{j}}$. We will report $\hat{\mu}_{g}$ given an assumption that $\mbox{logit}(R_{j})$ is on average unbiased. Since most decisions will be based upon the significance of differences $\hat{\mu}_{g^{\prime}}-\hat{\mu}_{g}$ we feel that this assumption is reasonable. We therefore model bias with a restricted prior of the form:

[TABLE]

which we show in Section 3.3 is suitable for MCMC posterior exploration.

Model parameters $\eta_{g,j}$ represent deviation about parameter $S_{j}$. Since $\eta_{g,j}\in\{-\infty,\infty\}$, Gaussian priors $N(0,\tau^{2}_{\eta})$ might seem acceptable. Certainly, if all $\eta_{g,j}$ deviated from [math] together, then the combination product $e^{\eta_{g,j}}S_{j}$ could support many values of $S_{j}$ to retain the same overall precision. Inferentially, this weakens our estimation ability of $S_{j}$. Instead, we define $\sum_{g}\eta_{g,j}=0$ for every $j$. Thus $S_{j}$ represents an average precision over all groups $g$. The prior for $\eta_{g,j}$ will be $\sum_{g}\eta_{g,j}=0$, but that $\eta_{g,j}\sim N(0,\tau^{2}_{\eta})$ marginally.

Constrained priors limit the parameter space to a lower-dimensional manifold. In Gaussian Gibbs-regression settings, such restrictions are commonly imposed (Gelfand et al., 1992). Vines et al. (1996) willingly accept a 10x slowdown in iteration speed to estimate regressions with constrained parameters for a mixed effects model, and observed reduced autocorrelation in the samples. In non-Gaussian settings,  Zhu et al. (2011); Shi et al. (2009); Zho et al. (2010) show advances in manifold exploration and suggest new ways to explore these parameter spaces. We demonstrate a slice-sampler suitable for our linear constraints in Section 3.3.

Due to orthogonality in design between $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ and $\mu_{g}$, we can approximate that for $J$ fixed and $n_{g}\rightarrow\infty$:

[TABLE]

where $\sigma^{2}_{u}\approx\frac{\alpha_{g}+2}{(\alpha_{g}\mu_{g}+1)(\alpha_{g}(1-\mu_{g})+1)}$ and $\sigma^{2}_{e}$ behaves as $2/S_{j}$, are logit-scale linear mixed regression model group and noise variances respectively. With unconstrained priors, leaving $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ confounding, it is uncertain that $\hat{\mu}_{g}$ would converge.

Table 4 summarizes our parameters and prior choices. We next describe choices in MCMC algorithm that enables investigation of this model, including these constrained priors.

3.3 Parameter Estimation through Slice Sampling

If $\bm{\Theta}$ represents the unknown model parameters, the product of our model-based likelihood, $P(\mathbf{Y}|\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P},\bm{\Theta})\times P(\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}|\bm{\Theta})$ and the priors, $p(\bm{\Theta})$, gives an unnormalized function proportional to the desired posterior $\mathcal{P}(\bm{\Theta}|\mathbf{Y})$:

[TABLE]

We use the technique of slice-sampling (Roberts et al., 2012; Neal, 2003) to achieve MCMC random draws from the posterior distribution $\mathcal{P}(\overset{{}_{\shortrightarrow}}{\mu},\overset{{}_{\shortrightarrow}}{\alpha},\ldots|\mathbf{Y})$. Let group $g\in\{1,\ldots,G\}$, and let $\backslash g$ be $\mathcal{G}-\{g\}$, or the set $\{1,2,\ldots,g-1,g+1,\ldots,G\}$. In slice-sampling, a new parameter $\alpha_{j}^{(t+1)}$ is drawn given current $\alpha_{g}^{(t)}$ as well as fixed \overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\shortrightarrow}}{\alpha}^{(t)}_{\mbox{\footnotesize{\backslash g}}},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)}, by considering the marginal, one dimensional function: f_{\alpha_{g}}(\alpha_{g}|\overset{{}_{\shortrightarrow}}{\alpha}^{(t)}_{\mbox{\footnotesize{\backslash g}}},\overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)}), which is the marginal unnormalized posterior at a given point $\alpha_{g}^{(t+1)}$ with other parameters fixed. In this case, given $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}$ and $\overset{{}_{\shortrightarrow}}{\mu}$, $\overset{{}_{\shortrightarrow}}{\alpha}$ is independent of the $\mu_{\mbox{\tiny{All}}}$, $\alpha_{\mbox{\tiny{All}}}$, $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$, $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$, $\bm{\eta}$, $\xi_{S}$, $\chi_{S}$, $u_{R}$, $\tau^{2}_{\eta}$ parameters.

To perform slice-sampling, first a uniform value $U_{1}$ is drawn from

[TABLE]

Then the slice sampler seeks left and right from point $\alpha_{g}^{(t)}$ to find the points $\alpha^{-}$ below and $\alpha^{+}$ above s.t. f_{\alpha_{g}}(\alpha^{-}|\overset{{}_{\shortrightarrow}}{\alpha}^{(t)}_{\mbox{\footnotesize{\backslash g}}},\overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)}) $=$ f_{\alpha_{g}}(\alpha^{+}|\overset{{}_{\shortrightarrow}}{\alpha}^{(t)}_{\mbox{\footnotesize{\backslash g}}},\overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)}) $=$ $U_{1}$. To complete the draw, a second uniform is drawn from the interval $[\alpha_{g}^{-},\alpha_{g}^{+}]$. This draw serves as the new draw of $\alpha_{g}^{(t+1)}$. Looping this procedure for all parameters serves as the MCMC scheme. The algorithmic order can be understood to reflect the order of unobserved parameters in Table 4.

This procedure is well-known and has an advantage over Metropolis-Hastings procedures in not requiring tuning or a proposal density. Slice-sampling is the backbone to the all-purpose Bayesian integration software “JAGS” (Plummer, 2003; Murphy, 2007). For our purposes, we reimplemented this process explicitly in R SHLIB-compiled C (Developers et al., ) for two reasons. First, we sought to gain efficiency in coding $f(\cdot|\cdot)$ densities by separating out conditionally independent parameters: for instance, conditional on $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}$, the posteriors for $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ and $\overset{{}_{\shortrightarrow}}{\mu}$ are independent. Efficiencies could be gained by using BLAS (Lawson et al., 1979) vector functions: conditional on $\mathbf{Y}$, the posterior for $P_{i}$ is proportional to the vector dot products $\sum_{j}S_{j}e^{\eta_{g(i),j}}R_{j}\log_{e}Y_{ij}$ and $\sum_{j}S_{j}e^{\eta_{g(i),j}}(1-R_{j})\log_{e}(1-Y_{ij})$, which are tabulated faster with block memory management.

But, more crucially, we needed to modify the slice-sampler method to allow for multivariate constrained prior distributions on $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ and $\bm{\eta}$. Recalling that we require $\sum_{g}\eta_{g,j}=0$, we draw new samples for $\eta_{g,j}^{(t+1)}$ not by slice-sampling along the univariate density $f_{\eta_{g,j}}(\eta_{g,j}|\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{(t)},\bm{\eta}_{,/j}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}^{(t)},\mathbf{Y})$ but along the full multivariate posterior density as defined along a vector, $\overset{{}_{\shortrightarrow}}{\eta}_{,j}^{(t)}(g,\Delta)$, which is perturbed by a $\frac{\Delta}{G-1}(Ge_{g}-\vec{1})$, where $e_{g}$ is the standard $g$th basis vector in $\mathbb{R}^{G}$ and $\vec{1}$ is a G-length vector of 1’s. In other words:

[TABLE]

In this case our notation considers perturbation of a column vector $\overset{{}_{\shortrightarrow}}{\eta}_{,j}\in\mathbb{R}^{G}$ which represents all $\eta_{g,j}$ values for a single tissue-gene $j$ fixed and the crosses $g\in\{1,2,\ldots,G\}$.

To perturb $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$, first define a logit basis reparameterization $L_{j}\equiv\log_{e}\frac{R_{j}}{1-R_{j}}$, and take a perturbation, $\frac{\Delta}{J-1}(Je_{j}-\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{1})$, in the $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}$ logit space.

It is important to note a crucial difference in the densities used in slice-sampling of $\bm{\eta}$ versus sampling of $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L})$. One cannot slice-sample over arbitrary reparameterizations without changing measure. For instance, in a univariate case, if $e^{\omega}=\theta$, then, if $q_{\theta}(\theta)$ was the target density then $q_{\omega}(\omega)=q_{\theta}(e^{\omega})\times e^{\omega}$. It would be inappropriate to find the two $\omega^{+},\omega^{-}$ such that $f_{\theta}(e^{\omega^{+}})=f_{\theta}(e^{\omega^{-}})=\mbox{Unif}(0,f_{\theta}(e^{\omega^{(t-1)}}))$ and then sample $\omega^{(t)}=\mbox{Unif}(\omega^{-},\omega^{+})$. One would see a dramatically downward-biased distribution for $\theta(\omega)$. For transformations $\theta=\omega-5$ or $\theta=4\omega$, slice-sampling from $f_{\theta}(\theta(\omega))$ or $f_{\theta}(\theta)$ produces acceptably distributed samples because the linear transformation does not change the relative measure along $d\theta$.

In the multivariate setting, consider parameterization $\omega\in\mathbb{R}^{p^{\prime}}$ that transforms linearly onto the desired parameter space $\theta\in\mathbb{R}^{p}$ with $p^{\prime}\leq p$, such as $\theta=\mathbf{A}\omega+\vec{b}$. Let $\mathbf{A}_{1}$ be the top $p^{\prime}\times p^{\prime}$ square matrix of $\mathbf{A}$ and $\mathbf{A}_{2}$ be the remainder $(p-p^{\prime})\times p^{\prime}$ rectangular matrix of $\mathbf{A}$. If matrix $\mathbf{A}_{1}$ is invertible, then $\bm{\theta}_{1:p^{\prime}}$ or the first $p^{\prime}$ coefficients of $\bm{\theta}$ determine $\bm{\theta}_{(p^{\prime}+1):p}$ which are the last $p-p^{\prime}$ coefficients of $\bm{\theta}$.

If $q_{\bm{\theta}}(\bm{\theta})$ were a target density of $f_{\bm{\theta}}(\bm{\theta})$ of unknown integration constant, subject to constraints, then $q_{\bm{\theta}}(\bm{\theta})=\frac{f_{\bm{\theta}}(\bm{\theta}_{1:p^{\prime}},\bm{\theta}_{(p^{\prime}+1):p})}{\int f_{\bm{\theta}}(\bm{\theta}_{1:p^{\prime}}^{\prime},\bm{\theta}_{(p^{\prime}+1):p}^{\prime})d\theta_{1}^{\prime}\ldots d\theta_{p^{\prime}}^{\prime}}$. But also,

[TABLE]

In Equation 8, the Jacobian $|\mathbf{A}_{1}|$ is part of the integration constant. Slice-sampling will automatically evaluate this constant.

In the case of $\overset{{}_{\shortrightarrow}}{\eta}_{,j}$, the rotation $\mathbf{A}_{1}$ is the matrix of $-1/(G-1)\vec{1}\vec{1}^{T}+\frac{G}{G-1}\mathbf{I}$ in $\mathbb{R}^{G-1}$. If $\mathbf{A}_{2}=1/(G-1)\vec{1}$ and $\vec{b}=\vec{0}$ then the reparameterized space $\overset{{}_{\shortrightarrow}}{\eta}_{,j}^{\prime}\in\mathbb{R}^{G-1}$ generates $\overset{{}_{\shortrightarrow}}{\eta}_{,j}=\{\mathbf{A}_{1}\overset{{}_{\shortrightarrow}}{\eta}_{,j}^{\prime},\mathbf{A}_{2}\overset{{}_{\shortrightarrow}}{\eta}_{,j}^{\prime}\}$ and we need not change measure.

In the case of $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$, which includes a nonlinear transformation $\vec{L}\rightarrow\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$, the Jacobian matrix is of more consequence;

[TABLE]

The change of variables from $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ to $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}$ measure is

[TABLE]

As with $\overset{{}_{\shortrightarrow}}{\eta}_{,j}$, a linear transformation from $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}^{\prime}\in\mathbb{R}^{J-1}$, such that $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}_{1:(J-1)}=\mathbf{A}_{1}\vec{L}^{\prime}$ and $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}_{(p^{\prime}+1):p}=\mathbf{A}_{2}\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{L}^{\prime}$, maintains the $\sum_{j}L_{j}=0$ constraint.

We note that not all $Y_{ij}$ are observed, sometimes because SNPs do not differentiate between crosses, and in other cases because of random experimental missingness. Imputation of the data missing at random is naturally accommodated in a Bayesian setting. After every iteration, the MCMC sampler draws new samples $Y_{ij}^{\mbox{\tiny{miss}}}$ missing data based upon current estimates of $P_{i}$ and $S_{j},R_{j},\eta_{g,j}$. Using these new draws of $Y^{\mbox{\tiny{miss}}}$ to make the data complete, samples of $P_{i},S_{j},R_{j},\eta_{g,j}|\mathbf{Y}^{\mbox{\tiny{complete}}}$ are drawn. When a gene is entirely excluded for a certain cross because of identical alleles, we do not impute $Y_{ij}^{\mbox{\tiny{miss}}}$ for that cross, and treat $\eta_{g,j}$ as essentially non-existent.

Because we vectorize and separate certain operations, computational costs are proportional to $\mathcal{O}(NJ)$ to update all $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ parameters, but $\mathcal{O}(NJ^{2})$ to update all $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ parameters and $\mathcal{O}(JGN)$ to update all $\bm{\eta}$ parameters, as well as $\mathcal{O}(G)$ to update $\overset{{}_{\shortrightarrow}}{\mu},\overset{{}_{\shortrightarrow}}{\alpha}$ parameters. Thus estimation of $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$ seems to be the most inefficient step in the computation, where most of the computation costs come in calculation of beta values $\frac{\Gamma(R_{j}P_{i}S_{j}e^{\eta_{g,j}}+1)\Gamma((1-R_{j})(1-P_{i})S_{j}e^{\eta_{g,j}}+1)}{\Gamma((1-R_{j}-P_{j}+2R_{j}P_{j})S_{j}e^{\eta_{g,j}}+1)}$, which do not allow for the efficient separation of parameters.

Gibbs samples can be used to infer posterior mean, standard-deviation, and quantiles which we use as evidence to simultaneously test biologically relevant hypotheses, such as $\mathcal{H}_{\mu}\mbox{: }\mu_{g}=.5$, $\mathcal{H}_{S}\mbox{: }S_{j_{1}}=S_{j_{2}}$, and $\mathcal{H}_{\eta}\mbox{: }\sum_{g}(\eta_{g,j}-\bar{\eta}_{j})^{2}<10$. In particular, tests for hypotheses of the form $\mathcal{H}_{\mu_{1},\mu_{2}}\mbox{: }\mu_{g_{1}}=\mu_{g_{2}}$ are desired, and we hope to use posterior values to both defend and reject hypotheses of this type.

We can achieve roughly 50K samples in a day, and the Gelman-Rubin (Gelman and Rubin, 1992) convergence measures reported in Section 5 using the complete data are quite low, so we conclude that the chain is optimized sufficiently for the purpose of testing our biological hypotheses.

3.4 Allele Effect Modeling

Previous literature (Ruvinsky, 2009; Kambere and Lane, 2009; Wang et al., 2010) has suggested the alleles in Xce behave in a manner akin to a “weight-biased coin”. Consider a cross between mother carrying Xce allele “$a$” and the father with allele “$b$”. One could imagine putting a coin together with a mass $M_{a}$ for the heads face, and a mass $M_{b}$ for the tails face. If a physical model suggested \mu(a\mbox{ \times }b)=M_{a}/(M_{a}+M_{b}), we could measure $M_{c}$ in a third allele and predict \mu(a\mbox{ \times }c)=M_{a}/(M_{a}+M_{c}). This model is consistent with a hypothesis that Xce varies through copy number variation (CNV). Longer sequences with more copies attract nucleosomes and transcription factors to attach to the X-chromosome and deactivate it, suggesting that longer sequences at Xce are weaker alleles.

Now that we have declared multiple inbred strains of mice to carry $\mbox{\emph{Xce}}^{a}$, $\mbox{\emph{Xce}}^{b}$, $\mbox{\emph{Xce}}^{c}$, …alleles, we can implement this model within our framework. Denote the previous specification, where $\mu_{g}$ are given i.i.d. priors, our “Independence Assumption” (IA) model. After analysis, we can compare IA estimates, where no structure was assumed for $\mu_{g}$, with estimates from this weight-biased coin (WBC) model.

Consider $d(g)$ be the allele carried by the mother or “dam” for cross $g$, and $s(g)$ be the allele carried by the father or “sire” inbred-strain for cross $g$. We model $\mu_{g}$ as a function of additive allele parameters $a_{1},a_{2},\ldots,a_{K_{a}}$ for $K_{a}$ total alleles and parent-of-origin allele-specific parameters $m_{1},m_{2},\ldots,m_{K_{m}}$ for a total number of parent-of-origin alleles $K_{m}$. In our case we propose:

[TABLE]

Having a larger $a_{d(g)}$ dam allele than the $a_{s(g)}$ sire allele pulls $\mu_{g}$ in the direction $\mu_{g}>.5$. Parent-of-origin effects $m_{d(g)}$ perturb this result based upon whether strains serve as dam or sire.

Again, we impose constraints $\sum_{k=1}^{K_{a}}a_{k}=0$ and $\sum_{k=1}^{K_{m}}m_{k}=0$. Vector-constrained slice sampling then allows for posterior sampling of these model parameters, along with the previously mentioned $\bm{\eta},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}$.

4 Verifying Method Performance through Simulation

Gibbs sampler models are computationally slow, and one can rarely explore all regimes to test a Bayesian estimator. It would be wrong to draw parameters from prior distributions to simulate datasets; to demonstrate frequency performance of posterior estimates we must use fixed parameter values. We simulate datasets of size $n=200$ individuals, $J=6$ tissue gene measurements, making $G=5$ groups with $40$ samples each. We fix $\overset{{}_{\shortrightarrow}}{\mu}^{\text{truth}}=(.25,.45,.5,.65,.75)^{T}$, which is an ordered set of realistic values for the group means, where adjacent $\mu_{g}$ values are separated by .2, .05, .15, and .1. We choose tissue-gene precisions $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{\text{truth}}=(200,200,100,100,50,50)$ and set $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}^{\text{truth}}$ such that $\mbox{logit}(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}^{\text{truth}})=(-3,-2,-1,0,1,5)/8$. Note that the most precise gene-expression values will be downward biased. We will keep the true secondary precision matrix $\bm{\eta}^{\text{truth}}=\bm{0}$, however, the Bayesian estimator will still fit $\bm{\eta}$ as a free $5\times 6$ matrix. We set $\overset{{}_{\shortrightarrow}}{\alpha}^{\text{truth}}=(50,50,50,50,50)^{T}$. We draw a $P_{i}$ for every individual $i$ based upon group membership $g(i)$, from which $Y_{ij}$ values are sampled, before 66% of $Y_{ij}$ are removed to represent data missingness (requiring at least one measurement per individual).

Our estimation draws 2000 Gibbs samples from the Bayesian Posterior, with fit time approximately 10 minutes, and disk usage approximately 100mb. We estimate posterior mean and credibility intervals for all model parameters. We ran simulations on the UNC KillDevil Cluster, which allowed roughly 200 cores in parallel, fitting 1000 replications of all of our simulation experiments (including those described later in this section) in roughly a day. We compare Bayesian estimates for $\mu_{g}$ to the sample mean of groups $g$: \bar{Y}_{g}=\frac{1}{\|Y^{\mbox{\tiny{obs}}}_{ij}\|_{0}}\sum_{\mbox{\tiny{g(i)=g}},\mbox{\tiny{j observed}}}Y_{ij}^{\mbox{\tiny{obs}}}, where, confidence intervals are based upon the $t$-distribution using standard deviation estimated from the data.

Simulation results in Table 5, show that sample means $\bar{Y}_{g}$ do not estimate the vector $\mu$ well, but that Bayesian credibility intervals have near-95% coverage of true $\mu_{g}$ values. For this $n=200$ dataset, with these values of $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ and $\vec{\overset{{}_{\shortrightarrow}}{\alpha}}$, Bayesian 95% credibility intervals will be of approximate width 6%. This credibility width is sufficient to differentiate almost all of the group $\mu_{g}^{(t)}$. Only estimates for $\mu_{2}=.45$ and $\mu_{3}=.5$ had overlapping posteriors. We use a rejection criterion $2\times\mbox{min}(\mathcal{P}(\mu_{3}>\mu_{2},\mu_{2}>\mu_{3}))<.05$. That is, we identify two $\mu_{j}$ parameters as different if the two-sided tail posterior probability that their order is reversed is less than $.05$. For the comparison of $\mu_{2},\mu_{3}$, this test has a power of 60.4%. Comparisons of group means for $\mu_{g}-\mu_{g^{\prime}}\geq.1$ have near 100% power.

Arguably, the sample mean, $\bar{Y}_{g}$, is not meant to estimate $\mu_{g}$, which is a quantity of the model. For a fairer comparison, consider the target of a population gene-expression mean $\bar{Y}^{\mbox{\tiny{pop}}}_{g}$ of group $g$. This would be the population average gene expression of the six genes for all possible mice of cross $g$. If we simulate 100K individuals or more from the model, we can calculate a satisfactory approximation \bar{Y}^{\mbox{\tiny{pop}}}_{g}=\frac{1}{J\times 100K}\sum_{\mbox{\tiny{100K simulated pups i^{\prime}}}}Y_{i^{\prime}j}. The sample mean $\bar{Y}_{g}$ of 40 pups should be an estimate for $\bar{Y}_{g}^{\mbox{\tiny{pop}}}$. In Table 6, we see that, in estimation of observable $\bar{Y}_{g}^{\mbox{\tiny{pop}}}$ that the Bayesian $95$% credibility intervals tend to be conservative and over-cover $\bar{Y}_{g}^{\mbox{\tiny{pop}}}$. The sample mean’s nominal $95$% confidence intervals offer only $90$% coverage for $\mu_{3},\mu_{4},\mu_{5}$. The bootstrap has wider bootstrap intervals than the Bayesian credibility intervals, but have less coverage. The Bayesian credibility intervals are .01 narrower, but have over-conservative $97$% coverage of $\bar{Y}_{g}^{\mbox{\tiny{pop}}}$.

From the Bayesian model, we can also generate an estimate for $\bar{Y}_{g}^{\mbox{\tiny{pop}}}$ by simulating data from the posterior predictive distribution for future $Y_{i^{\prime}j}$. For each step $(t)$ in the Gibbs Sampler, we can simulate 100K mice from the model using the current state of $(\overset{{}_{\shortrightarrow}}{\mu}^{(t)},\overset{{}_{\shortrightarrow}}{\alpha}^{(t)})$ and calculate expected gene expressions using $(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{(t)},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}^{(t)},\bm{\eta}^{(t)})$, taking the average \bar{Y}_{g}^{\mbox{\tiny{sim (t)}}}. These samples, \bar{Y}_{g}^{\mbox{\tiny{sim (t)}}}, treated as draws from the posterior distribution for $\bar{Y}^{\mbox{\tiny{pop}}}_{g}$.

We compare estimates for $\bar{Y}^{\mbox{\tiny{pop}}}$ using the Bayesian posterior predictive distribution, the sample mean, as well as the bootstrap mean including 95% bootstrap intervals, in Table 6.

Table 7 presents average performance of the Bayesian model in estimating the fixed $\overset{{}_{\shortrightarrow}}{\mu}^{\text{truth}},\overset{{}_{\shortrightarrow}}{\alpha}^{\text{truth}},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{\text{truth}},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{\text{truth}},\bm{\eta}^{\text{truth}}$ as well as fitting randomly simulated draws for $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}$. Parameter $\overset{{}_{\shortrightarrow}}{\alpha}$ tends to be under-covered, as it is a hierarchy parameter meant to describe variation in $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}$, which are unobserved latent data in the model. On this simulation, where $\vec{\mu}$ was fixed, the hyper parameters hyper parameters such as $\tau^{2}_{\eta},\alpha_{\mbox{\tiny{all}}}$ are undefined.

4.1 Simulation from an Alternate Model

Arguably, the beta-distribution model may not reflect the true mechanism generating the data. We simulate data an alternate linear mixed-effects model on the logit scale:

[TABLE]

In the model 12, observed $Y_{ij}$ are Gaussian on the logit scale, with a parameter $\mu_{g}$ to represent cross means, $\delta_{i}$ is the $i^{\mbox{\tiny{th}}}$ individual’s average deviance from the population mean, and each tissue-gene $j$ supplies a different sized bias and measurement noise. We use $\mbox{Bias}_{j}=\log_{e}\frac{R_{j}^{\text{truth}}}{1-R_{j}^{\text{truth}}}$ using values of $R_{j}^{\text{truth}}$ from the simulation above, and use $\sigma^{2}_{j}=1/S_{j}^{\text{truth}}$. The $\delta_{i}$ are simulated from a $N(0,1/\alpha_{g}^{\text{truth}}=1/50)$ distribution.

We fit our original Bayesian model to data generated from the alternate model, and then simulate data from the posterior predictive distribution to give the Bayesian model’s estimate for $\bar{Y}^{\mbox{\tiny{pop}}}_{g}$. In Table 8, we compare, as in Table 6, performance of the three estimators in this case where the model assumption for the Bayesian method is incorrect. The Bayesian method has the shortest intervals and yet is conservative, while the sample mean and bootstrap confidence intervals can under-cover.

4.2 “Weight-Biased Coin” Model

We simulate data from the “weight-biased coin” model (WBC) to show our performance in measuring $\overset{{}_{\shortrightarrow}}{a}$. We fix 6 alleles with additive strengths $\overset{{}_{\shortrightarrow}}{a}^{\text{truth}}=\{-5,-1,0,1,2,3\}^{T}/8$, but there will be no parent-of-origin effect, $\overset{{}_{\shortrightarrow}}{m}^{\text{truth}}=\{0,0,0,0,0,0\}^{T}$. We produce $n=320$ individuals per simulation, which will be in groups of size $40$ from a random selection of 8 crosses from the $6\times 6-6=30$ possible crosses. We simulate a study design that guarantees each allele is sampled at least once, but do not guarantee that an allele is featured both as a dam and a sire. In the real experiment, allele membership was not known or anticipated until the crosses were performed. We use simulation to judge whether a random selection of crosses of this type can properly differentiate $\alpha_{k},\alpha_{k^{\prime}}$. We use the same $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R}^{\text{truth}},\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}^{\text{truth}},\bm{\eta}^{\text{truth}},\overset{{}_{\shortrightarrow}}{\alpha}^{\text{truth}}$ from before and verify that confidence intervals cover with desired accuracy.

Table 9 shows that, while coverage is a proper 95%, the confidence intervals are quite wide, approximately $.6$ or $\pm.3$ for every allele effect. This is disheartening; adjacent allele effects are rarely distinguished in the posterior. HPDs for $\overset{{}_{\shortrightarrow}}{a}$ are affected by uncertainty for $m_{k}$ reciprocal effects, which is large when a strain has not served at least once as both dam and sire. For this reason, we retest in Table 10 in the same simulation framework, but where $\overset{{}_{\shortrightarrow}}{m}$ is set to zero. This produces narrower HPDs by 2/3 allowing adjacent $a_{j},a_{j^{\prime}}$ to be distinguished.

Parent-of-origin effects do, however, appear to be present in the real data, indicating nonzero $\overset{{}_{\shortrightarrow}}{m}$. We conclude that a study would need more than 8 crosses to study these alleles and demonstrate in Table 11 that 16 of 30 crosses is sufficient. In the real dataset we collected data from 12 out of 20 allele combinations We used our simulations to encourage experimenters to perform additional reciprocal crosses.

5 Data Analysis

The Calaway et al. (2013) dataset comprises 660 mice, with 24 tissue-gene combinations, for brain, liver, and kidney. Missingness is such that for 15840 possible measurements, only 2393 could be observed; thus 84.9% of potential $\mathbf{Y}$ is unobserved. In addition, 34 mice, with 10 brain-derived RNA-seq gene measurements are provided from a set of WSB, PWK, CAST wild-derived $\text{F}_{1}$ crosses. These crosses were measured on a different set of 10 genes and no pyrosequencing measurements were taken. Using both our IA and WBC models for $\overset{{}_{\shortrightarrow}}{\mu}$, we implement Gibbs samplers using the observed $\mathbf{Y}$ to draw from the posterior of $\overset{{}_{\mbox{\hskip 1.42262pt}\shortrightarrow}}{P}$, $\overset{{}_{\shortrightarrow}}{\mu}$, and $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$. In Figure LABEL:fig:ShrinkPlot, we demonstrate how the posterior mean estimates for $P_{i}$ compare to observed $Y_{ij}$ values, using a shrinkage diagram (Efron and Morris, 1977). Since estimates of $P_{i}$ are influenced by $\mu_{g}$ values, fitted values, $\hat{P}_{i}$, are different from the unweighted arithmetic mean, $\sum_{j}Y_{ij}/J_{i}$, of the observed gene expressions. Gene measurements with larger $S_{j}$ values carry higher weight in the estimate of $P_{i}$.

For three MCMC chains of length 7500 computed in parallel over the course of a few hours, with i.i.d. $\mbox{Beta}(5,5)$-distributed starting values, the maximum 95% quantile Gelman-Rubin convergence diagnostic (GRD) is 1.11 for the 32 $\vec{\mu}$ parameters with a multivariate potential scale reduction factor (PSRF) of 1.08 and a mean 5-lag autocorrelation .267 (Gelman and Rubin, 1992; Brooks and Gelman, 1998). In contrast, the maximum 95% quantile GRD is 1.06 for the 34 $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ parameters with $\mbox{Gamma}(10)/10$ distributed starting values, and a multivariate PSRF of 1.04, but the mean 5-lag autocorrelation is a much slower .784. We conclude that mixing of the chains is sufficient for analysis of this dataset.

The most important scientific objectives are accomplished by estimation of $\overset{{}_{\shortrightarrow}}{\mu}_{g}$. Most strains were anticipated to carry $\mbox{\emph{Xce}}^{b}$, and so we used three pieces of evidence to converge this suspicion. We established these candidates carried alleles stronger than $\mbox{\emph{Xce}}^{a}$ by showing that $\mu_{g}>.5$ when these candidates served as dam when crossed with known $\mbox{\emph{Xce}}^{a}$ carriers, we established these candidates carried alleles weaker than $\mbox{\emph{Xce}}^{c}$ by showing that $\mu_{g}<.5$ when they served as dam when crossed with $\mbox{\emph{Xce}}^{c}$ carriers, and finally we tried to establish that $\mu_{g}$ was close to $.5$ when these candidates were crossed with known $\mbox{\emph{Xce}}^{b}$ carriers. For crosses these candidates served as sire, the ordering of $\mu_{g}$ is reversed. We rejected undetermined candidates as allele members when crosses with preexisting allele holders demonstrated $\mu_{g}\neq.5$. We used a criterion of $2\times\mbox{min}(\mathcal{P}(\mu_{g}>.5|\mathbf{Y}),\mathcal{P}(\mu_{g}<.5|\mathbf{Y}))$, that is, twice the smallest tail posterior probability on either side of $.5$. An example of posterior estimation for $\mu_{g}$ in crosses with ALS is given in Figure 3, and a summary table of crosses of interest is in Table 12. These give estimates for $\mu_{g}$ and confidence measures, for crosses of the unknown strains with known “$\mbox{\emph{Xce}}^{a}$”, “$\mbox{\emph{Xce}}^{b}$”, “$\mbox{\emph{Xce}}^{c}$” carriers. When multiple crosses between two alleles have been conducted, we report $p_{g}=2\mbox{min}(P(\mu_{g}>.5),P(\mu_{g}<.5))$ for the lowest posterior value among those crosses.

When multiple crosses, including reciprocals, against an allele carrier were performed, we present the estimates with the smallest $p_{g}$. A “*” suggests that we have made a scientific decision to associate this strain with this allele. **: ALS was tentatively called a “$\mbox{\emph{Xce}}^{b}$” since the ALS $\times$B6 cross had $\mu_{g}=0.46$ $(0.43,0.51)$ with $p_{g}=0.092$ for $30$ animals, despite the reciprocal cross showing such a strong parent-of-origin effect (also for $30$ animals).

Operating on this principle we confirmed that LEWES, PERA, TIRANO, SJL, WLA, WSB, and ZALENDE carry $\mbox{\emph{Xce}}^{b}$. Tentatively, we believe ALS to be a $\mbox{\emph{Xce}}^{b}$ carrier, since it is not of type $\mbox{\emph{Xce}}^{a}$ or $\mbox{\emph{Xce}}^{c}$. While B6 $\times$ALS cross showed bias from $.5$, the ALS $\times$B6 seemed much closer to a cross between identical carriers. PWK appears to have an allele with strength near that of the $\mbox{\emph{Xce}}^{b}$ allele, although previous literature hypothesized that it contains a $\mbox{\emph{Xce}}^{c}$ as does CAST.

5.1 Analysis of Weight-biased Coin Model

Table 13 displays logit-scale estimated allele effects based upon our decisions for carriers of $\mbox{\emph{Xce}}^{a}$, $\mbox{\emph{Xce}}^{b}$, $\mbox{\emph{Xce}}^{c}$, as well as yet-unidentified placeholders for PWK ($\mbox{\emph{Xce}}^{p}$), and PANCEVO ($\mbox{\emph{Xce}}^{v}$). We observed that parent-of-origin effect is not necessarily attached to Xce allele, so we separated strains into potential parent-of-origin carriers based upon phylogenetic relationships between strains. We model nonzero $m_{k}$ parent-of-origin only for candidates that serve both as dam and sire, else we fix $m_{k}=0$. .

The WBC model measures B6 to have a large .46 parent-of-origin effect $m_{\mbox{\tiny{{B6}}}}$, and $m_{\mbox{\tiny{{WSB}}}}$ for WSB is large and negative. Wild-derived mus musculus representative PWK appears to carry an additive allele with effect strength between $\mbox{\emph{Xce}}^{a}$ and $\mbox{\emph{Xce}}^{b}$. Results for the weight-biased coin (WBC) model were not considered essential to the Calaway et al. (2013) paper, but this model can be incorporated with genetic sequence experiments to search for parent-of-origin loci.

5.2 Comparison of Weight-Biased Coin versus Independence Assumption

Comparing posterior means $\hat{\mu}_{g}$ from the WBC and the original IA model for $\mu_{g}$, the mean $\|\hat{\mu}_{g}^{\mbox{\tiny{WBC}}}-\hat{\mu}_{g}^{\mbox{\tiny{IA}}}\|_{2}$ is .12 and the correlation $\mbox{Corr}(\hat{\mu}_{g}^{\mbox{\tiny{WBC}}},\hat{\mu}_{g}^{\mbox{\tiny{IA}}})$ is $.63$. The cross “HF” which is a WSB $\times$CAST cross measured using RNA-seq shows the maximum difference with $\hat{\mu}_{g}^{IA}=0.42\mbox{ }(.297,\mbox{ }.55)$ and $\hat{\mu}_{g}^{WBC}=0.16225\mbox{ }(.053,\mbox{ }0.289)$.

Figure 4 demonstrates the differences between WBC and IA models. Points sizes indicate the number of individuals of each cross, varying in size from 2 to 43. Larger crosses contribute more to the likelihood, and the WBC model will sacrifice fit for the smaller crosses. For the WBC model, strains used only as dam or sire have fixed $m_{k}=0$. The green points in Figure 4 are between strains which both have non-zero $\hat{m}_{k}$, whereas the red points include one $\hat{m}_{k}$ contribution, and the purple points have none.

HPD intervals in the WBC and IA models do not overlap for only 6 crosses: HF (WSB $\times$CAST), AG (AJ $\times$ PWK), BZ (B6 $\times$ ZALENDE), SB (SJL $\times$ B6), BT (B6 $\times$TIRANO) and BO (B6 $\times$PANCEVO). Additionally, in crosses GH (PWK $\times$WSB), HG (WSB $\times$PWK), HA (WSB $\times$AJ), GA (PWK $\times$AJ), FH (CAST $\times$WSB) the WBC model seems to underfit the observed result, though these crosses are smaller and HPD intervals overlap.

Wide disparities between $\hat{\mu}^{\mbox{\tiny{WBC}}}$ and $\hat{\mu}^{\mbox{\tiny{IA}}}$ may suggest crosses have been insufficiently sampled, or gene measurements might be biased, Xce alleles are mislabeled, or that the logit-linear model is insufficient. But in general, the WBC model successfully predicts $\mu_{g}$.

5.3 Precision parameters

Tables 14 and 15 presents estimated $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$ precision parameters for the IA model. Some RNA-seq gene measurements (taken only in the brain), have similar $S_{j}$ levels as pyrosequencing. Ideally, we would have liked to measure both pyrosequencing and RNA-seq on the same tissue, same mouse, same cross same gene. But we were only able to procure RNA-seq from a set of PWK-WSB-CAST crosses sequenced to serve multiple experiments in addition to our own. In the pyrosequencing measurements, a gene seems to be more precise in the brain than in the kidney or liver. Estimates for $\hat{S}_{j}$ range from 300 (suggesting a s.d. of about $\pm.03$ near $P_{i}=.5$) to some RNA-seq measurements that only display $\hat{S}_{j}=8$ (an s.d. of about $\pm.17$) accuracy.

Based upon the Pòlya Urn abstraction described earlier, we can interpret these values suggesting multiple hundreds of cells in the brain at the point of X-inactivation, relative to 50 cells in the kidney or liver. Measurements of Xist, Hprt1, Acs14 measurements appeared in few crosses, and credibility for $S_{j}$ of these genes range from $10^{-13}$ to $300$, meaning that these genes were not measured in enough crosses provide reliable contribution.

6 Discussion and Conclusions

We have developed a Bayesian hierarchical model, robust and powerful enough to confirm our hypotheses, which uses allele-specific gene-expression to estimate whole-body X-inactivation. In simulations, we showed that Bayesian 95% HPD credibility intervals reliably cover the truth in simulated datasets similar in size and structure to our experimental dataset, and that these intervals improve coverage relative to asymptotic-$t$ and bootstrap-based estimators, even when data is generated from an alternate model. Our “independence-assumption” (IA) model, where $\mu_{g}$’s are given i.i.d. priors, demonstrated and confirmed allele membership of 10 previously untested strains, which was the critical analysis of the Calaway et al. (2013) experiment. Our weight-biased coin (WBC) model, where $\mu_{g}$ is a function of parental strains, including additive allele effects $\overset{{}_{\shortrightarrow}}{a}$ and reciprocal effects $\overset{{}_{\shortrightarrow}}{m}$, produced similar results to the IA prior, suggesting this scientific model might be used to predict $\mu_{g}$ in future crosses.

Modeling bias and precision of gene-specific measurements using free ($\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{S}$) and constrained $(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{R},\bm{\eta})$ parameters, we developed techniques in constrained slice-sampling that required both linear and non-linear transformation of the parameter set. Our model could accommodate allele-specific expression measured both by pyrosequencing and RNA-seq.

Here we must acknowledge limitations in experimental design that led to Bayesian recourse. Quantitative trait loci (QTL) are regions within the genome whose variation is highly correlated, and conjectured to be causal, with a phenotype — in our case X-inactivation skewing. Typical experiments to locate QTL, such as Advanced Intercrosses (Darvasi and Soller, 1995), in which individuals of known genotype are mated, rely on random recombination of the genome. But precise measurement of Xce bias requires replicates in the form of cloned individuals. Calaway et al. (2013)’s methodology resulted from fortuitous selection of already-sequenced inbred strains that were found to exhibit sequences of similarity and difference within the Chadwick et al. (2006) interval. Establishing that two strains have heavy difference in a sub-region, yet their offspring exhibit no X-inactivation skewing, gives credence to a belief that the sub-region of difference is not a candidate for Xce. Establishing that two strains, identical at a sub-region, have offspring that exhibit X-inactivation skewing, rejects a sub-region in favor of other regions where these two strains do differentiate genetically.

Prior to experiment we did not know all of which crosses to perform, did not know which sub-region we wished to justify, and did not know which strains would differentiate. As early strains in the experiment were established as carriers for known alleles “$a$”, “$b$”, “$c$”, we were further motivated to test other strains that continued to break down the candidate region. Such experimental choices of new test strains were done both in the light and in the dark of statistical analysis and estimates; as such, our stopping rule cannot be claimed to be a sample-independent choice. The “null” sampling space in this case is ill-posed: we did not randomly sample from a population, but instead chose sub-populations to explore based upon research intuition. Sometimes samples were obtained because they were readily available or cheaper to obtain. But, at other times, crosses were ordered because they were hoped to lead to a new allele. Sometimes additional samples were obtained to improve statistical confidence on a cross, adding to samples already collected. From a frequentist perspective, it is right to be skeptical of analysis on such a dataset. Since we cannot map out the decision metric that led to collection of samples, we cannot postulate the alternative sample space for most hypothesis tests, and cannot report a $p$-value.

Bayesian exchangeability  (de Finneti, 1937; Bernardo and Val, 1996; Lindley and Phillips, 1976; Hewitt et al., 1955; Jackman, 2009; Diaconis and Freedman, 1980) often justifies posterior analysis of sequential designs, and here we show the sampling mechanism is an ignorable design (Thompson and Seber, 1996). Consider indicators $I_{i,1},I_{i,2},\ldots,I_{i,G}$ which represent whether individual $i$ is in group $g\in\{1,2,\ldots,G\}$, such that $\sum_{g}I_{i,g}=1$. Let $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i}$ vector be assigned with a transition probability f^{i}_{\text{trans}}(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i}|Y_{\mbox{\tiny{<i}}},\mathbf{I}_{\mbox{\tiny{<i}}}), which is a multinomial dependent only upon previous observations and sampling choices. Conditional on $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i},\overset{{}_{\shortrightarrow}}{\mu},\overset{{}_{\shortrightarrow}}{\alpha}$, the $P_{i}$ are independent from other individuals. The full posterior, including arbitrary sequential design, is

[TABLE]

We see from expanded equation  13 above that because of conditional independence given $P_{i}$, and because all vectors $\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i}$ are observed and known, then the transition density, f_{\text{trans}}^{i}(\overset{{}_{\mbox{\hskip 1.42262pt}\shortleftarrow}}{I}_{i}|Y_{\mbox{\tiny{<i}}},\mathbf{I}_{\mbox{\tiny{<i}}}), which is only based upon observables, cannot influence the posterior. Because experimenters had no preference before observations of allele identities for the individuals chosen, we can accept that the stopping rule is proper (though certainly not optimal) in the definition of a proper stopping rule given in  Parmigiani and Inoue (2009). We would have eventually stopped collecting samples, no matter the true value of $\mu_{g}$.

We note, per a criticism in Lindley and Smith (1972), that our $\overset{{}_{\shortrightarrow}}{\mu}$ parameters are not actually exchangeable based upon science. Even if we do not know the value of $\mu_{g}$, we do know that if $\mu_{g},\mu_{g^{\prime}},\mu_{g^{\prime\prime}}$ were a sequence of crosses of a mother strain with unknown allele “UNK” to fathers of allele carriers “a”, “b”, “c”, then there should be a decreasing sequence: $\mu_{g}$ $\geq$ \mu_{g^{\prime}}$$\geq $\mu_{g^{\prime\prime}}$. In the “IA” model, we did not model allele-order hypothesis with our priors and relied upon the data to reveal this sequence. The WBC model did enforce this relation, and reached similar estimates of $\hat{\mu}_{g}$.

It is difficult to statistically argue for the hypothesis $\mu_{g^{\prime}}=\mu_{g}$, in contrast to proving $\mu_{g^{\prime}}\neq\mu_{g}$. We have relied upon a criterion of overlapping posterior. It is possible that the difference between two cross means, $\mu_{g}-\mu_{g^{\prime}}=\Delta$, might be small, non-zero, but imperceptible, such as $\Delta\approx.0000001$. In which case, we would inevitably decide incorrectly to treat $\mu_{g},\mu_{g^{\prime}}$ as fundamentally equivalent crosses between same alleles. We assume that alleles differ in effect enough to be statistically observable.

As high-throughput sequencing becomes cheaper, Bayesian algorithms of our design must be rejected in favor of more efficient methods. The Bayesian algorithm has complexity $\mathcal{O}(NJ^{2})$, due to Gibbs sampling of individual $S_{j}$ and $R_{j}$, making it unsuitable for whole-transcriptome datasets, where a complete genome has $J\approx 8000$. At the time of research, RNA-seq of a single brain-tissue sample might cost $1,000 and require multiple weeks of setup, sequencing, and bioinformatic analysis, but pyrosequencing a few genes, to greater precision, of the same sample costs$\approx$$20. A fully-analyzed RNA-seq sample might measure an individual’s $P_{i}$ average brain expression to a $\pm.1\%$ region of 95% credibility, but pyrosequencing of three genes would leave a $\pm 2\%$ region. But to measure cross mean $\mu_{g}$, it is better to have $P_{i}$ measured imprecisely in many individuals in multiple tissues rather than have high precision on any single-tissue $P_{i}$.

We note that point-estimate performances of the sample mean estimator and bootstrap are nearly identical to the posterior mean. Were it not for under-coverage in their confidence intervals, and a conceptual disconnect between these procedures and our parametric model, these methods might have served as suitable replacements to the Bayesian Gibbs sampling.

To the extent we have introduced techniques more generally applicable to statistics, our constrained priors show promise in other naturally unidentifiable problems with small order $J$. We are fortunate to be in a data-setting where a Bayesian estimator is statistically robust and viable, and can use posterior sampling to test a parametric model for its comparable nuances.

7 Acknowledgements

This project was supported by National Institutes of Health (NIH) grants R01GM104125 (ABL, WV), R35GM127000 (WV), and P50MH090338 and P50HG006582 (FPMdV, JDC). We also thank UNC Information Technology Services for computational support. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of this manuscript