Backtrace on a COLO829 sample

[brainstorm]

Hello @jhrf, thanks for this piece of software!

I have been running it for a few samples now, but there's one (non-empty, 22KiB BAM file with only telomeric regions from COLO829) which made it crash this way:

[Error] telomerecat stopped unexpecedtly, sorry!
Traceback (most recent call last):
  File "parabam/core.pyx", line 49, in parabam.core.CmdLineInterface.handle (parabam/core.c:3184)
  File "/home/brainstorm/dev/telomerecat/telomerecat/bam2length.py", line 42, in run_cmd
    simulator_n=self.cmd_args.simulator_runs)
  File "/home/brainstorm/dev/telomerecat/telomerecat/bam2length.py", line 80, in run
    correct_f2a=correct_f2a)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 954, in run
    trim)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 1021, in __get_read_types__
    vital_stats)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 551, in get_read_counts
    self.__paths_to_error_profile__(read_stat_paths)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 570, in __paths_to_error_profile__
    error_profile = self.get_error_profile(read_counts, random_counts)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 586, in get_error_profile
    thresh)
  File "/home/brainstorm/dev/telomerecat/telomerecat/telbam2length.py", line 614, in __get_significantly_enriched__
    thresh = np.percentile(hi_thresh, 95)
  File "/home/brainstorm/.miniconda/lib/python2.7/site-packages/numpy-1.14.0-py2.7-linux-x86_64.egg/numpy/lib/function_base.py", line 4291, in percentile
    interpolation=interpolation)
  File "/home/brainstorm/.miniconda/lib/python2.7/site-packages/numpy-1.14.0-py2.7-linux-x86_64.egg/numpy/lib/function_base.py", line 4033, in _ureduce
    r = func(a, **kwargs)
  File "/home/brainstorm/.miniconda/lib/python2.7/site-packages/numpy-1.14.0-py2.7-linux-x86_64.egg/numpy/lib/function_base.py", line 4405, in _percentile
    x1 = take(ap, indices_below, axis=axis) * weights_below
  File "/home/brainstorm/.miniconda/lib/python2.7/site-packages/numpy-1.14.0-py2.7-linux-x86_64.egg/numpy/core/fromnumeric.py", line 159, in take
    return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)
  File "/home/brainstorm/.miniconda/lib/python2.7/site-packages/numpy-1.14.0-py2.7-linux-x86_64.egg/numpy/core/fromnumeric.py", line 52, in _wrapfunc
    return getattr(obj, method)(*args, **kwds)
IndexError: cannot do a non-empty take from an empty axes.

[jhrf]

Hey, thanks for taking the time to log an issue for this error.

Telomerecat assumes input from a WGS experiment. When creating the error profile it uses some of the non-telomere reads to inform the model. Because you have used a telomere only BAM file it is failing in this assumption. Telomerecat should definitely handle this error better than it currently does. I will look to handle this error better in a future version. However it is also worth considering that the input you have provided is different from the one expected.

From a scientific perspective it's very interesting that you've got targeted telomere sequencing data. Telomerecat could probably be altered to accept these data, if you wish to discuss this then please email me at [email protected].

[brainstorm]

Thanks @jhrf for your feedback. Indeed that was our suspicion from looking at the paper (thanks to @MrsLaviniaG for pointing that out to me). We re-ran with the whole dataset and it seems to output reasonable numbers, although we found some possible floating point issues (hopefully just minor) in one of the runs, see 2.9960000000000004 below:

$ cat telomerecat_length_1519994622.csv
Sample,F1,F2,F4,Psi,Insert_mean,Insert_sd,Read_length,Initial_read_length,F2a,F2a_c,Length
sample,49507,7545,3066,2.9960000000000004,346.0,84.266,112,112,4479,4479,2928.2