将numpy.searchsorted方法应用到使用numpy.loadtxt从文本文件加载的数组中

Question

我目前正在研究生物信息学项目，我需要解决以下问题。

我有一个包含两列的文本文件“chr1.txt”：染色体上的位置和布尔变量True或False。

0虚假
10000真正
10001真正
10005虚假
10007真正
10011虚假
10013真正
10017虚假
10019虚假
10023虚假
10025真正
10029真
10031 false
10035真
10037假
....
该数据意味着，从0到10000的区域重复或（=不可映射 - >假），从10000到10005是唯一的（=可映射 - >真），从10005到10007再次重复，等等。该文件在248'946'406位置结束，并具有15'948'271行。要找到问题的一般解决方案，我想限制文件到上面可以看到的行。

我想加载这个文本文件到由两列组成的numpy数组。对于我用numpy.loadtxt：

import numpy as np  
with open('chr1.txt','r') as f: 
     chr1 = np.loadtxt(f, dtype={'names':('start','mappable'), 
     'formats':('i4','S1')}) 

这里是输出：

In [39]: chr1 
Out[39]: 
array([(0, b'f'), (10000, b't'), (10001, b't'), (10005, b'f'), 
     (10007, b't'), (10011, b'f'), (10013, b't'), (10017, b'f'), 
     (10019, b'f'), (10023, b'f'), (10025, b't'), (10029, b't'), 
     (10031, b'f'), (10035, b't'), (10037, b'f')], 
     dtype=[('position start', '<i4'), ('mappable', 'S1')]) 

，因为我想第二列被识别为布尔类型这看起来并不完美，但我没有找到办法。

从明年我要扔在一个随机数位10000和10037.之间

In [49]: np.random.randint(10000,10037) 
Out[49]: 10012 

现在我想的numpy.searchsorted方法适用于我的阵列的第一列以找出是否我的基因组在该位置是唯一可映射的。所以我想在这种情况下输出为5（我的数组中的元素（10011，b'f'）的索引）。如果我试图提取数组仅包括第一列 - 位置，我得到一个错误：

In [21]: chr1[:,0] 
--------------------------------------------------------------------------- 
IndexError        Traceback (most recent call last) 
<ipython-input-21-a63d052f1c5d> in <module>() 
----> 1 chr1[:,0] 

IndexError: too many indices for array 

我想这是因为我的数组并没有真正有两列

In [40]: chr1.shape 
Out[40]: (15,) 

那么我怎么才能提取只有位置和使用我现有的数组应用搜索排序方法？我是否应该以不同的方式将我的文本文件加载到数组中，以便确实有两列，第一列是整数类型，第二列 - 布尔值？

​​

然后我会看一个创建索引的第二个参数是真还是假的，并能够做出结论，如果位置是可映射的区域内。

非常感谢您的帮助！

Answer 1

我们可以提取对应于position start的数据，其中chr1['position start']和第二个字段的数据类似。我们将得到有效的布尔数组，并与't'进行比较。

因此，我们将有一个方法中，像这样 -

indx = chr1['position start'] 
mask = chr1['mappable']=='t' 
rand_num = np.random.randint(10000,10037) 
matched_indx = np.searchsorted(indx, rand_num)-1 

if mask[matched_indx]: 
    print "It is mappable!" 
else: 
    print "It is NOT mappable!" 

1）获取的数据和掩模/布尔数组 -

In [283]: chr1 # Input array 
Out[283]: 
array([( 0, 'f'), (10000, 't'), (10001, 't'), (10005, 'f'), 
     (10007, 't'), (10011, 'f'), (10013, 't'), (10017, 'f'), 
     (10019, 'f'), (10023, 'f'), (10025, 't'), (10029, 't'), 
     (10031, 'f'), (10035, 't'), (10037, 'f')], 
     dtype=[('position start', '<i4'), ('mappable', 'S1')]) 

In [284]: indx = chr1['position start'] 
    ...: mask = chr1['mappable']=='t' 
    ...: 

In [285]: indx 
Out[285]: 
array([ 0, 10000, 10001, 10005, 10007, 10011, 10013, 10017, 10019, 
     10023, 10025, 10029, 10031, 10035, 10037], dtype=int32) 

In [286]: mask 
Out[286]: 
array([False, True, True, False, True, False, True, False, False, 
     False, True, True, False, True, False], dtype=bool) 

2）获取的随机数，并使用searchsorted和使用IF-ELSE部分 -

In [297]: rand_num = 10012 # np.random.randint(10000,10037) 

In [298]: matched_indx = np.searchsorted(indx, rand_num)-1 

In [299]: matched_indx 
Out[299]: 5 

In [300]: if mask[matched_indx]: 
    ...:  print "It is mappable!" 
    ...: else: 
    ...:  print "It is NOT mappable!" 
    ...:  
It is NOT mappable!