In [1]:

    

import numpy as np
np.random.seed(42)
import sys
import cProfile
import h5py
sys.path.insert(0, '../..')
%reload_ext memory_profiler
%reload_ext autoreload
%autoreload 1
import allel; print(allel.__version__)
%aimport allel.stats.selection


    

    




0.21.0.dev0

In [2]:

    

%matplotlib inline
import matplotlib.pyplot as plt


    

In [3]:

    

callset = h5py.File('/data/coluzzi/ag1000g/data/phase1/release/AR3/haplotypes/main/hdf5/ag1000g.phase1.ar3.haplotypes.3R.h5',
                    mode='r')
callset


    

    
Out[3]:





<HDF5 file "ag1000g.phase1.ar3.haplotypes.3R.h5" (mode r)>

In [4]:

    

n_variants = 500000
n_samples = 100


    

In [5]:

    

h = allel.GenotypeArray(callset['3R']['calldata/genotype'][:n_variants, :n_samples]).to_haplotypes()
h


    

    
Out[5]:





HaplotypeArray((500000, 200), dtype=int8)



0
1
2
3
4
...
195
196
197
198
199




0
0
0
0
0
0
...
0
0
0
0
0


1
0
0
0
0
0
...
0
0
0
0
0


2
0
0
0
0
0
...
0
0
0
0
0


3
0
0
0
0
0
...
0
0
0
0
0


4
0
0
0
0
0
...
0
0
0
0
0



...

In [6]:

    

pos = callset['3R']['variants/POS'][:n_variants]
pos


    

    
Out[6]:





array([   1252,    1262,    1271, ..., 1927216, 1927221, 1927229], dtype=int32)

In [7]:

    

ac = h.count_alleles(max_allele=1)
is_seg = ac.is_segregating() & (ac.min(axis=1) > 1)
h_seg = h.compress(is_seg, axis=0)
pos_seg = pos.compress(is_seg)
ac_seg = ac.compress(is_seg, axis=0)
np.count_nonzero(is_seg)


    

    
Out[7]:





106315

In [8]:

    

%%time
nsl = allel.stats.nsl(h_seg, use_threads=False)


    

    




CPU times: user 6.58 s, sys: 0 ns, total: 6.58 s
Wall time: 6.57 s

In [10]:

    

%%time
nsl_threaded = allel.stats.nsl(h_seg, use_threads=True)


    

    




CPU times: user 6.74 s, sys: 0 ns, total: 6.74 s
Wall time: 3.41 s

In [11]:

    

nsl


    

    
Out[11]:





array([ 2.06801504,  3.38963917,  0.3487763 , ...,  0.08136895,
        0.93199028,  1.31336717])

In [12]:

    

nsl_threaded


    

    
Out[12]:





array([ 2.06801504,  3.38963917,  0.3487763 , ...,  0.08136895,
        0.93199028,  1.31336717])

In [13]:

    

np.count_nonzero(np.isnan(nsl)), np.count_nonzero(~np.isnan(nsl))


    

    
Out[13]:





(0, 106315)

In [14]:

    

np.count_nonzero(np.isinf(nsl)), np.count_nonzero(~np.isinf(nsl))


    

    
Out[14]:





(0, 106315)

In [15]:

    

plt.figure(figsize=(16, 4))
plt.plot(pos_seg, nsl, linestyle=' ', marker='o', mfc='none')
plt.grid(axis='y')
plt.xlabel('position (bp)')
plt.ylabel('unstandardised NSL');


    

In [16]:

    

plt.figure(figsize=(16, 6))
plt.plot(ac_seg[:, 1], nsl, linestyle=' ', marker='o', mfc='none')
plt.xlabel('alternate allele count')
plt.ylabel('unstandardised NSL score')
plt.grid(axis='y');


    

In [17]:

    

nsl_standardized, ac_bins = allel.stats.standardize_by_allele_count(nsl, ac_seg[:, 1])


    

In [18]:

    

