Skip to content

Error trying to get Fst for window_by_position #1343

Description

@alxsimon

Hi,
I keep hitting this issue that I don't manage to debug.
The window_by_variant and window_by_genome work fine.
I suspect this could be due to some windows not having any variant? But in that case I'd simply expect a nan value for such windows.

ds_fst = (
    ds
    .pipe(lambda ds: sg.window_by_position(ds, size=50_000))
    .pipe(sg.divergence)
    .pipe(sg.Fst)
)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[5], line 4
      1 ds_fst = (
      2 	ds
      3 	.pipe(lambda ds: sg.window_by_position(ds, size=50_000))
----> 4 	.pipe(sg.divergence)
      5 	.pipe(sg.Fst)
      6 )

File ~/miniforge3/envs/popgenpy/lib/python3.12/site-packages/xarray/core/common.py:872, in DataWithCoords.pipe(self, func, *args, **kwargs)
    869     kwargs[target] = self
    870     return f(*args, **kwargs)
--> 872 return func(self, *args, **kwargs)

File ~/coding/sgkit/sgkit/stats/popgen.py:233, in divergence(ds, cohort_allele_count, nansum, merge)
    230 sumfn = np.nansum if nansum else np.sum
    232 if has_windows(ds):
--> 233     div = window_statistic(
    234         d,
    235         sumfn,
    236         ds.window_start.values,
    237         ds.window_stop.values,
    238         dtype=d.dtype,
    239         axis=0,
    240     )
    241     new_ds = create_dataset(
    242         {
    243             variables.stat_divergence: (
   (...)    247         }
    248     )
    249 else:

File ~/coding/sgkit/sgkit/window.py:570, in window_statistic(values, statistic, window_starts, window_stops, dtype, chunks, new_axis, **kwargs)
    568     # new chunks are same except in first axis
    569     new_chunks = tuple([tuple(windows_per_chunk)] + list(desired_chunks[1:]))  # type: ignore
--> 570 return da.map_overlap(
    571     blockwise_moving_stat,
    572     values,
    573     dtype=dtype,
    574     chunks=new_chunks,
    575     depth=depth,
    576     boundary=0,
    577     trim=False,
    578     new_axis=new_axis,
    579 )

File ~/miniforge3/envs/popgenpy/lib/python3.12/site-packages/dask/array/overlap.py:722, in map_overlap(func, depth, boundary, trim, align_arrays, allow_rechunk, *args, **kwargs)
    717 args = [
    718     overlap(x, depth=d, boundary=b, allow_rechunk=allow_rechunk)
    719     for x, d, b in zip(args, depth, boundary)
    720 ]
    721 assert_int_chunksize(args)
--> 722 x = map_blocks(func, *args, **kwargs)
    723 assert_int_chunksize([x])
    724 if trim:
    725     # Find index of array argument with maximum rank and break ties by choosing first provided

File ~/miniforge3/envs/popgenpy/lib/python3.12/site-packages/dask/array/core.py:866, in map_blocks(func, name, token, dtype, chunks, drop_axis, new_axis, enforce_ndim, meta, *args, **kwargs)
    850     out = blockwise(
    851         apply_and_enforce,
    852         out_ind,
   (...)    863         **kwargs,
    864     )
    865 else:
--> 866     out = blockwise(
    867         func,
    868         out_ind,
    869         *concat(argpairs),
    870         name=name,
    871         new_axes=new_axes,
    872         dtype=dtype,
    873         concatenate=True,
    874         align_arrays=False,
    875         adjust_chunks=adjust_chunks,
    876         meta=meta,
    877         **kwargs,
    878     )
    880 extra_argpairs = []
    881 extra_names = []

File ~/miniforge3/envs/popgenpy/lib/python3.12/site-packages/dask/array/blockwise.py:272, in blockwise(func, out_ind, name, token, dtype, adjust_chunks, new_axes, align_arrays, concatenate, meta, *args, **kwargs)
    270 elif isinstance(adjust_chunks[ind], (tuple, list)):
    271     if len(adjust_chunks[ind]) != len(chunks[i]):
--> 272         raise ValueError(
    273             f"Dimension {i} has {len(chunks[i])} blocks, adjust_chunks "
    274             f"specified with {len(adjust_chunks[ind])} blocks"
    275         )
    276     chunks[i] = tuple(adjust_chunks[ind])
    277 else:

ValueError: Dimension 0 has 1114 blocks, adjust_chunks specified with 4456 blocks

sgkit version: latest from https://github.com/alxsimon/sgkit
dask version: '2024.8.0'

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions