API: xsmc package

class xsmc.XSMC(ts: dataclasses.InitVar, focal: int, panel: List[int], theta: float = None, rho_over_theta: float = 1.0, w: int = None, robust: bool = False)

Sample from, or maximize, the posterior distribution \(p(X_{1:N} | \mathbf{Y})\) of the sequentially Markov coalescent.

Parameters

  • ts – Tree sequence containing the data.

  • focal – Leaf node in ts corresponding to focal haplotype.

  • panel – Leaf node(s) in ts corresponding to panel haplotypes.

  • theta – Population-scaled mutation rate. If None, Watterson’s estimator is used.

  • rho_over_theta – Ratio of recombination to mutation rates.

  • w – Window size. Observations are binned into windows of this size. Recombinations are assumed to occur between adjacent bins, but not within them. If None, try to calculate a sensible default based on rho.

  • robust – If True, use robust model in which bins are classified as either segregating or nonsegregating, as in PSMC. Otherwise, allow for any number of mutations per bin.

Notes

The height of segments returned by this class are expressed in coalescent units. To convert to generations, rescale them by 2 * self.theta / (4 * mu), where mu is the biological mutation rate.

sample(k: int = 1, seed: int = None) → List[xsmc.segmentation.Segmentation]

Sample path(s) from the posterior distribution.

Parameters

  • k – Number of posterior path samples to draw, default 1.

  • seed – Random seed used for sampling.

Returns

A list of k posterior samples for the positional min-TMRCA of focal with panel.

Notes

If sampling many paths at once, it is more efficient to set k > 1 than to call sample() repeatedly.

viterbi(beta: float = None) → xsmc.segmentation.Segmentation

Compute the maximum a posteriori (a.k.a. Viterbi) path in haplotype copying model.

Parameters

beta – Penalty parameter for changepoint formation. See manuscript for details about this parameter.

Returns

A segmentation containing the MAP path.

class xsmc.segmentation.Segmentation(segments: List[xsmc.segmentation.Segment], panel: List[int])

A path through a haplotype copying model.

Parameters

  • segments – list of segments composing the path.

  • panel_size – size of panel used to construct this path.

draw(axis: matplotlib.axis.Axis = None, **kwargs) → None

Plot this segmentation as a step function.

Parameters

  • axis – Axis on which to draw plot. If None, current axis is used.

  • kwargs – Additional arguments passed to axis.plot().

Notes

Only plots segment heights. Segment haplotype is currently ignored.

classmethod from_ts(ts: tskit.TreeSequence, focal: int, panel: List[int]) → Segmentation

Return a segmentation consisting of the genealogical MRCA of focal among panel.

Parameters

  • ts – Tree sequence containing the data.

  • focal – Focal leaf node in ts.

  • panel – Panel leaf nodes in ts over which to compute GMRCA.

Returns

Computed segmentation.

Notes

The returned segmentation is not necessarily unique. If there is >1 GMRCA for a given IBD segment, one is chosen arbitrarily.

property panel

Alias for field number 1

rescale(x: float) → xsmc.segmentation.Segmentation

Return a new segmentation where each height is multiplied by x.

property segments

Alias for field number 0

to_pp() → scipy.interpolate.interpolate.PPoly

Return a piecewise polynomial representation of this segmentation.

Notes

Only represents segments heights. Identity of haplotype panel is currently ignored.