plt.figure(figsize=(16, 6))
plt.plot(ac_seg[:, 1], nsl_standardized, linestyle=' ', marker='o', mfc='none')
plt.xlabel('Alternate allele count')
plt.ylabel('Standardised NSL score')
plt.grid(axis='y');


    

In [19]:

    

plt.figure(figsize=(16, 4))
plt.plot(pos_seg, nsl_standardized, linestyle=' ', marker='o', mfc='none')
plt.grid(axis='y')
plt.xlabel('Position (bp)')
plt.ylabel('Standardised NSL score');


    

In [20]:

    

plt.figure(figsize=(16, 4))
plt.plot(pos_seg, np.abs(nsl_standardized), linestyle=' ', marker='o', mfc='none')
plt.grid(axis='y')
plt.xlabel('Position (bp)')
plt.ylabel('|Standardised NSL score|');


    

In [21]:

    

cProfile.run('allel.stats.nsl(h_seg[:50000], use_threads=False)', sort='time')


    

    




         92 function calls in 3.404 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        2    3.401    1.701    3.401    1.701 stats.pyx:803(nsl01_scan_int8)
        1    0.002    0.002    3.403    3.403 selection.py:562(nsl)
        1    0.000    0.000    3.404    3.404 {built-in method builtins.exec}
        6    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.array}
        1    0.000    0.000    3.403    3.403 <string>:1(<module>)
        2    0.000    0.000    0.000    0.000 ndarray.py:1810(__getitem__)
        4    0.000    0.000    0.000    0.000 stringsource:985(memoryview_fromslice)
       12    0.000    0.000    0.000    0.000 stringsource:341(__cinit__)
        5    0.000    0.000    0.000    0.000 numeric.py:414(asarray)
        2    0.000    0.000    3.401    1.701 {allel.opt.stats.nsl01_scan_int8}
        8    0.000    0.000    0.000    0.000 stringsource:643(memoryview_cwrapper)
        1    0.000    0.000    0.000    0.000 ndarray.py:1776(__new__)
        3    0.000    0.000    0.000    0.000 ndarray.py:1783(__array_finalize__)
       12    0.000    0.000    0.000    0.000 stringsource:368(__dealloc__)
        3    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        4    0.000    0.000    0.000    0.000 stringsource:962(__dealloc__)
        2    0.000    0.000    0.000    0.000 ndarray.py:1765(_check_input_data)
        1    0.000    0.000    0.000    0.000 {method 'view' of 'numpy.ndarray' objects}
        4    0.000    0.000    0.000    0.000 stringsource:545(__get__)
        1    0.000    0.000    0.000    0.000 <frozen importlib._bootstrap>:996(_handle_fromlist)
        8    0.000    0.000    0.000    0.000 stringsource:649(memoryview_check)
        4    0.000    0.000    0.000    0.000 stringsource:507(__getbuffer__)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        3    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
        1    0.000    0.000    0.000    0.000 {method 'setdefault' of 'dict' objects}

In [22]:

    

ihs_min0 = allel.stats.ihs(h_seg, pos_seg, min_ehh=0.0, include_edges=True)


    

In [23]:

    

plt.figure(figsize=(8, 8))
plt.plot(ihs_min0, nsl, marker='o', linestyle=' ', mfc='none')
plt.xlabel('Unstandardized IHS (min_ehh=0)')
plt.ylabel('Unstandardized NSL')
plt.grid(axis='both');


    

In [24]:

    

ihs = allel.stats.ihs(h_seg, pos_seg, min_ehh=0.05, include_edges=True)


    

In [25]:

    

plt.figure(figsize=(8, 8))
plt.plot(ihs, nsl, marker='o', linestyle=' ', mfc='none')
plt.xlabel('Unstandardized IHS (min_ehh=0.05)')
plt.ylabel('Unstandardized NSL')
plt.grid(axis='both');


    

In [26]:

    

