Loopless programming (for calculating methylation types)

DNA Methylation is used in plants as a means of epi-genetic regulation. Bi-sulfite sequencing is a method used to determine the methylation pattern of a given set of DNA. Methylation occurs at 'C'ytosines. In plants, there are 3 types of methylation, determined by the nucleotides that follow the 'C':
  • CG: a 'G' follows the 'C'

  • CHG: anything but a 'G' follows the 'C' and a 'G' follows that

  • CHH: no 'G' in the 2 subsequent nucleotides.

So, programmatically, it's a trivial matter to calculate the type of methylation that can occur at each cytosine in a given sequence. Though, once the edges and edge-cases and reverse-complement are handled, it becomes a few lines of code full of loops and if's. For python (and most scripting languages) that's slow and loopy and iffy. The rest of this code explains how I utilized numpy to make a fast methylation-type calculator without loops or ifs.

Loopless


Given a numpy array `seq`, it's possible to find all the 'C's with:
np.where(seq == 'C')
. From there one can calculate the methylation type without looping or if'ing with:

which potentially takes more memory but is extremely fast, and (IMO) nicely shows how to take advantage of the things numpy is good at using the slicing and indexing.
From there, I put that (slightly modified to differentiate between + and - methylation) into a function named _calc_methylation, and call it from this function:

which does some work to handle the ends of the sequence and reverse-complementing. The result, as shown in the docstring, is a numpy array where the values indicate the type of methylation as:

0: not a C or G
1: CG at a 'C' (+ strand )
2: CHG at a 'C' (+ strand )
3: CHH at a 'C' (+ strand )
4: CG at a 'G' (- strand )
5: CHG at a 'G' (- strand )
6: CHH at a 'G' (- strand )

The full implemenation is here. (even the reverse-complement is done without loops or ifs)

Comments

Unknown said…
I like your approach. In a similar spirit this little function finds the first index of every k-tuple in a numpy array. http://python.pastebin.com/f46f7dae7
7:46 PM
brentp said…
@diffusing thoughts . cool. could you show an example usage? i'm not quite sure i follow.
11:04 AM
Unknown said…
It allows to find occurrences of words in a sequence e.g. all start codons in a sequence of letters
>>> sequence = np.array('ATGCGCGTAGCTATGAGAGCATCGAT', dtype='c')
>>> find(sequence, 'ATG')
(array([ 0, 12]),)

The first ATG starts at index 0 the second at index 12. The output from "find" can be used to index an array.

>>> sequence[find(sequence, 'ATG')]
array(['A', 'A'],
dtype='|S1')
12:20 PM
brentp said…
ah, i see. is there a benefit over the
pure python version
?
i guess yours could more easily handle multiple ktup's with numpy broadcasting.
1:13 PM
Unknown said…
The main benefit is that my function generalizes to ndarray i.e. the input could represent an alignment i.e. be a rank-2 numpy array. The other benefit is that it works on numpy arrays as the "find" method is broken in char arrays.

In [1]: "AGCG".find('GC')
Out[1]: 1
In [2]: arr = np.array('AGCG', dtype='c').view(np.char.chararray)
In [3]: arr.find('GC')
Out[3]: array([-1, -1, -1, -1])

char array does not do the right thing because it takes each element from arr e.g. 'A' and calls it's method e.g. 'find' withe the given argument e.g. 'GC', which fails "G".find("GC") returns -1
5:02 AM
Anonymous said…
This comment has been removed by a blog administrator.
11:41 PM
Post a Comment

Popular posts from this blog

python interval tree

EDIT: added a couple points inline. I'm obsessed with trees lately -- of the CS variety, not the plant variety. Although we are studying poplar, so I'll be using trees to study trees . I'd tried a couple times to implement an interval tree from scratch following the Wikipedia entry . Today I more or less did that in python. It's the simplest possible form. There's no insertion (though that's trivial to add), it just takes a list of 'things' with start and stop attributes and creates a tree with a .find() method. The wikipedia entry that baffled me was about storing 2 copies of each node's intervals --one sorted by start and the other by stop. I didn't do that as I think in most cases it won't improve lookup time. I figure if you have 1 million elements and a tree of depth 16, then you have on average 15 intervals per node (actually fewer since there are the non-leaf nodes). So I just brute force each of those nodes and move to the next. I th...
Read more

levenshtein in cython

EDIT(2): fix markup for <char *> casts ... fix malloc (see comments. thanks Bao). NOTE: using a kwarg for limit slows things down. setting that to a required arg and using calloc for m2 speed things up to nearly as fast as the pylevenshtein. Well, it seems to be popular to code up the levenshtein . I actually have a use for this and wanted to practice some Cython , so I've written a version. I used bearophile's recipe , wikibooks , this (from k4st on reddit) and this for reference. It follows bearophile's code closely, using only O(m) space instead of O(mn). cdef extern from "stdlib.h": ctypedef unsigned int size_t size_t strlen(char *s) void *malloc(size_t size) void free(void *ptr) int strcmp(char *a, char *b) cdef inline size_t imin(int a, int b, int c): if a if c return c return a if c return c return b cpdef int levenshtein(char *a, char *b, int limit=100): cdef int m = strlen(a),...
Read more

synteny mapping

Image
Living in Synteny I've been working on automating synteny mapping between any pairs of genomes. Synteny is where there's a stretch of DNA or genes in some order on chromosomeA of organismX and due to a shared evolutionary history, you can find a similar stretch of genes in order on chromosomeB of organismY. Often there are small losses and inversions, but between closely related organisms like man and mouse , there's still a lot of synteny. Plants can undergo polyploidy , following which, a species can have 2 entire copies of its genome. Over time, much of the duplicated cruft is lost, and the homologous chromosomes diverge, but if the divergence is not too great, it's still possible (actually common) to find synteny within the genome of a single organism--as well as between organisms. I've written my own algorithm to find synteny which uses python sets, and numpy array slicing to do the heavy lifting. It is quite clever [wink]. And it _almost_ works but it'...
Read more