wrong sequence length after keep_interval()

17 views
Skip to first unread message

Miguel Navascués

unread,
Oct 1, 2020, 10:03:27 AM10/1/20
to msprim...@googlegroups.com
Hi,

I am using SLiM and pyslim+msprime (for recapitation) for simulating
whole human genomes. I want to apply msprime independently to smaller
chuncks of genome (e.g. chromosome arms) and wanted to use
keep_interval() as suggested by Peter Ralph
(https://groups.google.com/forum/#!topic/slim-discuss/UIwmOPuv9Ww).

However, I am not sure I am doing it right, the resulting tree sequence
seems to be of the same size of the original sequence:

>>> import msprime
>>> tree_sequence = msprime.simulate(sample_size=6, length=50e3, Ne=1000)
>>> tree_sequence.sequence_length
50000.0
>>> tree_sequence.keep_intervals(np.array([[0,10]])).sequence_length
50000.0

Shouldn't this give me a sequence length of 10?

I am using
Python 3.6.9
tskit 0.2.3


--
Miguel NAVASCUÉS

UMR CBGP, INRAE
Centre de Biologie pour la Gestion des Populations
755 avenue du campus Agropolis
CS30016
34988 Montferrier-sur-Lez cedex (France)

phone: +33499623370
fax: +33499623345
e-mail: miguel.navascues AT inrae.fr

Peter Ralph

unread,
Oct 1, 2020, 10:52:56 AM10/1/20
to Miguel Navascués, msprime-users
Hi, Miguel - it's true, what you expect would also make sense, and
this is something we debated when writing this method. However, the
way this method works is to *maintain* the coordinates, but just
remove all the *information* (trees and mutations) outside of the
given intervals.
(see https://github.com/tskit-dev/tskit/pull/261#issuecomment-511207084
for our rationale). To do what you expect, you should follow it up
with an ltrim and/or rtrim:
https://tskit.readthedocs.io/en/latest/python-api.html?#tskit.TreeSequence.ltrim

This should be explained in the documentation for keep_intervals:
apologies for the confusion!

- peter
> --
> You received this message because you are subscribed to the Google Groups "msprime-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to msprime-user...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/msprime-users/f704096a-3a68-2dd5-1886-d08fbca489f5%40inrae.fr.
Reply all
Reply to author
Forward
0 new messages