ihs_standardized, _ = allel.stats.standardize_by_allele_count(ihs, ac_seg[:, 1])


    

In [27]:

    

plt.figure(figsize=(8, 8))
plt.plot(ihs_standardized, nsl_standardized, marker='o', linestyle=' ', mfc='none')
plt.xlabel('Standardized IHS (min_ehh=0.05)')
plt.ylabel('Standardized NSL')
plt.grid(axis='both');


    

In [28]:

    

loc_variants = slice(4000000, 9000000, 1)
n_samples = 50

h = allel.GenotypeArray(callset['3R']['calldata/genotype'][loc_variants, :n_samples]).to_haplotypes()
h


    

    
Out[28]:





HaplotypeArray((5000000, 100), dtype=int8)



0
1
2
3
4
...
95
96
97
98
99




0
0
0
0
0
0
...
0
0
0
0
0


1
0
0
0
0
0
...
0
0
0
0
0


2
0
0
0
0
0
...
0
0
0
0
0


3
0
0
0
0
0
...
0
0
0
0
0


4
0
0
0
0
0
...
0
0
0
0
0



...

In [29]:

    

pos = callset['3R']['variants/POS'][loc_variants]
pos


    

    
Out[29]:





array([16532852, 16532853, 16532854, ..., 42555205, 42555210, 42555213], dtype=int32)

In [30]:

    

ac = h.count_alleles(max_allele=1)
is_seg = ac.is_segregating() & (ac.min(axis=1) > 10)
h_seg = h.compress(is_seg, axis=0)
pos_seg = pos.compress(is_seg)
ac_seg = ac.compress(is_seg, axis=0)
np.count_nonzero(is_seg)


    

    
Out[30]:





211040

In [43]:

    

def plot_score_gap(score, pos, ylim=(-10, 10)):

    fig = plt.figure(figsize=(16, 4))

    ax = fig.add_subplot(111)
    ax.plot(pos, score, linestyle=' ', marker='o', mfc='none', markersize=2)
    ax.grid(axis='y')
    ax.set_xlabel('position (bp)')
    ax.set_ylabel('score')
    ax.set_ylim(*ylim)

    ax = ax.twinx()
    x = (pos[:-1] + pos[1:]) / 2
    y = np.diff(pos)
    ax.plot(x, y)
    ax.set_ylabel('gap size (bp)')
    ax.autoscale(axis='x', tight=True);


    

In [44]:

    

nsl = allel.stats.nsl(h_seg)
nsl_standardized, _ = allel.stats.standardize_by_allele_count(nsl, ac_seg[:, 1])


    

In [45]:

    

nsl.shape, nsl_standardized.shape


    

    
Out[45]:





((211040,), (211040,))

In [46]:

    

plot_score_gap(nsl, pos_seg)


    

In [47]:

    

plot_score_gap(nsl_standardized, pos_seg)


    

In [48]:

    

plot_score_gap(np.abs(nsl_standardized), pos_seg, ylim=(0, 10))


    

In [49]:

    

accessibility = h5py.File('/data/coluzzi/ag1000g/data/phase1/release/AR3/accessibility/accessibility.h5', mode='r')
is_accessible = accessibility['3R']['is_accessible'][:]
is_accessible


    

    
Out[49]:





array([False, False, False, ..., False, False, False], dtype=bool)

In [52]:

    

ihs_access_adjusted = allel.stats.ihs(h_seg, pos_seg, min_ehh=0.05, max_gap=-1, clip_gap=-1, 
                                      is_accessible=is_accessible)
ihs_access_adjusted, _ = allel.stats.standardize_by_allele_count(ihs_access_adjusted, ac_seg[:, 1])


    

In [57]:

    

plot_score_gap(np.abs(nsl_standardized), pos_seg, ylim=(0, 10))
plt.title('|NSL|');


    

In [58]:

    

plot_score_gap(np.abs(ihs_access_adjusted), pos_seg, ylim=(0, 12))
plt.title('|IHS|');


    

In [ ]: