title: using numpy to randomly distribute DNA sequence reads over genomic features
body: <p>Hi I have written a script that randomly shuffles read sequences over the gene they were mapped to.
  This is useful if you want to determine if a peak that you observe over your gene of interest is statistically significant. I use this code to calculate False Discovery Rates for peaks in my gene of interest.
  Below the code:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>Below the code to fill the large matrix (a):</p>
  
  <pre><code>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 &lt; 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
  </code></pre>
  
  <p>Then generate a heat map to see if the distribution is roughly even:</p>
  
  <pre><code>plt.imshow(a)
  plt.show()
  </code></pre>
  
  <p>This generates the desired result but it is very slow because of the many for loops.
  I tried to do fancy numpy indexing but I constantly get the "too many indices error".</p>
  
  <p>Anybody have a better idea of how to do this?</p>