Gh1325 eva igv

[kebarr]

Closes issues 1325 (the main motivation for starting this branch) and these discovered whilst working on fixes: 1357, and 1356.

Re-implement the to_vcf function so that snp_calls from native datasets can be used for visualisation with IGV. Uses Biopython’s bgzf option, the issue with this is that it does not enable the output VCF to be indexed without using external tools. This decision was made because neither pysam nor tabix can be used in Windows without WSL so this was the only way I could find of exposing this functionality to Windows users. This does mean that only small regions can be viewed through the IGV

This implementation only outputs heterozygous calls, this can be made an option if preferable. Missingness, is no longer tested for but can be re-added. The input has changed to require a single sample ID, as currently the only use case for this method is for IGV. These are breaking changes compared to the previous to_vcf method, but this was not used anywhere in the code base, so will not affect any other methods. It also now requires the output filename to end with “vcf.gz,”

Tests added in test_vcf_exporter and demo added to igv.ipynv - rather than replace the current IGV methods, the local VCF format has been added as an alternative to grabbing one from the external URL. For the IGV notebook to run, the VCFs must be written to a path within the “notebooks” directory, this is explained in the notebook.

Validation of output VCFs was carried out by comparing to VCF linked in wgs_data_catalogue. The results are not identical because the variant calling pipeline used a site list frozen at the time it was called, whereas this just pulls SNPs from the Zarr dataset for the specific sample. They were confirmed using tools from ”bedtools” and “bcftools” intersect for individual chromosomes and the full ag3 genome. There are n sites in active disagreement.

Prior to PR, full test suite was run, and coverage, reported as 98% for to_vcf and 96% for igv.

[kebarr]

Just waiting for other tests to finished, but one failure:

"FAILED tests/anoph/test_snp_frq.py::test_allele_frequencies_advanced_with_sample_query[af1_sim] - ValueError: variant_position_aa_change must not be empty"

from tests(3.11,==2.0.2) appears to be related to code I haven't touched. Still waiting on the others. @jonbrenas @ahernank is this expected?

Will update this comment if any further fail in similar ways.

[jonbrenas]

Thanks @kebarr. The other tests were cancelled on failure. I restarted the test that failed, I assume it is a transient error that resulted from the randomly selected region containing no non-synonymous mutations. We may want to update the tests at some point to treat such a lack of data as a warning instead of an error, but it probably requires a bit of thinking.

[kebarr]

Please find attached screenshots obtained from the new igv.ipynb - showing agreement between the SNP calls pulled from the Zarr array and the VCF. They look different because homozygous calls are omitted from my VCF, including them would make the resource constraints of using an unindexed VCF even more severe, to the point of risking making the entire endeavor worthless.

I haven't yet had the opportunity to do proper profiling to determine exactly how many base pairs can be supported for viewing in IGV- it will be machine dependent, so showing vairation in RAM/CPU wrt base pairs in selected region is probably the best approach (though not ideal as SNPs will not be uniformly distributed), enabling users to determine what can be supported on their machine. It appears that RAM will be the largest determining factor for visualisation with IGV. I assume (based on absolutely no evidence) clock speed, number CPUs etc, impacts pulling out the SNP calls more so I'll separate snp_calls_to_vcf and IGV viewing for performance testing. Once we have something we're happy with for this type of testing it should be more generally useful.

I have also attached a very basic validation script to show how I have checked my sno_calls_to_vcf output to the reference VCF. If you would prefer me to convert this into a more intelligible report or make it more flexible (i.e. user can specify region or dataset) please let me know.

[validate_vcf.sh](https://github.com/user-attachments/files/33164889/validate_vcf.sh) [make_vcf_from_snp_calls.py](https://github.com/user-attachments/files/33164894/make_vcf_from_snp_calls.py)

[jonbrenas]

Thanks @kebarr. Should we try to find some time next week for you to walk through what you have done?

In the meantime, it looks like a lot of notebooks that don't have much to do with what you were doing were added to this PR. Could you do some clean up?

[kebarr]

@jonbrenas I've booked a meeting for Monday.

The changes to the notebooks are a side-effect of me adding nbstripout to my pre-commit-hooks, because it was difficult looking through massive diffs whilst working on the notebook side. I'm in travel training all day but will figure out how to remove those from this PR tomorrow, and look into having a "personal" pre-commit hook with it or another way of making it easier.

[kebarr]

I did some profiling of the two steps required for my approach vs the previous approach. The previous approach includes indexed VCFs so lower memory consumption is expected. It appears, and makes sense that, the updated version uses more memoory when more SNPs are present in the region. Which would make guiding users much more difficult than just saying "don't load regions more than X bps."

I've also attached the notebooks used for profiling + script used for making papermill run them in batches. Working on a more generalisable script.

test_capture_mem_output.ipynb
10_samples_10_increments_igv_views.csv
10_samples_10_increments_snp_calls.csv
papermill_test_script.py

10_samples_10_increments_default_vcf.csv
test_capture_metrics_default.ipynb