search

alxsimon / local_pcangsd

Local PCA with low-coverage data
https://alxsimon.github.io/local_pcangsd/_build/html/index.html
GNU General Public License v3.0
4 stars 1 forks source link

LinAlgError with pca_window() #2

Closed TeresaPegan closed 1 year ago

TeresaPegan commented 1 year ago

Hello! Thanks so much for developing this script, it has been great for me to be able to used lostruct on my low coverage data!

I am working with a bunch of different species and for some of them, I have edited their beagle files to remove certain individuals. I run into an error with local_pcangsd with these edited beagles, and I was wondering if you have any ideas! The error is triggered by pca_window() and I've been having trouble piecing apart the code to troubleshoot it because it involves the window-based zarr methods, which I am pretty unfamiliar with.

I don't get this error with my un-edited beagles, and I have not yet found any differences between the beagles that work and those that do not, other than the fact that I removed some columns (particular individuals) from the ones that don't work.

Let me know if you have any thoughts on what I could try! Thanks so much again! -Teresa

first couple lines of a beagle file that works:

marker  allele1 allele2 Ind0    Ind0    Ind0    Ind1    Ind1    Ind1    Ind2    Ind2    Ind2    Ind3    Ind3    Ind3    Ind4    Ind4    Ind4    Ind5    Ind5    Ind5    Ind6    Ind6    Ind6    Ind7    Ind7    Ind7    Ind8    Ind8    Ind8    Ind9    Ind9    Ind9    Ind10   Ind10   Ind10   Ind11   Ind11   Ind11   Ind12   Ind12   Ind12   Ind13   Ind13   Ind13   Ind14   Ind14   Ind14   Ind15   Ind15   Ind15   Ind16   Ind16   Ind16   Ind17   Ind17   Ind17   Ind18   Ind18   Ind18   Ind19   Ind19   Ind19   Ind20   Ind20   Ind20   Ind21   Ind21   Ind21   Ind22   Ind22   Ind22   Ind23   Ind23   Ind23   Ind24   Ind24   Ind24   Ind25   Ind25   Ind25   Ind26   Ind26   Ind26   Ind27   Ind27   Ind27   Ind28   Ind28   Ind28   Ind29   Ind29   Ind29   Ind30   Ind30   Ind30   Ind31   Ind31   Ind31   Ind32   Ind32   Ind32   Ind33   Ind33   Ind33   Ind34   Ind34   Ind34   Ind35   Ind35   Ind35   Ind36   Ind36   Ind36   Ind37   Ind37   Ind37   Ind38   Ind38   Ind38   Ind39   Ind39   Ind39   Ind40   Ind40   Ind40   Ind41   Ind41   Ind41   Ind42   Ind42   Ind42   Ind43   Ind43   Ind43   Ind44   Ind44   Ind44   Ind45   Ind45   Ind45   Ind46   Ind46   Ind46   Ind47   Ind47   Ind47   Ind48   Ind48   Ind48   Ind49   Ind49   Ind49   Ind50   Ind50   Ind50   Ind51   Ind51   Ind51   Ind52   Ind52   Ind52   Ind53   Ind53   Ind53
CM020336.1_336  3   1   0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.000044    0.333333    0.666622    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333
CM020336.1_361  1   3   0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.000044    0.333333    0.666622    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333

first couple lines of a beagle that does not work:

