Genomic arrays

A recurring task in genomics is mapping some quantity onto a coordinate system: coverage along a chromosome, the density of a feature, a score accumulated across overlapping annotations. The operation is always the same — take an interval, add a value to every position inside it, repeat a few million times.

The class below does that and nothing else. It allocates a zeroed NumPy array the length of the chromosome and exposes one method for adding a value across a range. NumPy handles the slice arithmetic, so the accumulation stays fast even when the array is a quarter of a billion positions long. You need NumPy installed to run it.

One warning worth internalising: this follows Python’s indexing rules. Positions are zero-based and the last coordinate is not included — the interval is open on the right. If your annotations came from a one-based, closed-interval format such as GFF, convert them before you get here rather than after you have plotted the result.

import numpy as np


class GenomicArray():
    def __init__(self, size):
        self.size = size
        self.store = np.zeros(self.size)

    def chunk_add(self, s, e, val):
        self.s = s
        self.e = e
        self.val = val
        self.store[self.s:self.e] += self.val
        return self.store


# length of array, e.g. for human chromosome 1
chr1_length = 247249719

# initialise chr1_array
chr1_array = GenomicArray(chr1_length)

chr1_array.chunk_add(3, 5, 10)
chr1_array.chunk_add(4, 10, 20)

# see that the values changed
print(chr1_array.store[0:20])

Running that gives 10 at position 3, 30 at position 4 where the two intervals overlap, and 20 from position 5 through 9. Position 10 stays at zero — the right-open rule doing its job, and the thing that will catch you out if you have spent the morning in a one-based format.

Why not a dictionary

The obvious alternative is a dictionary keyed on position, and for sparse data it is the better choice. But once features overlap across a significant fraction of a chromosome, the dictionary loses on both memory and speed, and you give up the ability to slice, window and summarise with NumPy directly.

A flat array of float64 across human chromosome 1 costs about two gigabytes. If that is more than you want to spend, cast the store to a smaller dtype when you construct it.