在python中有效地取稀疏数据的移动平均值并过滤阈值以上

Question

我的脚被一些基因组分析弄湿了，有点卡住了。我有一些非常稀疏的数据，需要找到移动平均线超过某个阈值的地方，将每个点标记为 1 或 0。数据是唯一类型的，因此我无法使用可用的程序进行分析。

每个点代表人类基因组上的一个点（碱基对）。对于每个数据集，有 200,000,000 个潜在点。数据本质上是一个约 12000 个索引/值对的列表，其中所有其他点都假定为零。我需要做的是在整个数据集上取一个移动平均值，并返回平均值高于阈值的区域。

我目前正在按顺序从数据集中读取每个点，并围绕我找到的每个点构建一个数组，但这对于大窗口大小来说非常慢。有没有更有效的方法来做到这一点，也许是用 scipy 或熊猫？

编辑：下面杰米的魔术代码效果很好（但我还不能投票）！我非常感激。

Answer 1

您可以使用 numpy 对整个事物进行矢量化。我构建了这个随机数据集（大约）12,000 个介于 0 和 199,999,999 之间的索引，以及一个同样长的介于 0 和 1 之间的随机浮点列表：

indices = np.unique(np.random.randint(2e8,size=(12000,)))
values = np.random.rand(len(indices))

2*win+1然后，我围绕每个构建一个总窗口大小的索引数组indices，以及一个相应的数组，表示该点对移动平均值的贡献量：

win = 10

avg_idx = np.arange(-win, win+1) + indices[:, None]
avg_val = np.tile(values[:, None]/(2*win+1), (1, 2*win+1))

剩下的就是计算重复指数并将对移动平均线的贡献加在一起：

unique_idx, _ = np.unique(avg_idx, return_inverse=True)
mov_avg = np.bincount(_, weights=avg_val.ravel())

您现在可以获得移动平均线超过 0.5 的指数列表，如下所示：

unique_idx[mov_avg > 0.5]

至于性能，首先将上面的代码变成一个函数：

def sparse_mov_avg(idx, val, win):
    avg_idx = np.arange(-win, win+1) + idx[:, None]
    avg_val = np.tile(val[:, None]/(2*win+1), (1, 2*win+1))
    unique_idx, _ = np.unique(avg_idx, return_inverse=True)
    mov_avg = np.bincount(_, weights=avg_val.ravel())
    return unique_idx, mov_avg

以下是一些窗口大小的时序，用于开头描述的测试数据：

In [2]: %timeit sparse_mov_avg(indices, values, 10)
10 loops, best of 3: 33.7 ms per loop

In [3]: %timeit sparse_mov_avg(indices, values, 100)
1 loops, best of 3: 378 ms per loop

In [4]: %timeit sparse_mov_avg(indices, values, 1000)
1 loops, best of 3: 4.33 s per loop