2

我有一个名为“数据”的结构化数组，同一个条目有几个分数。为了这个问题，我将“数据”减少到以下两列。

queryid    bitscore
gene1      500
gene1      480
gene1      440
gene2      900
gene2      300

我想要做的是提取相同的queryid 的最高值，即bitscore 至少比最高bitscore 低10% 的任何常见条目。

例如，只有前 2 个条目“gene1”应该被保留，因为第三个条目的比特分数低于 500 的 10%。对于基因 2，只有第一个应该被保留（这个很容易）。

queryid    bitscore
gene1      500
gene1      480
gene2      900

当我做一个这样的循环时：

for i in range(0, lastrow-1, 1):
if data[i]['queryid'] == data[i+1]['queryid']:
    if data[i+1]['bitscore'] < data[i]['bitscore']-(0.01*data[i]['bitscore']):
       data[i+1]['queryid'] = 'DELETE'

data = data[data[:]['queryid'] != 'DELETE']

所有“gene1”条目将被保留，因为 440 在 480 的 10% 范围内。

我可以将最高值添加到可以保留作为参考的另一列，但我想检查你们中是否有人对此有更好的想法......

4

3 回答 3

3

如果你能够使用 pandas，它就变成了一个单行问题：

from pandas import DataFrame
import numpy as np

# Taken from Theodros
data = zip(('gene1',) * 3 + ('gene2',) * 2, [500, 480, 440, 900, 300])
dtype = [('queryid', 'S6'), ('bitscore', 'i4')]
struct_arr = np.array(data, dtype=dtype)

# Create pandas DataFrame from NumPy struct array
df = DataFrame.from_records(struct_arr)

# Filter the rows per group
df.groupby('queryid').apply(lambda x: x[x["bitscore"] >= x["bitscore"].max() * 0.9])

产生：

          queryid  bitscore
queryid                    
gene1   0   gene1       500
        1   gene1       480
gene2   3   gene2       900
于 2013-07-17T22:21:44.137 回答
2

使用逻辑索引可能比for循环快得多。像这样的东西怎么样：

def high_bitscores(a,qid,thresh=0.9):
    valid = a[a['queryid'] == qid]
    return valid[valid['bitscore'] >= valid['bitscore'].max()*thresh]

编辑：如果要返回通过此标准的所有元素，data则可以遍历其中的唯一queryid值data并更新一组布尔索引，以指定哪些元素通过了测试：

def all_high_bitscores(a,thresh=0.9):

    # set of indices for the elements in a that we're going to keep
    keep = np.zeros(a.size,np.bool)

    for qid in set(a['queryid']):
        idx = a['queryid'] == qid
        keep[idx] = a[idx]['bitscore'] >= a[idx]['bitscore'].max()*thresh

    return a[keep]
于 2013-07-16T23:56:13.827 回答
1

您可以使用itertools.groupby来迭代具有相同queryid. 然后你可以过滤掉那些低于阈值的记录。在此步骤中，您将保留 numpy.structured_array 数据结构并处理单个元组。通过定义filter为生成器，您可以从过滤后的输出动态构建新的结构化数组。

from itertools import groupby
from operator import itemgetter


def filter(struct_arr, threshold):
    for k, g in groupby(struct_arr, key=itemgetter(0)):
        ref = g.next()
        yield ref
        ref = ref[1]
        for e in g:
            if (float(e[1]) / ref) < threshold:
                break
            yield e

例子：

data = zip(('gene_1',) * 3 + ('gene_2',) * 2, [500, 480, 440, 900, 300])
dtype = [('query_id', 'S6'), ('bitscore', 'i4')]
struct_arr = np.array(data, dtype=dtype)


np.fromiter(filter(struct_arr), dtype=dtype)

这使

array([('gene_1', 500), ('gene_1', 480), ('gene_2', 900)], 
      dtype=[('query_id', 'S6'), ('bitscore', '<i4')])
于 2013-07-16T23:56:32.947 回答