Bayesian Manifold-Constrained-Prior Model for an Experiment to Locate Xce
Alan B. Lenarcic, John D. Calaway, Fernando Pardo-Manuel de Villena, William Valdar
BBayesian Manifold-Constrained-Prior Model for an Ex-periment to Locate
Xce
Alan B. Lenarcic
Department of Genetics, University of North Carolina at Chapel Hill, USA.
John D. Calaway
Department of Genetics, University of North Carolina at Chapel Hill, USA.
Fernando Pardo-Manuel de Villena
Department of Genetics, University of North Carolina at Chapel Hill, USA.
William Valdar
Department of Genetics, University of North Carolina at Chapel Hill, USA.
Summary.
We propose an analysis for a novel experiment intended to locate thegenetic locus
Xce (X-chromosome controlling element), which biases the stochasticprocess of X -inactivation in the mouse. X -inactivation bias is a phenomenon wherecells in the embryo randomly choose one parental chromosome to inactivate, butshow an average bias towards one parental strain. Measurement of allele-specificgene-expression through pyrosequencing was conducted on mouse crosses of anuncharacterized parent with known carriers. Our Bayesian analysis is suitable for thisadaptive experimental design, accounting for the biases and differences in precisionamong genes. Model identifiability is facilitated by priors constrained to a manifold.We show that reparameterized slice-sampling can suitably tackle a general class ofconstrained priors. We demonstrate a physical model, based upon a “weighted-coin”hypothesis, that predicts X -inactivation ratios in untested crosses. This model sug-gests that Xce alleles differ due to a process known as copy number variation, wherestronger
Xce alleles are shorter sequences.
1. Introduction X -inactivation bias is a phenomenon observed in the female mouse, first proposedin Cattanach and Isaacson (1967) to be caused by a gene Xce (X-controlling ele-ment), which has unknown location and mechanism. In mammals, cells of the earlyfemale (XX) embryo inactivate one of their parental chromosomes through a ran-dom choice, and pass this decision onto their daughter cells (Lyon, 1962; Gendreland Heard, 2011), as shown in Figure 1. Thus, in mammals, tissues have mosaicregions of parental chromosome preference; this can be seen, for example, in the furof 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 withparents of white and black hair follicle color, counting the ratio of colors in the hy-brid animal, and observing an average bias toward one color. In our related report, E-mail: [email protected]: [email protected]: [email protected]: [email protected] a r X i v : . [ s t a t . A P ] D ec Alan Lenarcic et. al
Calaway et al. (2013), we proposed to investigate for the location of
Xce , refiningan interval proposed in Chadwick et al. (2006), by ascertaining the allele carried by10 unclassified strains. For experimental reasons relating to cost and accuracy, thatstudy measured allele-specific expression in infant hybrids using the technology ofpyrosequencing.
1/ 0Day 0: FertilizationOnly maternal X expressed l lllllllll lll llll llll lllllllll lll llll lll l lllllll llll llll lll llllllllllll llll lll l lll lllll llll ll lllllll l lllll lll llllllllllll l lllllll llll l lll llll lllllll llll llll lll llllllllllll llll lll l lll lllll llll ll lllllll l lllll lll llllllllllll l lllllll llll l lll lll lll llll lll lll lllllllll ll l lll llll ll ll llll llll l lll lll l l ll lll llllllllllllll ll llll ll llllllllll l llll lll ll
248 252/Day 5: X− InactivationChoice of one X for every cell 10 Day 30: BirthOrganism with tissuedifferentiated regions of homologousX− inactivation choice
Fig. 1.
Schematic of X -inactivation in embryonic development (Lyon, 1961). Near day5 (Mak et al., 2004; Takagi et al., 1982), progenitor cells inactivate either their maternal orpaternal X -chromosome, turning it into a Barr Body(Barr and Bertram, 1949; Lee, 2011),later passing this choice to future daughter cells. The full organism at birth is a mosaic of re-gions of predominant maternal or paternal expression, under usual circumstances expectedat an overall 50-50 balance. This balance is disturbed if there is a survival disadvantage forcells expressing either X -chromosome, or if there was an Xce bias influencing the initial se-lection. In this paper we measure X -linked gene expression on multiple mouse pup tissues,with one goal of inferring the X -inactivation proportions initially set on day 5. 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 BayesianMCMC-based hierarchical model, based upon the beta-distribution. The modelrelates allele-specific proportions to a summary proportion reflecting the initialproportion of progenitor cells that conducted the initial choice around 5 days in theembryo 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 useposterior tail values to argue for the identity of the allele carried by an unknownstrain, based on its average bias against known allele-carriers.To make our model identifiable, we propose priors for sets of our parametersthat are constrained to a manifold condition. Through reparameterizing the slice-sample distribution, we can generate acceptable draws from the correct, constrainedposterior. Constrained priors become especially necessary when we attempt todesign a physical model for the X -inactivation bias observed in the crosses. Topredict 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 thefindings of the Calaway et al. (2013)
Xce experiment and demonstrate value forBayesian hierarchical modeling in future allele-specific-expression experiments. eta model for X-inactivation
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, forexample, 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 andIsaacson, 1967; Cattanach and Rasberry, 1994, 1991). We have related these alleleswith respect to the order of their “strength”. An AJ × B6 hybrid which is a × b will favor expressing the B6 X -chromosome on average in approximately 65% oftheir cells; an a × c cross produces approximately a 30-70 ratio; for a × e it is20-80. Even when the population average for a given hybrid is 50-50, there isconsiderable variation in individual “skew”, i.e. the proportion an individual mightover-express one parental X -chromosome. Skewed ratios of 75-25, 27-72, etc. areobserved in individuals, and the allele specific expression individual genes mightpresent additional imprecision or biases.Chadwick et al. (2006) proposed a location of Xce to somewhere within a ∼ X -chromosome. This was based upon identi-fication of the allele carried by strains of mice, originally conducted with techniqueslike shaving and counting hair follicles of hybrids with differently colored parents.But counting hair follicles imprecisely measures X -inactivation only in mice oldenough to be shaven and only in skin tissues. Given the costs involved in studyingembryos and laboratory expertise, Calaway et al. (2013) sought an efficient andreliable measure of aggregate whole-body expression taken from newborn tissue.Due to variation of X -inactivation skew among clones of identical genome, multipleindividuals must be measured to establish that a given “first-generation hybrid”(F ) cross shows a bias, favoring measurement techniques of low cost.Pyrosequencing, a “sequencing by synthesis method” introduced in Ronaghiet al. (1998), was viewed as a cost-effective technology to measure gene expres-sion for strains whose genome sequences were known but whose Xce alleles werenot. Calaway et al. (2013) chose commonly expressed X -located genes: Ddx26b , Fgd1 , Rragb , Mid2 , Zfp185 , Xist , Hprt1 , and
Acsl4 , which featured known SNPsthat differentiate two strains. Pyrosequencing reads target-sequences of length 20-30 base pairs including either allele of a SNP: Exponential range PCR amplifiesup to 1 micro-liter ( µl ) of the target, or approximately . × copies, which,after purification, are counted through synthesis. Pyrosequencers report measure-ments in proprietary fluorescence units, where, conservatively, the final ratio ofallele-specific gene expression is based upon more than million targeted reads.Table 2 gives a short reference background to the mouse types used in thisexperiment. There are three mouse sub-species: mus musculus domesticus fromEurope, 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 strainswere initially generated from domesticated toy mice. Wild-derived strains comefrom 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 . Alan Lenarcic et. al
Table 1.
Information on common inbred mouse strains used in this experiment.
Strain Abbrev. details Sub-species derived
Xce alleleA/J AJ Albino (hearing loss) domesticus lab a 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 ??
The Calaway et al. (2013) first established which
Xce allele was carried by a few un-classified strains. Based upon the observations, a provisional, semi-statistical clas-sification was assumed, and this motivated the choice to classify additional strains.First, strains SJL/J (
SJL ) and ALS/LtJ (
ALS ) were analyzed because their se-quences split the known region between
Xce a and 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 measure-ments. Available frozen F samples of crosses including DDK/Pas ( DDK ), ZAL-ENDE/Ei (
ZALENDE ), PERA/EiJ (
PERA ) (Nishioka, 1995) were acquired. Somesuccessful breedings with PANCEVO/EiJ (
PANCEVO ) (Kim et al., 2005) were con-ducted to see if this second species of mouse had a different
Xce profile. Onlya subset of possible crosses were performed, and not all reciprocal crosses weregenerated. Count of replicates per unique crosses varied from 4 to 50.For convenience, summarizing details for the inbred parental strains are givenin Table 2. A statistical analysis must respect that the experiment is adaptive andsequential in nature, that many crosses which might have aided in the statisticalestimate will be missing, and that we should provide a statistical criterion not onlyto assert that a given cross
SJL × B6 has an average “ . ± . ”, but that weshould also be able to express how certain we are that SJL is an a allele carrier.By choosing closely related potential Xce a and Xce b allele carriers, predominantlyfrom 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 ap-proximately 500 kb from the gene
Xist whose expression causes X -inactivation.Calaway et al. (2013) established that allele Xce a is rare in the wild and thatmost individuals within a sub-species share the same Xce allele. They found that
Xce -allele strengths can potentially increase or decrease during speciation, whichsuggests
Xce is mutated through copy number variation. eta model for X-inactivation Table 2.
Example section of pyrosequencing maternal expression ratiosdataset. pup cross dam sire BrainDdx26b BrainRragb KidneyDdx26b KidneyRragb . . .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 NA7-4 1Wl 129S1-7 WLA-7 0.648 NA NA NA16-3 AlAj ALS-16 AJ-873 0.456 .501 0.447 .5116-5 AlAj ALS-16 AJ-873 0.574 .432 0.532 .59 . . .. . .
A partial Table 2.2 shows the challenges of the Calaway et al. (2013) data. Everymouse is a member of a single cross, which is coded to correspond to a uniquehybrid combination of mother and father inbred strains. We encode every crosswith a short code, such as “1Wl” for a × WLA cross, where is the mater-nal strain. Note that a
WLA × hybrid, with WLA as maternal strain, would becoded “Wl1” and coded as a different reciprocal cross. The other columns are mea-surements of maternal proportion expressed for given tissue-gene combinations.Allele-specific proportion is measured in multiple tissue-genes, with considerablemissingness: “NA” for genes that were either dropped at random for an individ-ual mouse during experiment, or because they are not experimentally viable forthe cross when no SNP exists on that gene to differentiate maternal and paternalcopies. Since SNPs are single base-pair substitutions, between any three parentalstrains there is no single SNP that can be used to measure allele-specific expressionfor all three pairings of two strains.If we naively take sample means of all observed gene expression in individualsstemming from one type of cross, we can summarize those averages in a table likeTable 2.2, 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 × a cross, have often not beenmeasured. We inevitably hope for a statistical method that can establish the allelecarried 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 haveon order 500 read counts of any SNP region, offering less precision than pyrose-quencing. A Bayesian method for estimating precision and bias with RNA-seqcount data can be found in Graze et al. (2012). Both RNA-seq and pyrosequencingare unbiased at preserving true, initial allelic stoichiometry. But the non-linear re-lationship between expression and luminescence in microarrays can lead to biasedestimation of allele-specific proportion.Late in the Calaway et al. (2013) experiment, brain-tissue RNA-seq measure-ments were made available for 34 individuals from F crosses between wild-derived Alan Lenarcic et. al
Table 3.
Unweighted sample mean expression, maternal allele(rows), paternal allele (columns). a b c ?? L ?? P ?? S a . b . . . . c . ?? L . ?? P . . . ?? S . . . . , where: a AJ , b B6 , ALS c CAST ?? L LEWES ?? P PWK ?? S SJL sub-species representatives
WSB , PWK , and
CAST , about individuals of each cross,including reciprocal crosses. It was hoped that our statistical procedure designedfor pyrosequencing output could accommodate RNA-seq data. As we show, RNA-seq measurements serve the same role in our model as pyrosequencing, althoughRNA-seq measurements often have less precision when baseline expression of a geneis low.
3. Statistical Model for Gene-expression Ratios
Denote the proportion of maternal gene-expression measurements Y ij ∈ [0 . , . , where i ∈ { , . . . , n } is the index of the mouse specimen used in theexperiment, and j ∈ { , . . . , J } is an index of tissue-gene combination (examples:“kidney- Fgd1 ”, “brain-
Ddx26b ”). We are given Y ij as a fractional proportion ofmaternal 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 inmultiple genes, we hope to estimate a quantity, P i , representing a best whole-bodyproportion of active maternal X -chromosome.Even if mouse i has a whole-body proportion of P i , replicate mice from thesame cross i (cid:48) , i (cid:48)(cid:48) , . . . will have different proportions P i (cid:48) , P i (cid:48)(cid:48) , . . . . The cross-specificmean µ g denote the population mean for mice in cross/group g ∈ { , . . . , G } . Weseek 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 µ g . Only a subset of genes j canbe 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 inferencemodel follows. Our scientific objectives are served by a Bayesian hierarchical model. Observed dataare proportions Y ij , so a beta-distribution regression model (Ferrari and Cribari-Neto, 2004) can model unimodal proportions. For our problem, the mean and modefor the beta densities should be constrained within (0 , , excluding and . Thus,we choose parameterized beta distributions of the form Beta ( A + 1 , B + 1) where A , B ≥ . We first give the parameterized mathematical model, and then describe eta model for X-inactivation the purpose of the parameters in the model. Let:Observed: Y ij ∼ Beta ( P i R j S j e η g,j + 1 , (1 − P i )(1 − R j ) S j e η g,j + 1) ,Where Latent: P i (cid:124)(cid:123)(cid:122)(cid:125) Whole Body ∼ Beta ( µ g α g + 1 , (1 − µ g ) α g + 1) . (1)In the above equation, we have latent, unobserved ratios P i for i ∈ { , . . . , n } foreach individual and parameter vectors (cid:0) R, (cid:0) S ∈ R J , (cid:1) α, (cid:1) µ ∈ R G and a perturbationmatrix η ∈ R G × J . In our notation, (cid:0) R and (cid:0) S are row-vectors each of length J ,the number of columns of Table 2.2, while (cid:1) P , (cid:1) µ, (cid:1) α are column-vectors for the n individuals or G crosses. The matrix η is of dimension G × J . R j ∈ (0 , represents bias for tissue-gene combination j ∈ { , . . . , J } . Thevalue R j = . suggests that this tissue-gene is relatively unbiased as a predictor ofoverall X -inactivation. If a tissue-gene has R j = . , then we expect this to push Y ij , on average for all crosses, in the direction of maternal expression.In contrast, S j ∈ (0 , ∞ ) , represent precision of measurements for a giventissue-gene combination j . Larger values of S j imply that tissue-gene j is a moreprecise measurement of P i . To add further difference in precision for certain crosses,parameters η g,j ∈ {−∞ , ∞} model deviation in precision for crosses g oftissue-genemeasurements j as compared to other crosses. The precision for cross g on tissue-gene j will be S j × e η g,j . Because the Beta-distribution implies a Pòlya-urn modelfor cell proliferation, we can use measurement of S j to infer in some measure howmany progenitor cells exist in each organ during the X -inactivation event. Letpriors for S j be Gamma ( χ S , ξ S ) ∝ e − x/ξ S x χ S − with prior for hyper-parameters ( ξ S , χ S ) ∼ Gamma ( . , . .The +1 to both terms of Beta ( A + 1 , B + 1) ensure stability of the posterior,keeping the beta density ∝ x A +1 − (1 − x ) B +1 − from having infinite boundarydensity if P i , µ g ≈ or if P i , µ g ≈ . Since we restrict our population of genes andindividuals to those expressing between and of the maternal X (individualsexpressing more or less are presumed either males or XO females carrying only asingle 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 samplingchoices.Hyperparameters µ g , α g represent the average and variance for P i within a crossgroup g . µ g ∈ (0 , reflects the average P i proportion of cross g , and is themost important parameter for establishing Xce membership. Nuisance parameter α g ∈ (0 , ∞ ) represents the variability of P i about µ g . For large α g values, P i variation about µ g will be less. We assign (cid:1) µ a hierarchical prior: µ , µ , . . . , µ G ∼ i.i.d. Beta ( µ All α All + 1 , (1 − µ All ) α All + 1) . (2)For Equation 2, the global mean and precision parameters, µ All , α
All are assignedweak hyper-priors, Unif (0 , and Gamma ( . , . respectively. For vector (cid:1) α , whichrepresents the spread of observed P i around their group mean µ g , we assign an i.i.d.Gamma (1 , exponential prior. Alan Lenarcic et. al
Due to identifiability, canonical i.i.d. priors do not suit our purposes for parameters (cid:0) R and η . If all J bias parameters R j deviated from . together, such as if ¯ R = J (cid:80) j R j = . , the model would shift ˆ P i to the sign ( . − ¯ R ) direction of all Y ij . Buta global movement of ˆ P i would shift all ˆ µ g . Relative differences ˆ µ g − ˆ µ g might bepreserved, but we would not obtain a precise estimate of ˆ µ g , if confounded with (cid:0) R .We explicitly assume (because we have no way of discovering otherwise) thattissue-gene combinations are “on-average” unbiased: (cid:80) j logit ( R j ) = 0 , wherelogit ( R j ) = log e R j − R j . We will report ˆ µ g given an assumption that logit ( R j ) ison average unbiased. Since most decisions will be based upon the significance ofdifferences ˆ µ g (cid:48) − ˆ µ g we feel that this assumption is reasonable. We therefore modelbias with a restricted prior of the form: R j ∼ Marginally Beta ( u R , u R ) conditional on (cid:88) j log R j − R j = 0 . (3)which we show in Section 3.3 is suitable for MCMC posterior exploration.Model parameters η g,j represent deviation about parameter S j . Since η g,j ∈{−∞ , ∞} , Gaussian priors N (0 , τ η ) might seem acceptable. Certainly, if all η g,j deviated from together, then the combination product e η g,j S j could support manyvalues of S j to retain the same overall precision. Inferentially, this weakens ourestimation ability of S j . Instead, we define (cid:80) g η g,j = 0 for every j . Thus S j represents an average precision over all groups g . The prior for η g,j will be (cid:80) g η g,j =0 , but that η g,j ∼ N (0 , τ η ) marginally.Constrained priors limit the parameter space to a lower-dimensional mani-fold. In Gaussian Gibbs-regression settings, such restrictions are commonly im-posed (Gelfand et al., 1992). Vines et al. (1996) willingly accept a 10x slowdown initeration speed to estimate regressions with constrained parameters for a mixed ef-fects model, and observed reduced autocorrelation in the samples. In non-Gaussiansettings, Zhu et al. (2011); Shi et al. (2009); Zho et al. (2010) show advances inmanifold exploration and suggest new ways to explore these parameter spaces. Wedemonstrate a slice-sampler suitable for our linear constraints in Section 3.3.Due to orthogonality in design between (cid:0) R and µ g , we can approximate that for J fixed and n g → ∞ : √ n g (ˆ µ g − µ g ) = ⇒ N (cid:0) , µ g (1 − µ g )( σ u + σ e /J ) (cid:1) + O (cid:0) /S j (cid:1) (4)where σ u ≈ α g +2( α g µ g +1)( α g (1 − µ g )+1) and σ e behaves as /S j , are logit-scale linearmixed regression model group and noise variances respectively. With unconstrainedpriors, leaving (cid:0) R confounding, it is uncertain that ˆ µ g would converge.Table 3.2 summarizes our parameters and prior choices. We next describechoices in MCMC algorithm that enables investigation of this model, includingthese constrained priors. eta model for X-inactivation Table 4.
Summary of parameters involved in our beta-regression model.
Prior Purpose µ All
Beta (1 , Overall mean expression of all crosses α All
Gamma ( . , . Beta dispersion of (cid:1) µ about µ All (cid:1) µ Beta | µ All , α
All
Mean parameters µ g for each cross g ∈ { , . . . , G } (cid:1) α Gamma (1 , Dispersion parameters α g for each cross g (cid:1) a (cid:80) k a k = 0 In “weight-biased coin” model (later Section 3.4) for (cid:1) µ , ad-ditive effects (cid:1) m (cid:80) k m k = 0 In weight-biased coin model, parent-of-origin effects (cid:0) S Gamma ( ξ S , χ S ) For tissue-genes j ∈ { , , . . . J } , models the precision ofthat tissue ξ S , χ S Gamma ( . , . Group dispersion hyper-parameters for the (cid:0) S parameters (cid:0) R (cid:80) j logit R j = 0 Models the bias of tissue-gene j toward maternal or paternalexpression relative to a true individual-level cell count P i η (cid:80) g η g,j = 0 Measures differential in precision about (cid:0) S where certaincrosses g are less accurate for tissue jτ η Inv-Chi-Square (1)
Extra dispersion parameter for η in addition to constraint. u R Gamma (1 , Additional dispersion parameter for R j to have marginal R u R − j (1 − R j ) u R − proportional marginal prior density (cid:1) P Beta | (cid:1) µ, (cid:1) α For each individual i ∈ { , . . . n } , models mean maternal cellactive proportion Y obs Observed portion of n × J data matrix Y all of proportionsof maternal gene expression Y miss Beta | (cid:1) P , (cid:0)
S , (cid:0) R, η Missing or unobservable data in Y all that is imputed Alan Lenarcic et. al If Θ represents the unknown model parameters, the product of our model-basedlikelihood, P ( Y | (cid:1) P , Θ ) × P ( (cid:1) P | Θ ) and the priors, p ( Θ ) , gives an unnormalized func-tion proportional to the desired posterior P ( Θ | Y ) : P ( (cid:1) µ, (cid:1) α, µ All , (cid:0) S, (cid:0) R, η , ξ S , χ S , u R , τ η , (cid:1) P | Y ) ∝ Prior ( (cid:1) µ, (cid:1) α, . . . , u R , τ η ) × N (cid:89) i =1 Beta ( P i ; µ g α g + 1 , (1 − µ g ) α g + 1) × J (cid:89) j =1 Beta ( Y ij ; P i R j S j e η g,j + 1 , (1 − P i )(1 − R j ) S j e η g,j + 1) . (5)We use the technique of slice-sampling (Roberts et al., 2012; Neal, 2003) to achieveMCMC random draws from the posterior distribution P ( (cid:1) µ, (cid:1) α, . . . | Y ) . Let group g ∈ { , . . . , G } , and let \ g be G − { g } , or the set { , , . . . , g − , g + 1 , . . . , G } .In slice-sampling, a new parameter α ( t +1) j is drawn given current α ( t ) g as wellas fixed (cid:1) µ ( t ) , (cid:1) α ( t ) \ g , (cid:1) P ( t ) , by considering the marginal, one dimensional function: f α g ( α g | (cid:1) α ( t ) \ g , (cid:1) µ ( t ) , (cid:1) P ( t ) ) , which is the marginal unnormalized posterior at a givenpoint α ( t +1) g with other parameters fixed. In this case, given (cid:1) P and (cid:1) µ , (cid:1) α is inde-pendent of the µ All , α All , (cid:0) S , (cid:0) R , η , ξ S , χ S , u R , τ η parameters.To perform slice-sampling, first a uniform value U is drawn from U ∼ Uniform (cid:18) , f α g ( α ( t ) g | (cid:1) α ( t ) \ g , (cid:1) µ ( t ) , (cid:1) P ( t ) ) (cid:19) . (6)Then the slice sampler seeks left and right from point α ( t ) g to find the points α − below and α + above s.t. f α g ( α − | (cid:1) α ( t ) \ g , (cid:1) µ ( t ) , (cid:1) P ( t ) ) = f α g ( α + | (cid:1) α ( t ) \ g , (cid:1) µ ( t ) , (cid:1) P ( t ) ) = U . Tocomplete the draw, a second uniform is drawn from the interval [ α − g , α + g ] . Thisdraw serves as the new draw of α ( t +1) g . Looping this procedure for all parametersserves as the MCMC scheme. The algorithmic order can be understood to reflectthe order of unobserved parameters in Table 3.2.This procedure is well-known and has an advantage over Metropolis-Hastingsprocedures in not requiring tuning or a proposal density. Slice-sampling is the back-bone to the all-purpose Bayesian integration software “JAGS” (Plummer, 2003;Murphy, 2007). For our purposes, we reimplemented this process explicitly inR SHLIB-compiled C (Developers, Becker, Chambers, M., and Wilks, Developerset al.) for two reasons. First, we sought to gain efficiency in coding f ( ·|· ) densitiesby separating out conditionally independent parameters: for instance, conditionalon (cid:1) P , the posteriors for (cid:0) S and (cid:1) µ are independent. Efficiencies could be gained byusing BLAS (Lawson et al., 1979) vector functions: conditional on Y , the poste-rior for P i is proportional to the vector dot products (cid:80) j S j e η g ( i ) ,j R j log e Y ij and (cid:80) j S j e η g ( i ) ,j (1 − R j ) log e (1 − Y ij ) , which are tabulated faster with block memorymanagement. eta model for X-inactivation But, more crucially, we needed to modify the slice-sampler method to allowfor multivariate constrained prior distributions on (cid:0) R and η . Recalling that werequire (cid:80) g η g,j = 0 , we draw new samples for η ( t +1) g,j not by slice-sampling along theunivariate density f η g,j ( η g,j | (cid:0) R ( t ) , (cid:0) S ( t ) , η ( t ) ,/j , (cid:1) P ( t ) , Y ) but along the full multivariateposterior density as defined along a vector, (cid:1) η ( t ) ,j ( g, ∆) , which is perturbed by a ∆ G − ( Ge g − (cid:126) , where e g is the standard g th basis vector in R G and (cid:126) is a G-lengthvector of 1’s. In other words: (cid:1) η ( t ) ,j ( g, ∆) = (cid:18) η ( t )1 ,j − ∆ G − , η ( t )2 ,j − ∆ G − , . . . , η ( t ) g,j + ∆ , . . . , η ( t ) G,j − ∆ G − (cid:19) (7)In this case our notation considers perturbation of a column vector (cid:1) η ,j ∈ R G which represents all η g,j values for a single tissue-gene j fixed and the crosses g ∈ { , , . . . , G } .To perturb (cid:0) R , first define a logit basis reparameterization L j ≡ log e R j − R j , andtake a perturbation, ∆ J − ( Je j − (cid:0) , in the (cid:0) L logit space.It is important to note a crucial difference in the densities used in slice-samplingof η versus sampling of (cid:0) R ( (cid:0) L ) . One cannot slice-sample over arbitrary reparame-terizations without changing measure. For instance, in a univariate case, if e ω = θ ,then, if q θ ( θ ) was the target density then q ω ( ω ) = q θ ( e ω ) × e ω . It would be inappro-priate to find the two ω + , ω − such that f θ ( e ω + ) = f θ ( e ω − ) = Unif (0 , f θ ( e ω ( t − )) and then sample ω ( t ) = Unif ( ω − , ω + ) . One would see a dramatically downward-biased distribution for θ ( ω ) . For transformations θ = ω − or θ = 4 ω , slice-samplingfrom f θ ( θ ( ω )) or f θ ( θ ) produces acceptably distributed samples because the lineartransformation does not change the relative measure along dθ .In the multivariate setting, consider parameterization ω ∈ R p (cid:48) that transformslinearly onto the desired parameter space θ ∈ R p with p (cid:48) ≤ p , such as θ = A ω + (cid:126)b .Let A be the top p (cid:48) × p (cid:48) square matrix of A and A be the remainder ( p − p (cid:48) ) × p (cid:48) rectangular matrix of A . If matrix A is invertible, then θ p (cid:48) or the first p (cid:48) coefficients of θ determine θ ( p (cid:48) +1): p which are the last p − p (cid:48) coefficients of θ .If q θ ( θ ) were a target density of f θ ( θ ) of unknown integration constant, subjectto constraints, then q θ ( θ ) = f θ ( θ p (cid:48) , θ ( p (cid:48) +1): p ) (cid:82) f θ ( θ (cid:48) p (cid:48) , θ (cid:48) ( p (cid:48) +1): p ) dθ (cid:48) ...dθ (cid:48) p (cid:48) . But also, q ω ( ω ) = f θ ( A ω + (cid:126)b p (cid:48) , A ω + (cid:126)b ( p (cid:48) +1): p ) | A | (cid:82) f θ ( A ω (cid:48) + (cid:126)b p (cid:48) , A ω (cid:48) + (cid:126)b ( p (cid:48) +1): p ) | A | dω (cid:48) . . . dω (cid:48) p . (8)In Equation 8, the Jacobian | A | is part of the integration constant. Slice-samplingwill automatically evaluate this constant.In the case of (cid:1) η ,j , the rotation A is the matrix of − / ( G − (cid:126) (cid:126) T + GG − I in R G − . If A = 1 / ( G − (cid:126) and (cid:126)b = (cid:126) then the reparameterized space (cid:1) η (cid:48) ,j ∈ R G − generates (cid:1) η ,j = { A (cid:1) η (cid:48) ,j , A (cid:1) η (cid:48) ,j } and we need not change measure.In the case of (cid:0) R , which includes a nonlinear transformation (cid:126)L → (cid:0) R , the Jacobian Alan Lenarcic et. al matrix is of more consequence; (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂∂L e L e L ∂∂L e L e L . . . ∂∂L J e L e L ∂∂L e L e L ∂∂L e L e L . . . . . .. . . . . . . . . . . . ∂∂L J e L e L . . . . . . ∂∂L J e LJ e LJ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = J (cid:89) j =1 R j (1 − R j ) . (9)The change of variables from (cid:0) R to (cid:0) L measure is f (cid:0) L ( (cid:0) L | Y , (cid:0) S, (cid:1) P , η , u R ) = f (cid:0) R ( (cid:0) R ( (cid:0) L ) | Y , (cid:0) S, (cid:1) P , η , u R ) (cid:89) j R j ( L j ) (1 − R j ( L j )) . (10)As with (cid:1) η ,j , a linear transformation from (cid:0) L (cid:48) ∈ R J − , such that (cid:0) L J − = A (cid:126)L (cid:48) and (cid:0) L ( p (cid:48) +1): p = A (cid:0) L (cid:48) , maintains the (cid:80) j L j = 0 constraint. Fig. 2.
A two-dimensional example of the issues faced using a change in variablesfor slice-sampling. In the original plot, a Gaussian posterior density in 2 dimensions isconsidered, but it is desired to sample new ( X, Y ) along the red curve representing a ( X = e R e R , Y = . R + 1 . parameterization. Thus, converting the density of the ( R, T ) plane involves a Jacobian term, and the observed posterior is warped and non-Gaussian.Still, slice sampling can be performed for drawing a new R with T = 0 fixed, where in firsta random height of density is chosen underneath the original red dot starting point. Thena uniform variable is drawn from the band of values for R who have posterior density whohave equal or higher density than the selected height. We note that not all Y ij are observed, sometimes because SNPs do not differenti-ate between crosses, and in other cases because of random experimental missingness.Imputation of the data missing at random is naturally accommodated in a Bayesiansetting. After every iteration, the MCMC sampler draws new samples Y miss ij missingdata based upon current estimates of P i and S j , R j , η g,j . Using these new drawsof Y miss to make the data complete, samples of P i , S j , R j , η g,j | Y complete are drawn.When a gene is entirely excluded for a certain cross because of identical alleles, wedo not impute Y miss ij for that cross, and treat η g,j as essentially non-existent. eta model for X-inactivation Because we vectorize and separate certain operations, computational costs areproportional to O ( N J ) to update all (cid:1) P , (cid:0) S parameters, but O ( N J ) to update all (cid:0) R parameters and O ( JGN ) to update all η parameters, as well as O ( G ) to update (cid:1) µ, (cid:1) α parameters. Thus estimation of (cid:0) R seems to be the most inefficient step inthe computation, where most of the computation costs come in calculation of betavalues Γ( R j P i S j e ηg,j +1)Γ((1 − R j )(1 − P i ) S j e ηg,j +1)Γ((1 − R j − P j +2 R j P j ) S j e ηg,j +1) , which do not allow for the efficientseparation of parameters.Gibbs samples can be used to infer posterior mean, standard-deviation, andquantiles which we use as evidence to simultaneously test biologically relevant hy-potheses, such as H µ : µ g = . , H S : S j = S j , and H η : (cid:80) g ( η g,j − ¯ η j ) < . Inparticular, tests for hypotheses of the form H µ ,µ : µ g = µ g are desired, and wehope 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 (Gelmanand Rubin, 1992) convergence measures reported in Section 5 using the completedata are quite low, so we conclude that the chain is optimized sufficiently for thepurpose of testing our biological hypotheses. 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-biasedcoin”. Consider a cross between mother carrying
Xce allele “ a ” and the fatherwith allele “ b ”. One could imagine putting a coin together with a mass M a forthe heads face, and a mass M b for the tails face. If a physical model suggested µ ( a × b ) = M a / ( M a + M b ) , we could measure M c in a third allele and predict µ ( a × 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 copiesattract nucleosomes and transcription factors to attach to the X -chromosome anddeactivate it, suggesting that longer sequences at Xce are weaker alleles.Now that we have declared multiple inbred strains of mice to carry
Xce a , Xce b , Xce c , . . . alleles, we can implement this model within our framework. Denote theprevious specification, where µ g are given i.i.d. priors, our “Independence Assump-tion” (IA) model. After analysis, we can compare IA estimates, where no structurewas assumed for µ 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 µ g as a function of additive allele parameters a , a , . . . , a K a for K a total alleles andparent-of-origin allele-specific parameters m , m , . . . , m K m for a total number ofparent-of-origin alleles K m . In our case we propose:logit ( µ g ) = ( a d ( g ) + m d ( g ) ) − ( a s ( g ) − m s ( g ) ) . (11)Having a larger a d ( g ) dam allele than the a s ( g ) sire allele pulls µ g in the direction µ g > . . Parent-of-origin effects m d ( g ) perturb this result based upon whetherstrains serve as dam or sire.Again, we impose constraints (cid:80) K a k =1 a k = 0 and (cid:80) K m k =1 m k = 0 . Vector-constrained slice sampling then allows for posterior sampling of these model pa- Alan Lenarcic et. al
Table 5.
Simulating 1000 datasets from a fixed (cid:126) (cid:1) µ vector in model True µ ˆ µ − µ -0.003 -0.001 -0.001 0.001 0Method Mean (cid:112) (ˆ µ − µ ) % HPD Width 0.056 0.063 0.063 0.061 0.056 % HPD Coverage 0.974 0.972 0.962 0.962 0.969Sample Mean (Bias) ˆ µ − µ (cid:112) (ˆ µ − µ ) % CI Width 0.054 0.061 0.06 0.056 0.051 % CI Coverage 0.619 0.957 0.933 0.786 0.487 rameters, along with the previously mentioned η , (cid:0) R .
4. Verifying Method Performance through Simulation
Gibbs sampler models are computationally slow, and one can rarely explore allregimes to test a Bayesian estimator. It would be wrong to draw parameters fromprior distributions to simulate datasets; to demonstrate frequency performance ofposterior estimates we must use fixed parameter values. We simulate datasets ofsize n = 200 individuals, J = 6 tissue gene measurements, making G = 5 groupswith samples each. We fix (cid:1) µ truth = ( . , . , . , . , . T , which is an ordered setof realistic values for the group means, where adjacent µ g values are separated by .2,.05, .15, and .1. We choose tissue-gene precisions (cid:0) S truth = (200 , , , , , and set (cid:0) R truth such that logit ( (cid:0) R truth ) = ( − , − , − , , , / . Note that the mostprecise gene-expression values will be downward biased. We will keep the truesecondary precision matrix η truth = , however, the Bayesian estimator will still fit η as a free × matrix. We set (cid:1) α truth = (50 , , , , T . We draw a P i forevery individual i based upon group membership g ( i ) , from which Y ij values aresampled, before 66% of Y ij are removed to represent data missingness (requiring atleast one measurement per individual).Our estimation draws 2000 Gibbs samples from the Bayesian Posterior, withfit time approximately 10 minutes, and disk usage approximately 100mb. We es-timate posterior mean and credibility intervals for all model parameters. We ransimulations on the UNC KillDevil Cluster, which allowed roughly 200 cores in par-allel, fitting 1000 replications of all of our simulation experiments (including thosedescribed later in this section) in roughly a day. We compare Bayesian estimatesfor µ g to the sample mean of groups g : ¯ Y g = (cid:107) Y obs ij (cid:107) (cid:80) g ( i ) = g , j observed Y obs ij , where,confidence intervals are based upon the t -distribution using standard deviation es-timated from the data.Simulation results in Table 4, show that sample means ¯ Y g do not estimate thevector µ well, but that Bayesian credibility intervals have near-95% coverage oftrue µ g values. For this n = 200 dataset, with these values of (cid:0) S and (cid:126) (cid:1) α , Bayesian95% credibility intervals will be of approximate width 6%. This credibility width issufficient to differentiate almost all of the group µ ( t ) g . Only estimates for µ = . eta model for X-inactivation Table 6.
Simulation performance in estimating population simulatedmean.
True ¯ Y pop g ˆ¯ Y g − ¯ Y pop g (cid:113) ˆ¯ Y g − ¯ Y pop g ) % HPD Width 0.049 0.055 0.055 0.053 0.05 % HPD Coverage 0.974 0.972 0.957 0.962 0.976Sample Mean ˆ¯ Y g − ¯ Y pop g (cid:113) ( ˆ¯ Y g − ¯ Y pop g ) % CI Width 0.047 0.052 0.052 0.047 0.042 % CI Coverage 0.914 0.936 0.894 0.9 0.899Bootstrap Mean ˆ¯ Y g − ¯ Y pop g -0.001 0 0 0.001 0Mean (cid:113) ( ˆ¯ Y g − ¯ Y pop g ) % BI Width 0.054 0.061 0.06 0.056 0.051 % BI Coverage 0.939 0.955 0.936 0.937 0.951 and µ = . had overlapping posteriors. We use a rejection criterion × min ( P ( µ >µ , µ > µ )) < . . That is, we identify two µ j parameters as different if the two-sided tail posterior probability that their order is reversed is less than . . For thecomparison of µ , µ , this test has a power of 60.4%. Comparisons of group meansfor µ g − µ g (cid:48) ≥ . have near 100% power.Arguably, the sample mean, ¯ Y g , is not meant to estimate µ g , which is a quan-tity of the model. For a fairer comparison, consider the target of a populationgene-expression mean ¯ Y pop g of group g . This would be the population average geneexpression of the six genes for all possible mice of cross g . If we simulate 100Kindividuals or more from the model, we can calculate a satisfactory approximation ¯ Y pop g = J × K (cid:80) i (cid:48) Y i (cid:48) j . The sample mean ¯ Y g of 40 pups should bean estimate for ¯ Y pop g . In Table 4, we see that, in estimation of observable ¯ Y pop g thatthe Bayesian % credibility intervals tend to be conservative and over-cover ¯ Y pop g .The sample mean’s nominal % confidence intervals offer only % coverage for µ , µ , µ . The bootstrap has wider bootstrap intervals than the Bayesian credi-bility intervals, but have less coverage. The Bayesian credibility intervals are .01narrower, but have over-conservative % coverage of ¯ Y pop g .From the Bayesian model, we can also generate an estimate for ¯ Y pop g by simulat-ing data from the posterior predictive distribution for future Y i (cid:48) j . For each step ( t ) in the Gibbs Sampler, we can simulate 100K mice from the model using the currentstate of ( (cid:1) µ ( t ) , (cid:1) α ( t ) ) and calculate expected gene expressions using ( (cid:0) S ( t ) , (cid:0) R ( t ) , η ( t ) ) ,taking the average ¯ Y sim ( t ) g . These samples, ¯ Y sim ( t ) g , treated as draws from the pos-terior distribution for ¯ Y pop g .We compare estimates for ¯ Y pop using the Bayesian posterior predictive distri-bution, the sample mean, as well as the bootstrap mean including 95% bootstrapintervals, in Table 4.Table 4 presents average performance of the Bayesian model in estimating thefixed (cid:1) µ truth , (cid:1) α truth , (cid:0) S truth , (cid:0) S truth , η truth as well as fitting randomly simulated draws Alan Lenarcic et. al
Table 7.
Bayesian estimator performance measuring other modelparameters (cid:126)α (cid:126)µ P (cid:126)S (cid:126)R η √ Average Squared Bias ˆ θ .
994 0 .
001 0 .
002 31 .
091 0 .
001 0 . √ MSE for ˆ θ .
383 0 .
022 0 .
046 52 .
129 0 .
019 0 . Mean HPD Width .
007 0 .
06 0 .
157 189 .
364 0 .
039 2 . Mean HPD Coverage .
882 0 .
968 0 .
931 0 .
969 0 .
927 0 . for (cid:1) P . Parameter (cid:1) α tends to be under-covered, as it is a hierarchy parameter meantto describe variation in (cid:1) P , which are unobserved latent data in the model. On thissimulation, where (cid:126)µ was fixed, the hyper parameters hyper parameters such as τ η , α all are undefined. Arguably, the beta-distribution model may not reflect the true mechanism gener-ating the data. We simulate data an alternate linear mixed-effects model on thelogit scale: logit ( Y ij ) = µ g ( i ) + δ i + Bias j + N (0 , σ j ) . (12)In the model 12, observed Y ij are Gaussian on the logit scale, with a parameter µ g torepresent cross means, δ i is the i th individual’s average deviance from the populationmean, and each tissue-gene j supplies a different sized bias and measurement noise.We use Bias j = log e R truth j − R truth j using values of R truth j from the simulation above, anduse σ j = 1 /S truth j . The δ i are simulated from a N (0 , /α truth g = 1 / 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 theBayesian model’s estimate for ¯ Y pop g . In Table 4.1, we compare, as in Table 4,performance of the three estimators in this case where the model assumption forthe Bayesian method is incorrect. The Bayesian method has the shortest intervalsand yet is conservative, while the sample mean and bootstrap confidence intervalscan under-cover. We simulate data from the “weight-biased coin” model (WBC) to show our per-formance in measuring (cid:1) a . We fix 6 alleles with additive strengths (cid:1) a truth = {− , − , , , , } T / , but there will be no parent-of-origin effect, (cid:1) m truth = { , , , , , } T . We produce n = 320 individuals per simulation, which will bein groups of size from a random selection of 8 crosses from the × − possible crosses. We simulate a study design that guarantees each allele is sampledat least once, but do not guarantee that an allele is featured both as a dam anda sire. In the real experiment, allele membership was not known or anticipateduntil the crosses were performed. We use simulation to judge whether a randomselection of crosses of this type can properly differentiate α k , α k (cid:48) . We use the same eta model for X-inactivation Table 8.
Simulation estimating ˆ¯ Y pop g from alternate linear mixed-effects True ˆ¯ Y pop g ˆ¯ Y g − ¯ Y pop g (cid:113) ( ˆ¯ Y g − ¯ Y pop g ) % HPD Width 0.025 0.031 0.031 0.029 0.025 % HPD Coverage 0.99 0.985 0.982 0.986 0.986Sample Mean ˆ¯ Y g − ¯ Y pop g (cid:113) ( ˆ¯ Y g − ¯ Y pop g ) % CI Width 0.031 0.038 0.038 0.034 0.028 % CI Coverage 0.947 0.956 0.946 0.951 0.948Bootstrap Mean ˆ¯ Y g − ¯ Y pop g (cid:113) ( ˆ¯ Y g − ¯ Y pop g ) % BI Width 0.033 0.041 0.042 0.037 0.03 % BI Coverage 0.918 0.935 0.939 0.945 0.945 Table 9.
Performance of ˆ a k estimates in Bayesian WBC model True a k -0.625 -0.125 0 0.125 0.25 0.375Bayes Mean ˆ a k − a k -0.003 0.006 0 -0.001 0 -0.043Mean (cid:112) (ˆ a k − a k ) % HPD Width 0.624 0.608 0.611 0.609 0.602 0.616 % HPD Coverage 0.998 0.998 0.999 1 0.999 0.997 (cid:0) R truth , (cid:0) S truth , η truth , (cid:1) α truth from before and verify that confidence intervals coverwith desired accuracy.Table 4.2 shows that, while coverage is a proper 95%, the confidence intervals arequite wide, approximately . or ± . for every allele effect. This is disheartening;adjacent allele effects are rarely distinguished in the posterior. HPDs for (cid:1) a areaffected by uncertainty for m k reciprocal effects, which is large when a strain hasnot served at least once as both dam and sire. For this reason, we retest in Table 4.2in the same simulation framework, but where (cid:1) m is set to zero. This producesnarrower HPDs by 2/3 allowing adjacent a j , a j (cid:48) to be distinguished.Parent-of-origin effects do, however, appear to be present in the real data, in-dicating nonzero (cid:1) m . We conclude that a study would need more than 8 crosses tostudy these alleles and demonstrate in Table 4.2 that 16 of 30 crosses is sufficient.In the real dataset we collected data from 12 out of 20 allele combinations Weused our simulations to encourage experimenters to perform additional reciprocalcrosses. Table 10.
Repeat of Table 4.2 where (cid:1) m = 0 is assumed True a k -0.625 -0.125 0 0.125 0.25 0.375Bayes Mean ˆ a k − a k (cid:112) (ˆ a k − a k ) % HPD Width 0.176 0.172 0.171 0.17 0.173 0.172 % HPD Coverage 0.965 0.964 0.956 0.961 0.96 0.964 Alan Lenarcic et. al
Table 11.
Simulation where (cid:1) m is estimated featuring 16 randomlysampled crosses True a k -0.625 -0.125 0 0.125 0.25 0.375Bayes Mean ˆ a k − a k -0.002 0.001 0.002 -0.001 0 0.001Mean (cid:112) (ˆ a k − a k ) % HPD Width 0.124 0.119 0.119 0.118 0.119 0.121 % HPD Coverage 0.967 0.965 0.954 0.959 0.971 0.968
5. Data Analysis
The Calaway et al. (2013) dataset comprises 660 mice, with 24 tissue-gene combi-nations, for brain, liver, and kidney. Missingness is such that for 15840 possiblemeasurements, only 2393 could be observed; thus 84.9% of potential Y is unob-served. In addition, 34 mice, with 10 brain-derived RNA-seq gene measurementsare provided from a set of WSB , PWK , CAST wild-derived F crosses. These crosseswere measured on a different set of 10 genes and no pyrosequencing measurementswere taken. Using both our IA and WBC models for (cid:1) µ , we implement Gibbs sam-plers using the observed Y to draw from the posterior of (cid:1) P , (cid:1) µ , and (cid:0) S . In Figure ?? ,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 µ g values, fitted values, ˆ P i , are different from the unweightedarithmetic mean, (cid:80) j Y ij /J i , of the observed gene expressions. Gene measurementswith 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 courseof a few hours, with i.i.d. Beta (5 , -distributed starting values, the maximum95% quantile Gelman-Rubin convergence diagnostic (GRD) is 1.11 for the 32 (cid:126)µ parameters with a multivariate potential scale reduction factor (PSRF) of 1.08 anda 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 (cid:0) S parameterswith Gamma (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 mixingof the chains is sufficient for analysis of this dataset.The most important scientific objectives are accomplished by estimation of (cid:1) µ g .Most strains were anticipated to carry Xce b , and so we used three pieces of evi-dence to converge this suspicion. We established these candidates carried allelesstronger than Xce a by showing that µ g > . when these candidates served as damwhen crossed with known Xce a carriers, we established these candidates carriedalleles weaker than Xce c by showing that µ g < . when they served as dam whencrossed with Xce c carriers, and finally we tried to establish that µ g was close to . when these candidates were crossed with known Xce b carriers. For crosses thesecandidates served as sire, the ordering of µ g is reversed. We rejected undeterminedcandidates as allele members when crosses with preexisting allele holders demon-strated µ g (cid:54) = . . We used a criterion of × min ( P ( µ g > . | Y ) , P ( µ g < . | Y )) , thatis, twice the smallest tail posterior probability on either side of . . An example ofposterior estimation for µ g in crosses with ALS is given in Figure 3, and a summarytable of crosses of interest is in Table 5. These give estimates for µ g and confidence eta model for X-inactivation measures, for crosses of the unknown strains with known “ Xce a ”, “ Xce b ”, “ Xce c ”carriers. When multiple crosses between two alleles have been conducted, we re-port p g = 2 min ( P ( µ g > . , P ( µ g < . for the lowest posterior value among thosecrosses.When multiple crosses, including reciprocals, against an allele carrier were per-formed, we present the estimates with the smallest p g . A “*” suggests that wehave made a scientific decision to associate this strain with this allele. **: ALS was tentatively called a “
Xce b ” since the ALS × B6 cross had µ g = 0 .
46 (0 . , . with p g = 0 . for animals, despite the reciprocal cross showing such a strongparent-of-origin effect (also for animals). ALS/LtJ
A/JC57BL/6JC57BL/6JCAST/EiJ
ALS* x CASTALS* x B6B6 x ALS*ALS* x AJ
ALS* x CAST
Ddx26b b b bbb bbb bbbb bbbb bb bbbb bbbb bb bbbk k kkk kkk kkk k kk kk kk kk kk k kkkv v vvv vvv vv vvv vv vv v v v vvv vv vvv
Fgd1 b b bbb bbb bbb bbb bb bb bb bbb bbb bb bbbk kkk kkk kk k k kkkk kk kk kk k kkkv v vvv vvv vvvv v vv vv v v v vv v vvv vv
ALS* x B6
Fgd1 bb bbb bb bbb bbbb bbbbb bbb bbbb bb bk kk k kkk k k kkk kkkk kk k k kk k kk k vv vv vvv v vv vvv vvvv vvv vvvv
Mid2 bb bbb bb b bbbb bb bb bb bb
B6 x ALS*
Fgd1 bb bb b bb bbbbbb b bbb bbbb bb bbb bbbb
Mid2 bb bb b bb bbb bbbb bbb bbb bbb bb b bbbb
Ddx26b kk kk k kk kkkkkkk kkk kkk k kk kk k kkkk
ALS* x AJ
Ddx26b b bbb bb bb bb bb bb bb b bb bb bbb bb bb bbk kk k kk k kkkkkkk k k kk kk kk kkk kk kkv vv v vvvvv vvv vv vvvv v vvv vv vv vv
Fgd1 b bbb bb bb bb bb bb bbb bb bb bbb bb bb bb k kk k kk kk kkk kkkk k kk kk kk kkk kk kkvvvv v vv vvvvv vv vv vv v vvv vv vv vv
Rragb b b bb b bbb b bbb bbb bb bbbbbb kkk kkk k kk kkkkkk k kk kk kk kkk kkkk v vvv vvvv vv vv vv vv vvv vvv
Fig. 3.
Left: for four crosses using strain
ALS , we plot posterior density for µ g in terms ofactive ALS . Blue represents crosses with known
Xce b carriers (dotted line for the reciprocalcross where B6 served as maternal strain), black for Xce a carriers, and red for Xce c . Center:we plot posterior 95% credibility intervals for µ g of these crosses. The posterior suggestsit is highly unlikely that µ g for the ALS × CAST cross is greater than .5, with mean .34. Theposterior for
ALS × AJ cross is concentrated above . , with mean .59. The ALS × B6 crosscarries posterior weight over µ g = . . The reciprocal B6 × ALS cross does not, suggestingsignificant parent-of-origin bias. On the right are observed gene expressions Y ij used toestimate µ g . Points are labeled “b”, “k”, “v” for brain, kidney and liver respectively. Operating on this principle we confirmed that
LEWES , PERA , TIRANO , SJL , WLA , WSB , and
ZALENDE carry
Xce b . Tentatively, we believe ALS to be a
Xce b carrier, sinceit is not of type Xce a or Xce c . While B6 × ALS cross showed bias from . , the ALS × B6 seemed much closer to a cross between identical carriers. PWK appears to havean allele with strength near that of the
Xce b allele, although previous literaturehypothesized that it contains a Xce c as does CAST . Table 5.1 displays logit-scale estimated allele effects based upon our decisions forcarriers of
Xce a , Xce b , Xce c , as well as yet-unidentified placeholders for PWK ( Xce p ),and PANCEVO ( Xce v ). We observed that parent-of-origin effect is not necessarilyattached to Xce allele, so we separated strains into potential parent-of-origin carriersbased 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 B6 ,and m WSB for
WSB is large and negative. Wild-derived mus musculus representative Alan Lenarcic et. al
Table 12.
Estimate ˆ µ g , tail probability ( p g ) against the hypothesis H : µ g = . . Xce a Xce b Xce c µ g p g µ g p g µ g p g ALS . . , .
63) 0 .
003 0 . . , . ** .
003 0 . . , .
43) 0 . SJL . . , .
61) 0 .
003 0 . . , . * .
142 0 . . , . < . LEWES . . , .
64) 0 .
003 0 . . , . * . WLA . . , .
59) 0 .
009 0 . . , . * . WSB . . , .
67) 0 .
022 0 . . , . * .
545 0 . . , .
33) 0 . PWK . . , .
76) 0 .
008 0 . . , .
67) 0 .
007 0 . . , .
28) 0 . TIRANO . . , . < .
001 0 . . , . * . ZALENDE . . , . < .
001 0 . . , . * . PERA . . , . * . PANCEVO . . , .
5) 0 .
082 0 . . , .
44) 0 . Table 13.
Logit-scale estimated allele effects
Allele Additive (cid:1) a Xce a -0.36(-0.52, -0.18) Xce b -0.01(-0.19, 0.19) Xce c Xce p -0.23(-0.42, -0.03) Xce v -0.67(-0.88, -0.45) , Maternal (cid:1) m AJ / ALS -0.07 (-0.5, 0.32) B6 CAST
LEWES
PWK
WSB -0.81 (-1.54, -0.08) eta model for X-inactivation PWK appears to carry an additive allele with effect strength between
Xce a and Xce b .Results for the weight-biased coin (WBC) model were not considered essential tothe Calaway et al. (2013) paper, but this model can be incorporated with geneticsequence experiments to search for parent-of-origin loci. Comparing posterior means ˆ µ g from the WBC and the original IA model for µ g ,the mean (cid:107) ˆ µ WBC g − ˆ µ IA g (cid:107) is .12 and the correlation Corr (ˆ µ WBC g , ˆ µ IA g ) is . . The cross“HF” which is a WSB × CAST cross measured using RNA-seq shows the maximumdifference with ˆ µ IAg = 0 .
42 ( . , . and ˆ µ W BCg = 0 . . , . .Figure 4 demonstrates the differences between WBC and IA models. Pointssizes 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 sacrificefit for the smaller crosses. For the WBC model, strains used only as dam or sirehave fixed m k = 0 . The green points in Figure 4 are between strains which bothhave non-zero ˆ m k , whereas the red points include one ˆ m k contribution, and thepurple points have none.HPD intervals in the WBC and IA models do not overlap for only 6 crosses:HF ( WSB × CAST ), AG ( AJ × PWK ), BZ ( B6 × ZALENDE ), SB (
SJL × B6 ), BT ( B6 × TIRANO ) and BO ( B6 × PANCEVO ). Additionally, in crosses GH (
PWK × WSB ), HG(
WSB × PWK ), HA (
WSB × AJ ), GA ( PWK × AJ ), FH ( CAST × WSB ) the WBC modelseems to underfit the observed result, though these crosses are smaller and HPDintervals overlap.Wide disparities between ˆ µ WBC and ˆ µ IA may suggest crosses have been insuffi-ciently 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 suc-cessfully predicts µ g . Tables 5.3 and 5.3 presents estimated (cid:0) S precision parameters for the IA model.Some RNA-seq gene measurements (taken only in the brain), have similar S j levelsas pyrosequencing. Ideally, we would have liked to measure both pyrosequencingand RNA-seq on the same tissue, same mouse, same cross same gene. But we wereonly able to procure RNA-seq from a set of PWK - WSB - CAST crosses sequenced to servemultiple 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. Estimatesfor ˆ S j range from 300 (suggesting a s.d. of about ± . near P i = . ) to someRNA-seq measurements that only display ˆ S j = 8 (an s.d. of about ± . ) accuracy.Based upon the Pòlya Urn abstraction described earlier, we can interpretthese 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 genesrange from − to , meaning that these genes were not measured in enoughcrosses provide reliable contribution. Alan Lenarcic et. al . . . . WBC vs. IA, Posterior Means
Independence Assumption, m ^ gIA W e i gh t − B i a s ed C o i n , m ^ g W B C l l l l l l l l l ll l l l l l AG GAHASB FHHF GH HGBTBZ BO ll l l
10 individuals40 individuals
Fig. 4.
Comparison of posterior means ˆ µ IA on x-axis and ˆ µ WBC on the y-axis. Size of pointrepresents the number of individuals used in a cross. Strains not used as both dam andsire could not be measured for a maternal effect (cid:126)m . Green points correspond to crossesbetween two strains that each contribute a m j effect. Red points included one strain whichhad fixed m j = 0 because it is unidentifiable. Purple points have both strains fixed m j = 0 .We plot confidence ellipses for the crosses that do not agree between WBC and IA. Thereare 6 crosses out of 34 where % HPD intervals between the two models do not overlap. Table 14. (cid:0) S precisions from pyrosequencing (brain/kidney/liver) measurements Fgd1 Ddx26b Rragb1 Mid2 brain . . . ( . , . ) ( . , . ) ( . , . ) ( . , . )kidney . . . . ( . , . ) ( . , . ) ( . , . ) ( . , . )liver . . . . ( . , ) ( . , . ) ( , . ) ( . , . ) Zfp182 Xist Hprt1 Acs14 brain . . . . ( . , . ) ( . × − , ) ( . × − , . ) ( . × − , . )kidney . . . . ( . , ) ( . × − , . ) ( . × − , . ) ( . × − , . )liver
349 28 . . . ( . , . ) ( . × − , . ) ( . × − , . ) ( . × − , . ) eta model for X-inactivation Table 15.
Precision in genes from RNA-seq Measurements (brainonly)
Gripap1 Zdhhc9 Ids Pls3 Arhgef9 . . . . . ( . , . ) ( . , . ) ( . , . ) ( . , . ) ( . , ) Zmym3 Sh3bgr1 Gprasp1 Gnl3l Tspyl5 . . . . . ( . , . ) ( . , . ) ( . , ) ( . , . ) ( . , )
6. Discussion and Conclusions
We have developed a Bayesian hierarchical model, robust and powerful enoughto confirm our hypotheses, which uses allele-specific gene-expression to estimatewhole-body X -inactivation. In simulations, we showed that Bayesian 95% HPDcredibility intervals reliably cover the truth in simulated datasets similar in sizeand structure to our experimental dataset, and that these intervals improve cov-erage relative to asymptotic- t and bootstrap-based estimators, even when data isgenerated from an alternate model. Our “independence-assumption” (IA) model,where µ g ’s are given i.i.d. priors, demonstrated and confirmed allele membershipof 10 previously untested strains, which was the critical analysis of the Calawayet al. (2013) experiment. Our weight-biased coin (WBC) model, where µ g is afunction of parental strains, including additive allele effects (cid:1) a and reciprocal effects (cid:1) m , produced similar results to the IA prior, suggesting this scientific model mightbe used to predict µ g in future crosses.Modeling bias and precision of gene-specific measurements using free ( (cid:0) S )and constrained ( (cid:0) R, η ) parameters, we developed techniques in constrained slice-sampling that required both linear and non-linear transformation of the parameterset. Our model could accommodate allele-specific expression measured both bypyrosequencing and RNA-seq.Here we must acknowledge limitations in experimental design that led toBayesian recourse. Quantitative trait loci (QTL) are regions within the genomewhose variation is highly correlated, and conjectured to be causal, with a pheno-type — in our case X -inactivation skewing. Typical experiments to locate QTL,such as Advanced Intercrosses (Darvasi and Soller, 1995), in which individuals ofknown genotype are mated, rely on random recombination of the genome. But pre-cise 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 anddifference within the Chadwick et al. (2006) interval. Establishing that two strainshave heavy difference in a sub-region, yet their offspring exhibit no X -inactivationskewing, gives credence to a belief that the sub-region of difference is not a candi-date for Xce . Establishing that two strains, identical at a sub-region, have offspringthat exhibit X -inactivation skewing, rejects a sub-region in favor of other regionswhere these two strains do differentiate genetically.Prior to experiment we did not know all of which crosses to perform, did notknow which sub-region we wished to justify, and did not know which strains would Alan Lenarcic et. al differentiate. As early strains in the experiment were established as carriers forknown alleles “ a ”, “ b ”, “ c ”, we were further motivated to test other strains thatcontinued to break down the candidate region. Such experimental choices of newtest strains were done both in the light and in the dark of statistical analysis andestimates; as such, our stopping rule cannot be claimed to be a sample-independentchoice. The “null” sampling space in this case is ill-posed: we did not randomlysample from a population, but instead chose sub-populations to explore based uponresearch intuition. Sometimes samples were obtained because they were readilyavailable or cheaper to obtain. But, at other times, crosses were ordered becausethey were hoped to lead to a new allele. Sometimes additional samples were ob-tained to improve statistical confidence on a cross, adding to samples already col-lected. From a frequentist perspective, it is right to be skeptical of analysis onsuch a dataset. Since we cannot map out the decision metric that led to collectionof samples, we cannot postulate the alternative sample space for most hypothesistests, and cannot report a p -value.Bayesian exchangeability (de Finneti, 1937; Bernardo and Val, 1996; Lindleyand Phillips, 1976; Hewitt et al., 1955; Jackman, 2009; Diaconis and Freedman,1980) often justifies posterior analysis of sequential designs, and here we show thesampling mechanism is an ignorable design (Thompson and Seber, 1996). Con-sider indicators I i, , I i, , . . . , I i,G which represent whether individual i is in group g ∈ { , , . . . , G } , such that (cid:80) g I i,g = 1 . Let (cid:0) I i vector be assigned with a tran-sition probability f i trans ( (cid:0) I i | Y < i , I < i ) , which is a multinomial dependent only uponprevious observations and sampling choices. Conditional on (cid:0) I i , (cid:1) µ, (cid:1) α , the P i are in-dependent from other individuals. The full posterior, including arbitrary sequentialdesign, is P ( Θ | Y ) ∝ n (cid:89) i =1 f i trans ( (cid:0) I i | Y < i , I < i ) × Beta ( P i ; (cid:88) g I i,g µ g α g + 1 , (cid:88) g I i,g (1 − µ g ) α g + 1) × (cid:89) j Beta ( Y ij ; P i R j S j e η g,j + 1 , (1 − P i )(1 − R j ) S j e η g,j + 1) × Priors ( S , R , µ , α , η ) . (13)We see from expanded equation 13 above that because of conditional inde-pendence given P i , and because all vectors (cid:0) I i are observed and known, then thetransition density, f i trans ( (cid:0) I i | Y < i , I < i ) , which is only based upon observables, cannotinfluence the posterior. Because experimenters had no preference before observa-tions of allele identities for the individuals chosen, we can accept that the stoppingrule is proper (though certainly not optimal) in the definition of a proper stoppingrule given in Parmigiani and Inoue (2009). We would have eventually stoppedcollecting samples, no matter the true value of µ g .We note, per a criticism in Lindley and Smith (1972), that our (cid:1) µ parametersare not actually exchangeable based upon science. Even if we do not know thevalue of µ g , we do know that if µ g , µ g (cid:48) , µ g (cid:48)(cid:48) were a sequence of crosses of a mother eta model for X-inactivation strain with unknown allele “UNK” to fathers of allele carriers “a”, “b”, “c”, then thereshould be a decreasing sequence: µ g ≥ µ g (cid:48) ≥ µ g (cid:48)(cid:48) . In the “IA” model, we did notmodel allele-order hypothesis with our priors and relied upon the data to reveal thissequence. The WBC model did enforce this relation, and reached similar estimatesof ˆ µ g .It is difficult to statistically argue for the hypothesis µ g (cid:48) = µ g , in contrastto proving µ g (cid:48) (cid:54) = µ g . We have relied upon a criterion of overlapping posterior.It is possible that the difference between two cross means, µ g − µ g (cid:48) = ∆ , mightbe small, non-zero, but imperceptible, such as ∆ ≈ . . In which case, wewould inevitably decide incorrectly to treat µ g , µ g (cid:48) as fundamentally equivalentcrosses between same alleles. We assume that alleles differ in effect enough to bestatistically observable.As high-throughput sequencing becomes cheaper, Bayesian algorithms of ourdesign must be rejected in favor of more efficient methods. The Bayesian algorithmhas complexity O ( N J ) , due to Gibbs sampling of individual S j and R j , makingit unsuitable for whole-transcriptome datasets, where a complete genome has J ≈ . 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 ≈ $20. A fully-analyzed RNA-seq sample might measure an individual’s P i averagebrain expression to a ± . region of 95% credibility, but pyrosequencing of threegenes would leave a ± region. But to measure cross mean µ g , it is better to have P i measured imprecisely in many individuals in multiple tissues rather than havehigh precision on any single-tissue P i .We note that point-estimate performances of the sample mean estimator andbootstrap are nearly identical to the posterior mean. Were it not for under-coveragein their confidence intervals, and a conceptual disconnect between these proceduresand our parametric model, these methods might have served as suitable replace-ments to the Bayesian Gibbs sampling.To the extent we have introduced techniques more generally applicable to statis-tics, our constrained priors show promise in other naturally unidentifiable problemswith small order J . We are fortunate to be in a data-setting where a Bayesian es-timator is statistically robust and viable, and can use posterior sampling to test aparametric model for its comparable nuances.
7. Acknowledgements
This project was supported by National Institutes of Health (NIH) grantsR01GM104125 (ABL, WV), R35GM127000 (WV), and P50MH090338 andP50HG006582 (FPMdV, JDC). We also thank UNC Information Technology Ser-vices for computational support. The funders had no role in study design, datacollection and analysis, decision to publish, or preparation of this manuscript Alan Lenarcic et. al
References
Barr, M. L. and E. G. Bertram (1949). A Morphological Distinction between Neu-rones of the Male and Female, and the Behaviour of the Nucleolar Satellite duringAccelerated Nucleoprotein Synthesis.
Nature 163 , 676–677.Bernardo, M. and U. D. Val (1996). The Concept of Exchangeability and its Ap-plications.
Far East Journal of Mathematical Science 4 , 111–121.Brooks, S. and A. Gelman (1998). General methods for monitoring convergenceof iterative simulations.
Journal of Computational and Graphical Statistics 7 ,434–455.Calaway, J., A. Lenarcic, J. P. Didion, J. R. Wang, J. B. Searle, L. McMillan,W. Valdar, and F. de Pardo (2013). Genetic architecture of skewed x inactivationin the laboratory mouse.
PLoS Genetics 9 (10).Cattanach, B. and J. Isaacson (1967). Controlling elements in the mouse X-chromosome.
Genetics 57 (2), 331–346.Cattanach, B. and C. Rasberry (1991). Identification of the Mus spretus Xce allele.
Mouse Genome 89 , 565–566.Cattanach, B. and C. Rasberry (1994). Identification of the Mus castaneus Xceallele.
Mouse Genome 92 , 114.Cattanach, B. M. and J. H. Isaacson (1965). Genetic control over the inactivation ofautosomal genes attached to the X-chromosome.
Molecular and General GeneticsMGG 96 (4), 313–323.Chadwick, L. H., L. M. Pertz, K. W. Broman, M. S. Bartolomei, and H. F. Willard(2006). Genetic Control of X Chromosome Inactivation in Mice: DeïňĄnition ofthe Xce Candidate Interval.
Genetics 173 (4), 2103–2110.Chang, T. (1983). Binding of cells to matrixes of distinct antibodies coated on solidsurface.
J Immunol Methods. 65 , 17–23.Darvasi, A. and M. Soller (1995). Advanced Intercross Lines, an ExperimentalPopulation for Fine Genetic Mapping. (3), 1199–1207.de Finneti, B. (1937). Foresight: Its Logical Laws, Its Subjective Sources.
Annalesde l’Institut Henri Poincare 7 .Developers, R., R. A. Becker, Chambers, J. M., and A. R. Wilks. R documentation:Build shared object/dll for dynamic loading. http://stat.ethz.ch/R-manual/R-devel/library/utils/html/SHLIB.html . Accessed: 2013-09-13.Diaconis, P. and D. Freedman (1980). Finite exchangeable sequences.
The Annalsof Probability 8 (4), 745–764.Efron, B. and C. Morris (1977). Stein’s Paradox in Statistics.
Scientific Ameri-can 236 , 119 – 127. eta model for X-inactivation Ferrari, S. and F. Cribari-Neto (2004, Aug). Beta Regression for Modelling Ratesand Proportions.
Journal of Applied Statistics 31 (7), 799–815.Gelfand, A. E., A. F. M. Smith, T.-m. Lee, A. F. M. Smith, and A. E. Gelfand(1992). Bayesian Analysis of Constrained Parameter and Truncated Data Prob-lems Using Gibbs Sampling.
Journal of the American Statistical Associa-tion 87 (418), 523–532.Gelman, A. and D. Rubin (1992). Inference from iterative simulation using multiplesequences.
Statistical Science 7 , 457–511.Gendrel, A.-V. and E. Heard (2011, Dec). Fifty years of X-inactivation research.
Development (Cambridge, England) 138 (23), 5049–55.Graze, R. M., L. L. Novelo, V. Amin, J. M. Fear, G. Casella, S. V. Nuzhdin, andL. M. McIntyre (2012, Jun). Allelic imbalance in Drosophila hybrid heads: exons,isoforms, and evolution.
Molecular biology and evolution 29 (6), 1521–32.Hall, D. A., J. Ptacek, and M. Snyder (2007). Protein Microarray Technology.
Mech Ageing Dev. 128 , 161–167.Hewitt, E., L. J. Savage, E. Hewitt, and L. J. Savage (1955). Symmetric Measures onCartesian Products.
Transactions of the American Mathematical Society 80 (2),470–501.Jackman, S. (2009).
Bayesian Analysis for the Social Sciences (1 ed.). Chichester,West Sussex, UK: John Wiley & Sons, Ltd.Kambere, M. B. and R. P. Lane (2009, Feb). Exceptional LINE density at V1Rloci: the Lyon repeat hypothesis revisited on autosomes.
Journal of molecularevolution 68 (2), 145–59.Kim, K., S. Thomas, I. B. Howard, T. A. Bell, H. E. Doherty, F. Ideraabdullah,D. A. Detwiler, and F. P.-m. Villena (2005). The genus Mus as a model forevolutionary studies Meiotic drive at the Om locus in wild-derived inbred mousestrains.
Biological Journal of the Linnean Society 84 , 487–492.Lawson, C. L., R. J. Hanson, D. Kincaid, and F. T. Krogh (1979). Basic linearalgebra subprograms for fortran usage.
ACM Trans. Math. Soft. (5), 308–323.Lee, J. (2011). Gracefully ageing at 50, X-chromosome inactivation becomes aparadigm for RNA and chromatin control.
Nature Reviews Molecular Cell Biol-ogy (12), 815–826.Lindley, D. and L. Phillips (1976). Inference For a Bernoulli Process ( a BayesianView )*.
The American Statistician 30 (3), 112–119.Lindley, D. V. and A. F. M. Smith (1972). Bayes Estimates for the Linear Model.
Journal of the Royal Statistical Society. Series B 34 (1), 1–41.Lyon, M. F. (1961). Gene Action in the X-chromosome of the Mouse (Mus musculusL.).
Nature 190 , 372 – 373. Alan Lenarcic et. al
Lyon, M. F. (1962, Jun). Sex chromatin and gene action in the mammalian X-chromosome.
American journal of human genetics 14 , 135–48.Mak, W., T. B. Nesterova, M. de Napoles, R. Appanah, S. Yamanaka, A. P. Otte,and N. Brockdorff (2004). Reactivation of the paternal x chromosome in earlymouse embryos.
Science 303 (5658), 666–669.Morin, R., M. Bainbridge, A. Fejes, M. Hirst, M. Kryzwinski, T. Pugh, H. McDon-ald, R. Varhol, S. Jones, and M. Marra (2008). Whole Transciptome ShotgunSequencing.
Biotechniques 45 (1), 81–94.Murphy, K. (2007). Software for Graphical models : a review.
Isba Bulletin De-cember , 1–3.Neal, R. (2003). Slice Sampling Author.
Annals of Applied Statistics 31 (3), 705–741.Nishioka, Y. (1995, Mar). The origin of common laboratory mice.
Genome /National Research Council Canada = Génome / Conseil national de recherchesCanada 38 (1), 1–7.Parmigiani, G. and L. Inoue (2009).
Decision Theory: Principles and Approaches (1 ed.). Wiley.Plummer, M. (2003). JAGS : A program for analysis of Bayesian graphical modelsusing Gibbs sampling JAGS : Just Another Gibbs Sampler.
DSC 2003 WorkingPapers .Roberts, G. O., J. S. Rosenthal, and S. Rosenthal (2012). Convergence of slicesampler Markov chains.
Journal of the Royal Statistical Society. Series B 61 (3),643–660.Ronaghi, M., M. Uhlén, and P. l. Nyrén (1998). A Sequencing Method Based onReal-Time Pyrophosphate.
Science 281 (5375), 363–365.Ruvinsky, A. (2009).
Genetics and Randomness . Boca Raton, FL: CRC Press.Sebat, J., B. Lakshmi, J. Troge, J. Alexander, J. Young, P. Lundin, S. Må nér,H. Massa, M. Walker, M. C. Chi, N. Navin, R. Lucito, J. Healy, J. Hicks, K. Ye,A. Reiner, T. C. Gilliam, B. Trask, N. Patterson, A. Zetterberg, and M. Wigler(2004). LARGE-SCALE COPY NUMBER POLYMORPHISM IN THE HUMANGENOME.
Science 305 (5683), 525–528.Shi, X., M. Styner, J. Lieberman, J. G. Ibrahim, W. Lin, and H. Zhu (2009, Jan).Intrinsic Regression Models for Manifold-Valued Data.
Journal of the AmericanStatistical Association 5762 , 192–199.Stankiewicz, P. and J. R. Lupski (2010). Structural Variation in the Human Genomeand its Role in Disease.
Annual Review of Medicine 61 , 437–455. eta model for X-inactivation Takagi, N., O. Sugawara, and M. Sasaki (1982, January). Regional and tempo-ral changes in the pattern of X-chromosome replication during the early post-implantation development of the female mouse.
Chromosoma 85 (2), 275–86.Thompson, S. K. and G. A. Seber (1996).
Adaptive Sampling . John Wiley & Sons.Vines, S. K., W. R. Gilks, and P. Wild (1996). Fitting Bayesian multiple randomeffects models.
Statistics And Computing 6 , 337–346.Wang, X., P. D. Soloway, and A. G. Clark (2010, Jan). Paternally biased X inacti-vation in mouse neonatal brain.
Genome biology 11 (7), R79.Zho, H., Y. Chen, J. Ibrahim, Y. Li, C. Hall, and W. Lin (2010). Intrinsic Regres-sion Models for Positive-Definite Matrices With Applications to Diffusion Tensorimaging.
Journal of American Statistical Association 104 (487), 1203–1212.Zhu, H., J. G. Ibrahim, and N. Tang (2011, May). Bayesian influence analysis: ageometric approach.