windowed_r_squared: IndexError for single-SNP windows

[pmonnahan]

Hi,

I have been trying to implement the windowed_r_squared function, and I believe that I've come across a bug related to windows containing a single SNP. Below is a snippet generated via pdb, showing that the error arises within the 'windowed_statistic' function. I noticed in the source code of 'windowed_statistic' that there is an 'if n==0' catch and am wondering if this needs to be generalized for r2 calculations where a pair of observations is required at a minimum.

I could increase window size or try to identify and remove these windows beforehand, but I thought I would check whether this is a known (or non-) issue on your end.

r2, windows, counts = allel.windowed_r_squared(pos, ct, size=winsize, start=start, stop=stop, step=int(winsize / 2), fill=-9)` (Pdb) s --Call-- /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/ld.py(161)windowed _r_squared() -> def windowed_r_squared(pos, gn, size=None, start=None, stop=None, step=None, (Pdb) /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/ld.py(214)windowed _r_squared() -> if isinstance(percentile, (list, tuple)): (Pdb) /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/ld.py(222)windowed _r_squared() -> def statistic(gnw): (Pdb) /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/ld.py(226)windowe$ _r_squared()

... After about 5 windows with error-free calculation, I encounter this:

return windowed_statistic(pos, gn, statistic, size, start=start, (Pdb) s wv = values[start_idx:stop_idx] (Pdb) n /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/window.py(366)windowed_statistic() s = statistic(wv) (Pdb) print(n) 1 (Pdb) print(wv) [[0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 2 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0]] (Pdb) wv.shape (1, 99) (Pdb) n IndexError: cannot do a non-empty take from an empty axes. /home/spectorl/pmonnaha/.conda/envs/py36/lib/python3.6/site-packages/allel/stats/window.py(366)windowed_statistic() -> s = statistic(wv)

Thanks!

[hardingnj]

Thanks @pmonnahan .

This looks to reproduce with:

import allel
import numpy as np

pos = np.random.choice(100_000, 100, replace=False)

gn = np.random.randint(0, 3, (100, 100))

r2, windows, counts = allel.windowed_r_squared(pos, gn, size=5000, start=1, stop=100_000, fill=-9)

And this line (equivalent) causes the error:

np.percentile(np.array([]), 50)

I guess a check for the shape of gnw is the best option here.

[pmonnahan]

Thanks @hardingnj for the reprex and marking for change. Also, thanks to you and others for such a great package.