嗨，我写了一个脚本，随机地随机排列它们定位到的基因上的阅读序列。如果您要确定在目的基因上观察到的峰是否具有统计学意义，这将很有用。我使用此代码来计算感兴趣基因中峰的错误发现率。代码下方：
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）：
for i in reads: # for read with read length i
r = np.random.randint(-i,featurelength-1,iterations) # generate random read start positions for the read i
for j in b: # for each row in a:
pos = r[j] # get the first random start position for that row
if pos < 0: # start position can be negative because a read does not have to completely overlap with the feature
a[j][:pos+i]+=1
else:
a[j][pos:pos+i]+=1 # add the read to the array and repeat
然后生成一个热图，以查看分布是否大致均匀：
plt.imshow(a)
plt.show()
这会产生所需的结果，但是由于有许多for循环，因此它非常慢。我试图做一些花哨的numpy索引，但是不断出现“索引过多错误”。
有人对如何执行此操作有更好的主意吗？
花式索引有点棘手，但仍然可能：
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
为了对此进行一些解构，我们是：
创建一个简单的索引数组作为从0到i的列向量： np.arange(i)[:,None]
将r（行向量）中的每个元素相加，然后广播以(i,iterations)在的列中制作一个具有正确偏移量的大小矩阵a。
[0,featurelength)通过将索引固定在范围内np.clip。
最后，我们a为每行（b）和相关列（idx）设置了花式索引。
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句