luntergroup / octopus

Bayesian haplotype-based mutation calling
MIT License
299 stars 37 forks source link

BadRegion error #214

Closed brentp closed 2 years ago

brentp commented 2 years ago

Describe the bug

Encountered an exception during calling 'BadRegion'. This means there is a bug and your results are untrustworthy.

Version

$ octopus --version
  Using htslib 1.12
  Copyright (C) 2021 Genome Research Ltd.
  License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
  This is free software: you are free to change and redistribute it.
  There is NO WARRANTY, to the extent permitted by law.
  octopus version 0.7.4
  Target: x86_64 Linux 5.4.0-72-generic
  SIMD extension: AVX2
  Compiler: GNU 9.3.0
  Boost: 1_74

Command Command line to install octopus:

conda

Command line to run octopus:

octopus -R $ref \                                                               
    -p Y=2 chrY=2 chrM=1 chrMT=1 MT=1 --threads ${task.cpus} --one-based-indexing \
    -T ${region} \                                                              
    --debug /hpc/compgen/users/bpedersen/${reg}.debug.log \                     
    --disable-denovo-variant-discovery \                                        
    -i $workDir/${reg}.${index}.crams.list \                                    
    --source-candidates-file $workDir/${reg}.${index}.vcfs.list \               
    -o ${output_path}  
octopus -R GRCh38_full_analysis_set_plus_decoy_hla.fa     -p Y=2 chrY=2 chrM=1 chrMT=1 MT=1 --threads 10 --one-based-indexing     -T chr2:1-100000000     --disable-denovo-variant-discovery     -i /hpc/compgen/users/bpedersen/genetics-first-rdsolve-data/work/chr2_1-100000000.1.crams.list     --source-candidates-file /hpc/compgen/users/bpedersen/genetics-first-rdsolve-data/work/chr2_1-100000000.1.vcfs.list     -o chr2_1-100000000.1.population.vcf.gz

the end of the output looks like this:

[2021-10-12 10:39:08] <INFO>   chr2:15009603            15.0%           2m 29s            14m 6s
  [2021-10-12 10:39:16] <INFO>   chr2:16009658            16.0%           2m 37s           13m 46s
  [2021-10-12 10:39:20] <INFO>               -             100%           2m 42s                 -
  [2021-10-12 10:39:20] <INFO> Removed 3 temporary files
  [2021-10-12 10:39:21] <EROR> A program error has occurred:
  [2021-10-12 10:39:21] <EROR> 
  [2021-10-12 10:39:21] <EROR>     Encountered an exception during calling 'BadRegion'. This means
  [2021-10-12 10:39:21] <EROR>     there is a bug and your results are untrustworthy.
  [2021-10-12 10:39:21] <EROR> 
  [2021-10-12 10:39:21] <EROR> To help resolve this error run in debug mode and send the log file to
  [2021-10-12 10:39:21] <EROR> https://github.com/luntergroup/octopus/issues.
  [2021-10-12 10:39:21] <INFO> ------------------------------------------------------------------------

Additional context

I am splitting by region and some regions are succeeding, but this one fails each time.

The debug log is too large to add here, but I can share with you via box or other if you send your email address.

dancooke commented 2 years ago

Thanks for the bug report. Could you try re-running with a single thread (makes the log easier to follow) and send the resulting debug log to dcooke@well.ox.ac.uk?

BTW, it shouldn't be directly causing this bug, but htslib 1.12 has a bug that can cause missing variants.

brentp commented 2 years ago

Thanks for the reply. This could take some time, but i will set it running. Ouch! didn't know about that bug. Looks like it's fixed in 1.13+?

dancooke commented 2 years ago

Yeah - fixed in 1.13.

brentp commented 2 years ago

I won't hijack tihs issue, except to note that the reason I had 1.12 was this:

$ conda env create -n oct13 --file octopus-env.yaml 
Collecting package metadata (repodata.json): done
Solving environment: / 
Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abort.
failed                                                                                                                                                          
Solving environment: - 
Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abort.
failed                                                                                                                                                          

UnsatisfiableError: The following specifications were found to be incompatible with each other:

Output in format: Requested package -> Available versions

Package htslib conflicts for:
htslib[version='>=1.13']
octopus=0.7.4 -> htslib[version='>=1.12,<1.13.0a0']
bcftools[version='>=1.13'] -> htslib[version='>=1.13,<1.14.0a0']The following specifications were found to be incompatible with your system:

  - feature:/linux-64::__glibc==2.17=0
  - feature:|@/linux-64::__glibc==2.17=0
  - bcftools[version='>=1.13'] -> libgcc-ng[version='>=9.3.0'] -> __glibc[version='>=2.17']
  - octopus=0.7.4 -> libgcc-ng[version='>=9.3.0'] -> __glibc[version='>=2.17']

Your installed version is: 2.17

Note that strict channel priority may have removed packages required for satisfiability.
$ cat octopus-env.yaml 
name: octopus
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  - bcftools>=1.13
  - htslib>=1.13
  - octopus=0.7.4
dancooke commented 2 years ago

Thanks - I'll take a look at the conda issue.

The bug is coming from the variant filtering step, could you rerun your command (can be with multithreading) but add the --disable-call-filter option. If you're able to send over the resulting unfiltered VCF, then I should be able to work it out from there. If you're not able to, then can you subset at chr2:16,273,306-16,331,770 - which looks to be the problematic region?

brentp commented 2 years ago

Hi, I can't share the VCF, but it completed with the new flags you suggested, then I tried to tabix it and get:

[E::hts_idx_push] Unsorted positions on sequence #1: 16332262 followed by 16332260
tbx_index_build failed: chr2_1-100000000.1.population.vcf.gz

indeed, if those positions are out of order while all others are ordered (I checked by changing and re-trying the indexing).

Those 2 variants and the one before look like this:

chr2    16332262        .       T       G       1382.15 .       AC=2;AN=4;DP=125;MQ=59;NS=3     GT:GQ:DP:MQ:PS:PQ       1|0:501:36:59:16332155:13       1|.:293:40:59:16332155:12       0|.:755:49:59:16332155:46
chr2    16332260        .       G       T       1049.19 .       AC=2;AN=4;DP=124;MQ=59;NS=3     GT:GQ:DP:MQ:PS:PQ       .|0:501:35:59:16332155:13       .|1:293:40:59:16332155:12       0|1:755:49:59:16332155:46

so (I'm sure if I can diagnose, then you'll be able to, but I can't help but guess), perhaps you are decomposing and normalizing and not re-sorting after... and the filtering is relying on them to be sorted.

dancooke commented 2 years ago

That'll be the problem - this has previously come up (https://github.com/luntergroup/octopus/issues/176). The final fix was merged into develop at 6f4cbc051962a31d636819e0f92e5f9f33e16103 - unfortunately that's after v0.7.4 so isn't on Bioconda. Happy to make a new release but building Bioconda releases has been a bit of a pain (https://github.com/bioconda/bioconda-recipes/pull/27169). I push develop to DockerHub regularly if that's an option for you.

brentp commented 2 years ago

cheers. would you consider a sandybridge-dev tag to dockerhub? if not, I can build, but that would be convenient.