marker  allele1 allele2 Ind0    Ind0    Ind0    Ind1    Ind1    Ind1    Ind2    Ind2    Ind2    Ind3    Ind3    Ind3    Ind4    Ind4    Ind4    Ind5    Ind5    Ind5    Ind6    Ind6    Ind6    Ind7    Ind7    Ind7    Ind9    Ind9    Ind9    Ind10   Ind10   Ind10   Ind11   Ind11   Ind11   Ind14   Ind14   Ind14   Ind16   Ind16   Ind16   Ind18   Ind18   Ind18   Ind19   Ind19   Ind19   Ind21   Ind21   Ind21   Ind22   Ind22   Ind22   Ind23   Ind23   Ind23   Ind24   Ind24   Ind24   Ind25   Ind25   Ind25   Ind26   Ind26   Ind26   Ind28   Ind28   Ind28   Ind29   Ind29   Ind29   Ind30   Ind30   Ind30   Ind31   Ind31   Ind31   Ind32   Ind32   Ind32   Ind33   Ind33   Ind33   Ind35   Ind35   Ind35   Ind36   Ind36   Ind36   Ind38   Ind38   Ind38   Ind39   Ind39   Ind39   Ind40   Ind40   Ind40   Ind41   Ind41   Ind41   Ind42   Ind42   Ind42   Ind43   Ind43   Ind43   Ind45   Ind45   Ind45   Ind47   Ind47   Ind47   Ind48   Ind48   Ind48   Ind49   Ind49   Ind49   Ind50   Ind50   Ind50   Ind51   Ind51   Ind51   Ind52   Ind52   Ind52
CM020336.1_336  3   1   0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.000044    0.333333    0.666622    0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333
CM020336.1_361  1   3   0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.000044    0.333333    0.666622    0.333333    0.333333    0.333333    0.666649    0.333333    0.000018    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333    0.333333

local_pcangsd code up to the point where it hits the error:

import os
import local_pcangsd as lp
import lostruct
from skbio.stats.ordination import pcoa
import matplotlib.pyplot as plt
import seaborn as sns
import dask
import pandas as pd
import numpy as np
import xarray as xr
import sys
from matplotlib.backends.backend_pdf import PdfPages

input = sys.argv[1] 
store = sys.argv[2]

print(input)
print(store)

lp.beagle_to_zarr(input, store, chunksize=100000)

print('finished beagle to zarr!')

ds = lp.load_dataset(store, chunks=100000) # open the Dataset
ds

ds = lp.window(ds, type='position', size=50000, min_variant_number = 500)
# the position argument makes it so that windows are given in their positions
ds

# record the window loci
windf = pd.DataFrame({'window_start' : ds.window_start.to_numpy(),
              'window_stop' : ds.window_stop.to_numpy(),
             })

print('starting pca')

#%%time
pca_zarr_store = lp.pca_window(
    ds,
    zarr_store=sys.argv[3], # where to store the result
    tmp_folder=sys.argv[4], # need a tmp folder, /tmp/tmp_local_pcangsd is default
    k=10, # number of PCs to retain
)

The error:

