利用numpy随机分配基因组特征上的DNA序列读数

Question

嗨，我写了一个脚本，随机地把读到的序列按照它们所映射到的基因进行排列。 如果你想确定你在感兴趣基因上观察到的一个峰值在统计学上是否显著，这是很有用的。我用这个代码来计算我感兴趣基因峰值的错误发现率。 以下代码：

import numpy as np
import matplotlib.pyplot as plt
iterations = 1000 # number of times a read needs to be shuffled
featurelength = 1000  # length of the gene
a = np.zeros((iterations,featurelength))  # create a matrix with 1000 rows of the feature length
b = np.arange(iterations)                 # a matrix with the number of iterations (0-999)
reads = np.random.randint(10,50,1000)     # a random dataset containing an array of DNA read lengths

下面的代码填充大矩阵（a）：

^{pr2}$

然后生成一个热图，看看分布是否大致均匀：

plt.imshow(a)
plt.show()

这将生成所需的结果，但由于for循环太多，因此速度非常慢。 我试图做花哨的纽比索引，但我经常得到“太多索引错误”。在

有人知道怎么做吗？在

Answer 1

花式索引有点棘手，但仍有可能：

for i in reads:
    r = np.random.randint(-i,featurelength-1,iterations)
    idx = np.clip(np.arange(i)[:,None]+r, 0, featurelength-1)
    a[b,idx] += 1

要稍微解释一下，我们是：

1. 创建一个简单的索引数组作为列向量，从0到i:np.arange(i)[:,None]

2. 将r（一个行向量）中的每个元素相加，该元素广播以使一个大小为(i,iterations)的矩阵具有正确的偏移量到a的列中。

3. 通过np.clip，将索引限制在[0,featurelength)范围内。

4. 最后，我们为每一行（b）和相关列（idx）设置索引a。