Local PCA with low-coverage data
Home Page: https://alxsimon.github.io/local_pcangsd/_build/html/index.html
License: GNU General Public License v3.0
Find me on Mastodon: https://ecoevo.social/@alxsim
Admin of ecoevo.social Mastodon server.
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
In pca_window rename zarr_store to store.
beagle is not the most useful output of ANGSD, bcf is nicer for downstream manipulation.
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.