Traceback (most recent call last):
  File "/home/tmpegan/WGLC_Boreal/Nov2022/lostruct_pcangsd.py", line 44, in <module>
    pca_zarr_store = lp.pca_window(
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/local_pcangsd/local_pcangsd.py", line 411, in pca_window
    zarr_list = result.compute(
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/base.py", line 288, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/base.py", line 570, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/threaded.py", line 79, in get
    results = get_async(
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/local.py", line 507, in get_async
    raise_exception(exc, tb)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/local.py", line 315, in reraise
    raise exc
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/local.py", line 220, in execute_task
    result = _execute_task(task, data)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/core.py", line 119, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/optimization.py", line 969, in __call__
    return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/core.py", line 149, in get
    result = _execute_task(task, cache)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/core.py", line 119, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/dask/array/core.py", line 456, in _pass_extra_kwargs
    return func(*args[len(keys) :], **kwargs)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/local_pcangsd/local_pcangsd.py", line 378, in blockwise_moving_stat
    res_ij = _pcangsd_wrapper(
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/local_pcangsd/local_pcangsd.py", line 195, in _pcangsd_wrapper
    vals, vectors = np.linalg.eig(C)
  File "<__array_function__ internals>", line 5, in eig
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/numpy/linalg/linalg.py", line 1317, in eig
    _assert_finite(a)
  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/numpy/linalg/linalg.py", line 208, in _assert_finite
    raise LinAlgError("Array must not contain infs or NaNs")
numpy.linalg.LinAlgError: Array must not contain infs or NaNs
alxsimon commented 1 year ago

Hi Teresa, First thank you for trying out this module, I didn't really expect someone would start to use it at this point. I'm still trying to improve it and it is a slow work in progress, I'm glad I can have some feedback!

I've encountered this kind of error before but in another context when trying to create simulated beagle files without adding uncertainty in the genotype likelihoods (1 0 0 for example).

I think the cause of the error is that in this line, the C covariance matrix sent back by pcangsd have NaNs. But I am not sure what the root cause is in pcangsd. I'll try to think how I can improve error handling so that we can get more information.

Would you be able to send me a minimal non-working example containing 2 or 3 windows?

In the meantime, it might also be interesting to run the beagle file with pcangsd directly to try to see if it gives similar issues.

Alexis

alxsimon commented 1 year ago

@Rosemeis would you have any lead on where I should start looking for issues of this type? I feel like I am missing a filtering step in the data I pass to pcangsd.covariance.emPCA but not sure where to start. Thanks!

alxsimon commented 1 year ago

I have a lead, I think I simply forgot to filter for MAF, so if any fixed markers happen this should lead to errors...

TeresaPegan commented 1 year ago

Ah, that would perfectly explain why it happens when I delete a subset of individuals! Markers that were polymorphic might be fixed in the reduced subset.

alxsimon commented 1 year ago

let me try to fix this and get back to you

alxsimon commented 1 year ago

Could you try out the maf branch with your data? I've added a min_maf argument to the pca_window function. I haven't found a good way to store which loci are used or not in the final result though, I'll need to come back to that later.

Rosemeis commented 1 year ago

Yes the MAF would be my main suspect, or potentially that all individuals have missing data (all 0.333). :-)

Best, Jonas

TeresaPegan commented 1 year ago

It works on the maf branch! Thanks so much!

alxsimon commented 1 year ago

Great, I'll try to find a way to store the loci list info and merge this in main.

TeresaPegan commented 1 year ago

Another really small bug that's probably not worth its own issue -- on line 570 of local_pcangsd.py, the window start indices have to be ds.window_start.values, not just ds.window_start, otherwise I get this error:

  File "/home/tmpegan/anaconda3/envs/local_pcangsd/lib/python3.9/site-packages/local_pcangsd/local_pcangsd.py", line 570, in get_window_center
    pos_start = ds.variant_position.values[ds.window_start]
IndexError: too many indices for array: array is 1-dimensional, but 16 were indexed

The definition of pos_stop is formatted correctly as ds.variant_position.values[ds.window_stop.values - 1] , so just a typo or something. Thanks! -Teresa

alxsimon commented 1 year ago

I've now merged the changes into main. There are some breaking changes though:

Don't hesitate to open a new issue if there are other problems.

TeresaPegan commented 1 year ago

Hi @alxsimon, I am going to defend my dissertation soon and I have a chapter that uses this repository (hopefully to be submitted as a paper not long after). I wanted to ask if you have a preference about how this repository is cited and whether you have any publications/preprints on it. I can at least link the repository. Thanks! -Teresa

alxsimon commented 1 year ago

Hi Teresa, for the time being you can just link the repository. I have some plan to publish it in some form or another (preprint, zenodo doi or final publication) but haven't put the time in it yet. Don't hesitate to fill me in on the advancement of your manuscript, this will be an incentive for me to finish things up ;)

TeresaPegan commented 6 months ago

Hi @alxsimon, I just wanted to let you know that I have a manuscript in review that references this repository. Keep me up to date if you have a preprint or forthcoming ms to watch out for! Thanks, -Teresa

alxsimon commented 6 months ago

Hi Teresa, unfortunately I won't have more time to publish anything on it. Referencing the zenodo archive is enough on my side. Congrats on the manuscript and happy this could be useful to someone. Best, Alexis