Skip to content

Deal with multiple contigs and sequence lengths #964

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Feb 25, 2025

Conversation

hyanwong
Copy link
Member

@hyanwong hyanwong commented Sep 4, 2024

Introduces a contig_id parameter to variant_data, which adds non matching contigs to the site mask, as described in #949. Fixes #949

Copy link

codecov bot commented Sep 4, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.38%. Comparing base (1aa0233) to head (d5dd3f4).
Report is 1 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #964      +/-   ##
==========================================
+ Coverage   93.31%   93.38%   +0.07%     
==========================================
  Files          18       18              
  Lines        6460     6475      +15     
  Branches     1096     1101       +5     
==========================================
+ Hits         6028     6047      +19     
+ Misses        294      291       -3     
+ Partials      138      137       -1     
Flag Coverage Δ
C 93.38% <100.00%> (+0.07%) ⬆️
python 95.77% <100.00%> (+0.09%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@hyanwong
Copy link
Member Author

hyanwong commented Sep 7, 2024

Also fixes #763 I think

@hyanwong hyanwong force-pushed the multi-contig branch 4 times, most recently from aa733d5 to 8742ee8 Compare September 10, 2024 17:58
@hyanwong hyanwong requested a review from benjeffery September 10, 2024 18:12
@benjeffery
Copy link
Member

benjeffery commented Feb 18, 2025

Thanks @hyanwong - I've been some thinking on API design here and am thinking a simpler solution is to error out if there is more than one contig in the used variants.
As for contig_length, I don't think this is the same conceptually as tskit sequence_length? As in tskit positions have to be less than sequence length, whereas there is no such requirement in vcf zarr, so we risk causing odd (to the user) errors downstream if we use contig_length. So how about we check if contig_length is greater than the positions and error early if it isn't, and allow the user to override with a sequence_length argument. If both are missing use the largest position+1.

@hyanwong
Copy link
Member Author

Sorry, brief reply while away, but this came from user feedback in tsinfer: people want to be able to take a vcf with multiple chromosomes in it (which is the common input format), select a chromosome, and run inference without too much extra thought (e.g. they are probably not the original generators of the vcf, and may not want to have to parse the vcf themselves). I assume this would include getting the seq len from the vcf if that is specified? But whatever makes it least error-prone for a naive user to run tsinfer on an externally provided vcf file, I guess?

@hyanwong
Copy link
Member Author

I think checking contig length if present and erroring out early if it is problematic seems fine?

@benjeffery
Copy link
Member

people want to be able to take a vcf with multiple chromosomes in it (which is the common input format), select a chromosome, and run inference without too much extra thought

I think this scheme enables that - you bio2zarr then write a contig-specific mask, which is a one-liner we can put in the tutorial. I don't want to add things like contig_id to the API when they overlap and interact with existing functionality (masking).

@andrewkern
Copy link
Member

Thanks @hyanwong - I've been some thinking on API design here and am thinking a simpler solution is to error out if there is more than one contig in the used variants. As for contig_length, I don't think this is the same conceptually as tskit sequence_length? As in tskit positions have to be less than sequence length, whereas there is no such requirement in vcf zarr, so we risk causing odd (to the user) errors downstream if we use contig_length. So how about we check if contig_length is greater than the positions and error early if it isn't, and allow the user to override with a sequence_length argument. If both are missing use the largest position+1.

Having the option to override the sequence_length is a feature that I'm looking for in a project that i'm doing right now. I've been using tsinfer to create tree sequences from windowed chunks of a VCF file that I will analyze downstream in a ML pipeline.

@benjeffery
Copy link
Member

How does this look @hyanwong

@hyanwong
Copy link
Member Author

LGTM. Fine to merge, I reckon.

@hyanwong
Copy link
Member Author

(p.s. #765 suggests that the error message which currently reads

"Sites belong to multiple contigs ({contig_names}). "
"Please restrict sites to one contig using the sites_mask argument."

Should have a code snippet to show users how to do this?

@benjeffery
Copy link
Member

Should have a code snippet to show users how to do this?

Good idea, added

@benjeffery
Copy link
Member

@jeromekelleher This is now ready for your review.

Copy link
Member

@jeromekelleher jeromekelleher left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM, spotted one minor issue

):
self._sequence_length = sequence_length
self._contig_index = None
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to define self._contig_id here also, or will get attribute error when not defined.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed, and I've added checking this to the tests.

@mergify mergify bot merged commit d8dafc7 into tskit-dev:main Feb 25, 2025
14 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

VariantData sequence_length not being taken from zarr contig_length field
4 participants