TrimOverlap: RuntimeError: Unexpected sample relationship

[wdecoster]

Hi,

I encountered the following error in the TrimOverlap step, using medaka v0.7.0-alpha.1:

Traceback (most recent call last):
  File "/home/wdecoster/anaconda3/envs/medaka/bin/medaka", line 10, in <module>
    sys.exit(main())
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/medaka.py", line 350, in main
    args.func(args)
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/variant.py", line 472, in variants_from_hdf
    vcf_writer.write_variants(variants, sort=True)
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/vcf.py", line 272, in write_variants
    variants = medaka.common.loose_version_sort(variants, key=lambda v: '{}-{}'.format(v.chrom, v.pos))
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/common.py", line 708, in loose_version_sort
    it = list(it)
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/variant.py", line 206, in decode_variants
    for s, _ in yield_trimmed_consensus_chunks(samples):
  File "/home/wdecoster/anaconda3/envs/medaka/lib/python3.6/site-packages/medaka/variant.py", line 69, in yield_trimmed_consensus_chunks
    raise RuntimeError('Unexpected sample relationship {} between {} and {}'.format(repr(rel), s1.name, s2.name))
RuntimeError: Unexpected sample relationship <Relationship.reverse_overlap: 'The end of s2 overlaps the start of s1.'> between chr1:198506700.10723-198510472.2 and chr1:198506700.1723-198506883.0
Failed to call variants from consensus chunks.

Please let me know if there is a file which I can share to help you debug this issue. The bam file I'm using is 80x coverage of the human genome...

Cheers,
Wouter

[cjw85]

Wouter,
If you can create a chr1-only bam and its not too big, I would like to have a look.

One curiosity, though not necessarily wrong and something medaka should handle, is that you have (at least one) read with a very long insert: the decimal places in the region string in your log above indicates insertion columns with respect to the reference.

[cjw85]

I'm just reviewing the changes we made between the pre-release and v0.7.0.

There was a commit which altered how chunks are sorted when they are read back from the HDF datastore. It looks to me in your case the chunks have indeed been read back incorrectly giving the Relationship.reverse_overlap error.

[wdecoster]

Could a structural variant explain the very long insert?

While attempting to share the data there are a few hurdles. One is "NA19240-chr1.cram has an invalid extension. Valid extension(s): zip, gz, 7z, tar, tgz, sqfs, hdf, fast5, xz." But okay, I can rename that one. But second: "I confirm that this data is not human genomic data, and that any necessary agreements and consents required to share this data with Oxford Nanopore are in place." -=> it's human, but I can share as it's from NA19240.

[cjw85]

I have run your data using the v0.7.0 release, It took around 30 mins to run the initial medaka consensus pass using a GPU and a further 10 minutes to create an unphased vcf. I believe the error you saw was therefore corrected between the pre-release and the release.

[wdecoster]

Alright then, I'll update again and see how it goes. Thanks!