Pattern matching with degenerate characters and don't-cares using bit-vector encoding

The pattern matching problem falls into two practical categories: exact and approximate pattern matching.

In pattern matching we have a text T of length n and a pattern X of length m ≤ n. In exact pattern matching we want to find X in T exactly as it is. In approximate pattern matching we want to find X in T but it doesn't need to match exactly. We define a threshold value k which is the maximum number of optimal changes that can be made to X to make it match a part of T. If more than k changes are required to make X match a part of T then we conclude that it is not in fact a factor of T and move on to look at the next position.

Exact pattern matching is easily solved using the naive method in worst-case O(nm) time. Of course, better algorithms exist to reduce this complexity. The KMP algorithm is famous but it is not the best one to use. Other algorithms include Shift-And/Shift-Or, Horspool, BNDM and BOM, of which the last two are the most effective. They achieve their efficiency using a combination of text skipping - moving past parts of the text that is impossible to hold the pattern - and bit-parallel manipulation. Bit-parallelism takes advantage of the computer processor's ability to quickly recall items from registers and cache and manipulate it in parallel using simple shift operations and logic gate calculations. Exact pattern matching is also quickly achieved using on-line indexing methods such as Aho-Corasick Automata or the suffix tree/array approach which speed up multiple pattern look up.

Approximate pattern matching traditionally relies on Dynamic Programming approaches which have an O(nm) time complexity. Improvements often rely on bit-parallelism, avoiding unnecessary processing by filtering out parts of the text and various clever tricks relying on low k.

In this blog post I will show how to simulate approximate pattern matching using the naive exact pattern matching algorithm. One presumes it is possible to use better exact pattern-matching algorithms to achieve faster results, but the goal of this post is not to make the ultimate algorithm. Instead, it is to explain how to use bit-parallelism to achieve pattern matching for degenerate and don't-care characters.

I am using the DNA/RNA alphabet for this example. DNA consists of the characters A, C, G, and T. RNA substitutes T for U, making its alphabet A, C, G, U. The characters represent nucleotide bases, for example, A stands for Adenine. We introduce a few extra degenerate characters and a don't care character (N) as defined in the IUPAC Alphabet standard:

A.................Adenine
C.................Cytosine
G.................Guanine
T.................Thymine
U.................Uracil
R.................A or G
Y.................C or T
S.................G or C
W.................A or T
K.................G or T
M.................A or C
B.................C or G or T
D.................A or G or T
H.................A or C or T
V.................A or C or G
N.................any base

I call a character like R or B degenerate, meaning it can represent more than one single alphabet character. I call N a don't care character because it can represent any base. Note that there are 4 core letters (A[T|U]CG), assuming that T is equivalent to U and 15 characters all together considering. So why do we have this special IUPAC standard? Well, one reason is because there are enzymes that recognise a segment of DNA with some flexibility. The restriction enzyme EcoRII's specificity subunit recognises and cuts DNA when it finds the site CCWGG. W could be A or T. It is useful for biologists to know where it makes its cuts in the genome. Another reason the IUPAC alphabet might be useful is to define DNA sequences where the exact identity of a base is unclear or may contain errors when it is first sequenced.

We can represent every IUPAC character using a binary representation different to the ASCII standard. Through this binary representation and clever use of bit masks it is possible to make R = A or R = G, or N = any base.

First I create an array containing 122 elements whose index match to the ASCII numerical value of all the letters above. For example 'A' is index 65, 'B' is index 66 and so on. I could have made this a smaller array taking only the range A - Y, but I do not like using subtract all the time in my code. Thus, most of the positions will hold the value 0, but for A, I define its value to be 1, 0001 in binary. In fact, I define the 4 primary nucleotide bases as such:

NUC['A'] = 1; //00000001
NUC['C'] = 2; //00000010
NUC['G'] = 4; //00000100
NUC['T'] = 8; //00001000
NUC['U'] = 8; //00001000

Notice that the position of the binary ones make each character distinct. For each distinct character a new bit is needed. In our case with an effective alphabet size of 4, it can all fit into the first nibble of a byte.

Now we design the degenerate characters by making them bit-masks of the characters they can represent.

NUC['R'] = 5;  //00000101 - a/g
NUC['Y'] = 10; //00001010 - c/t
NUC['S'] = 6;  //00000110 - g/c
NUC['W'] = 9;  //00001001 - a/t
NUC['K'] = 12; //00001100 - g/t
NUC['M'] = 3;  //00000011 - a/c
NUC['B'] = 14; //00001110 - c/g/t
NUC['D'] = 13; //00001101 - a/g/t
NUC['H'] = 11; //00001011 - a/c/t
NUC['V'] = 7;  //00000111 - a/c/g
NUC['N'] = 15; //00001111 - a/c/g/t

Before we can use this, I need to describe how to find an exact pattern in a text using the regular method. Let's say we have a text T = "GACCAGGAG" and a pattern X = "CCWGG" and we need to find X in T. But because the naive method cannot recognise the degenerate W character in T, it would not be sensible to use X as is. In fact, you would be forced to do two searches, once using CCAGG and the next time CCTGG. This is the naive search method:

i = 0
while i <= n - m:
    j = 0
    while j < m and T[i+j] == X[j]:
        j += 1
        if j == m:
            print "Match found at position: ", i
    i += 1

We use the variable i to keep track of the position in T and j the position in X. The outer while loop moves the pointer i through T one position at a time. So for every position in T, the inner while loop tries to move through X. It only continues to loop through X as long as the current character in T and subsequent characters match that of X until j is equal to the length of the pattern, in which case we record a match. If before j == m there is a mismatch, it breaks out of the inner loop and moves to the next position in T to keep looking for the pattern.

OK. If that is sufficiently clear, we can modify the above exact pattern matching strategy to match degenerate characters and don't cares based on the following principles.

Each character A[T|U]CG in the NUC array we defined above has a unique binary code with the position of its binary 1 differing. We use logical AND (&) to check if one character matches another.

  NUC['A']: 00000001
& NUC['A']: 00000001
--------------------
            00000001

  NUC['A']: 00000001
& NUC['C']: 00000010
--------------------
            00000000

So when A matches A, the resulting character is A and when it doesn't, it gives a different character (0). With degenerate characters we do:

  NUC['A']: 00000001
& NUC['W']: 00001001
--------------------
            00000001

  NUC['T']: 00001000
& NUC['W']: 00001001
--------------------
            00001000

  NUC['C']: 00000010
& NUC['W']: 00001001
--------------------
            00000000

As you can see, W matches either A or T to give the character it was tested with, and when compared with C, a character it does not represent, it does not match it. The don't care character N will match A, C, G or T/U because it's bit pattern is 00001111. So if we modify our naive algorithm above, we can change the condition of the inner loop so that if a character from X, whether it is A[T|U]CG or a degenerate character, when ANDed with a character from T results in the same character from T, then it is a match.

For ease, I create an array P to hold the binary representation of the characters in X as defined in the NUC array.

P = []
for j = 0 to m - 1:
    P.append(NUC[X[j]])

Then I can easily use it in the modified naive method, which this time supports degenerate and don't care characters:

i = 0
while i <= n - m:
    j = 0
    while j < m and (NUC[T[i+j]] & P[j]) == NUC[T[i+j]]:
        j += 1
        if j == m:
            print "Match found at position: ", i
    i += 1

So there it is. I searched with X = "CCWGG" and it found CCAGG in T. If I had CCTGG in T it would have found it as well. Also, if I didn't care about the value of one or two of the letters, then I could have used N in the pattern and it would have found the pattern I was looking for.

The naive method is limited to searching at O(nm) time complexity in the worst case, but if this technique is incorporated into a better exact string matching algorithm, it will be able to do it much quicker and with the added benefit of finding patterns using a degenerate character set which supports don't cares.

ISO-639-2 list of languages and their directions

I collated a list of the languages (ISO-639-2) and the direction they are most commonly written in and provided it as a TXT or CSV file for download here: https://github.com/webmasterar/languagesWithDirections

Examples:

ara,,ar,Arabic,arabe,rtl
chi,zho,zh,Chinese,chinois,ttb
eng,,en,English,anglais,ltr

This is mostly useful with the correctly displaying text on websites using the HTML bdo tag or CSS styles direction and/or writing-mode.

Latitude and Longitude to Timezone PHP script

I converted a Java code I found online that was really useful to PHP so anyone can use it now. It basically returns the timezone of a location given its latitude and longitude. Check out the example:


<?php
include 'TimezoneMapper.php';
$london = ['lat' => 51.504380, 'lng' => -0.127727];
echo TimezoneMapper::latLngToTimezoneString($london['lat'], $london['lng']);
//outputs: Europe/London

You can get the code here: https://github.com/webmasterar/TimezoneMapperPHP

A quick guide to compiling C/C++ code and writing Makefiles (Part 1: Compiling)

This is the first of a 3 series post. This post will be about compiling. The second post will be about Makefiles. The third post will be about compiling libraries.

C/C++ coding is difficult enough as it is to code and debug without having to spend more time figuring out how write a good Makefile to compile the damn code and fix the Makefile bugs! In this post I will talk about compiling using the command line tools. But usually you don't write these commands manually on the console to build a project and be able to run the application. You write the commands in the Makefile and then you just build the project using the Makefile.

A Makefile is used to store the commands needed to build or clean a project. The user may download your project and just has to write the make command on the console to build the project to create the executable binary file that will run on their machine. In Windows you get given binaries like virus.exe with the .exe executable extension that just works on pretty much all Windows computers as long as there are no architectural requirements (64bit versions of software will not run on 32bit systems). But this is not the case for Unix/Linux. With all the different distros available and their internal complexities, it is not likely that a binary executable built on FreeBSD will also run on Ubuntu. Also there is no quick way to state a file is executable because .exe is not used to denote an executable file in the Linux world. So, for Linux mainly, the source code is often packaged so you can dowload and extract it, and you are expected to compile the code using the make command to create the executable for your machine in order to run it.

Which tool you compile your code with varies by operating system and there is a choice of compilers for each system anyway. On Windows I use MingW and on Linux I use GCC. There are often different compilers for C or C++ code too. It should be noted that C++ compilers are capable of compiling C code but sometimes complain about things that are not acceptable anymore in C++. I recommend compiling C code with a C compiler and C++ code with a C++ compiler. You can mix C and C++ code, but if you do, be prepared to edit your C code.

Anyway, compiling. The C compiler tool in Linux is gcc and the C++ compiler is g++. If you use MinGW to build a project on Windows, it will know what you are trying to compile your code with, so it's cool to carry on using gcc and g++. There are two sequential steps to compiling.

1. (dependent on the project) Make sure you have linked in any required library files (.a/.so) and their header files (.h).
2. Compile the source code (.c/.cpp) files into object(.o) files.
3. Combine the object files into one executable file.

Let's start with step two first because I will write another article about compiling and linking libraries at a later time. The following command will compile the code into object files for a hypothetical project consisting of two source files and the header file associated with the important program: exampleProg.c, reusableCmd.c and reusableCmd.h.

gcc -I . -c exampleProg.c reusableCmd.c

This will create the files exampleProg.o and reusableCmd.o. The flag -I . means find the header files required to compile this program are in the current directory (.). If we needed to include other header files in another directory, we would add an additional I flag like so: -I ./../dir2/libX.

Then we have to compile those files into an executable file. On Linux we don't give it an extension usually, but on Windows we specify it ends with .exe. To compile our files on Linux we do:

gcc -I . -o exampleProg exampleProg.o reusableCmd.o

And there you have it. The code will produce the executable file exampleProg, which is ready to run. The C++ compiler works in much the same way too. It is best to add some additional flags to optimize the executable file or improve debugging. Here are some flags I use a lot in my projects:

-Wall = Show all code warnings
-lm = Load the Math library if you #include <math .h%gt;- you need to add this flag explicitly if you include this header!
-funroll-loops = optimize loops
-O3 = Use maximum processor and memory optimizations
-fomit-frame-pointer = Omit some of the debugging pointers included in the executable - use this only once your code works perfectly
-msse4.2 = Use SSE4.2 machine instructions where available to improve performance - this is especially useful on modern Intel processors.
-std=c99 = This sets the expected standard of the compilation. You can also use an even newer standard such as c11 if your compiler supports it.


In the next post I will talk about Makefiles.

Counting Sort (index-based sort) reinvented

So the past weekend I independently came up with what I would later find out was the counting sort and I went as far as implementing it before I found out it already existed. Awkward. I initially called it the quick dictionary sort and afterwards renamed it to the quick tally sort (qtsort) because a "dictionary" is a list of words and not individual characters and numbers which is what this algorithm is used for.

The counting sort is a linear time O(n) sorting algorithm because it does not work by an element value comparison technique, which has a proven best worst case complexity of O(n log n). Many sorting algorithms are O(n²) in the worst case, but satisfyingly the counting sort is linear. That's a revelation when it comes to sorting because it is so fast! The best thing about the counting sort is that it is so elegantly simple and can take up little extra memory too. The space complexity is also very small if the range of the elements is small O(n_max - n_min).

Let's go through an example with a character array to clarify. Imagine a string t made of a set of characters: t = {A,D,B,E,A}. We know that the ascii value of character A is 65, B is 66 and the rest of the alphabet is ordered sequentionally in the same way. We can create an array in the range of characters in t to hold all the unique values. The character A is the smallest value (65) in t and the character E is the largest (69).

We create an array T, which stands for tally, and it should hold the frequency each character occurs in text t. As arrays start at index 0, we can make any character in t, A through D, be represented by an index of our array by subtracting the minimum value from each characters. Also, our array T, will be of size = t_max - t_min + 1, or 5 = 69 - 65 + 1. To let the letter A be represented as the first index (0) we do: index = A - t_min. T looks like this when we initialize it and we haven't gone through the text yet:

 idx|val
+---+---+
| 0 | 0 |
| 1 | 0 |
| 2 | 0 |
| 3 | 0 |
| 4 | 0 |
+---+---+

Notice that there is position 2 in the array even when we don't have the character C in t. We go through each character in t and each time we increment the value at the calculated index of that character, so eventually our T array looks like this. We have two A's, and one each of B, D and E.

 idx|val
+---+---+
| 0 | 2 |
| 1 | 1 |
| 2 | 0 |
| 3 | 1 |
| 4 | 1 |
+---+---+

The T array is filled out using a simple loop which goes through each character of t and increments the counter, so by the end of the loop it has made a tally of how often each character occurs. This action is represented by this bit of psuedo code:

foreach char in t:
    T[char - t_min] = T[char - t_min] + 1

Now that we have filled the tally array, we can go through it in sequential order and wherever the value is 1 or more, we print that many of the character. This can be implemented with a simple for loop through the range of the items in T. With each index in the array, we work out which character it really represents by adding t_min back to it, so with the first element, c = idx + t_min, A= 65 = 0 + 65; and with the next index, B = 66 = 1 + 65. If the value at an index in T is 0, meaning no characters exist for it in t, then we skip it and go to the next index. Here is the psuedo code:

foreach idx,val in T:
    if val > 0:
        c = idx + t_min
        while val > 0:
            print c
            val = val - 1

This simple bit of code would then print out t in the correct order: AABDE. The algorithm runs linear in time O(n), but is limited by the range of the values, so if you had an huge range of numbers it would occupy a lot of memory for T and would need to go through all of T in order to print it out. Therefore the Counting sort is best applied to long strings or arrays of numbers within a small range of values. It is unfortunately unsuitable for use with other types of data such as floats or words. It's role lies in sorting lists of numbers, such as table indexes, and characters very quickly and using little memory. It is better to use it than other sorts for certain datasets and will yield a result faster than O(n log n) time.

You can check out a working example on github written in C and works only with characters currently, but can be modified with ease to work with larger integers. You don't need to specify t_min or t_max because the program works that out in linear time too! Here it is: https://github.com/webmasterar/qtsort

Understanding GROUP BY and when to use it (SQL)

This is a little guide I wrote for second year computer science students at KCL doing the database module. But anyone is welcome to benefit from this.

---

Imagine we have a list of students and each of these students attend a single college. We create the students table and colleges tables containing data like:

Table:Students
+---+----------+-----+
| 1 | Dave     | KCL |
| 2 | Gemma    | KCL |
| 3 | Ali      | UCL |
| 4 | Adam     | CWC |
| 5 | Claudia  | UCL |
+---+----------------+

Table:Colleges
+-----+------------------------------+
| UCL | University College London    |
| KCL | King's College London        |
| CWC | City of Westminister College |
| WLC | West London College          |
+-----+------------------------------+

Next, imagine the dean of KCL asks you: How many students attend each college?
So now you have to flex your SQL skills and you want to give him a list of college codes and the number of students attending them. Something like this:

+-----+---+
| WLC | 0 |
| CWC | 1 |
| KCL | 2 |
| UCL | 2 |
+-----+---+

But how do we achieve this? Let's give it a try and go through some common mistakes. First, let's do it the naive (newb) way:

SELECT
Colleges.code, COUNT(*) AS num_students
FROM
Colleges
JOIN
Students ON Colleges.code = Students.college_code;
Result:
+-----+---+
| CWC | 5 |
+-----+---+

Not what we wanted, right? Why does this happen? That's right, COUNT() is an aggregation function and it returns one row containing the count of all the students, and along with it the first row of the college code.

You realise you need to use GROUP BY:

SELECT
Colleges.code, COUNT(*) AS num_students
FROM
Colleges
JOIN
Students ON Colleges.code = Students.college_code
GROUP BY
Colleges.code;

+-----+---+
| CWC | 1 |
| KCL | 2 |
| UCL | 2 |
+-----+---+

Observe that because we use an INNER JOIN, WLC is not in the list. How do we correct this? That's right, you use a LEFT JOIN so all the colleges are listed.

So that's got us the correct result. Now, you think to yourself, what if he had asked me for a list of colleges with 2 or more students? So you consider using a condition in the WHERE clause:

SELECT
Colleges.code, COUNT(*) AS num_students
FROM
Colleges
JOIN
Students ON Colleges.code = Students.college_code
WHERE
num_students >= 2
GROUP BY
Colleges.code;

Result:
ERROR: Unknown column 'num_students' in 'where clause'.

OK. So we try again with COUNT(*) >= 2.

Result:
ERROR: Invalid use of group function.

So what is the correct solution?! You realise that you need to apply a condition that is related to the aggregated resultset that is returned by the query.
When you use aggregate functions like COUNT and GROUP BY, you would need to use "HAVING" like so:

SELECT
Colleges.code, COUNT(*) AS num_students
FROM
Colleges
JOIN
Students ON Colleges.code = Students.college_code
GROUP BY
Colleges.code
HAVING
num_students >= 2;

+-----+---+
| KCL | 2 |
| UCL | 2 |
+-----+---+

Just what we wanted here!

Anyway, you get back to our original task - get a list of colleges and the number of students that attend them. You know how to answer it but you want to list the colleges ordered by num_students.
You use LEFT JOIN and ORDER BY to achieve this using the query:

SELECT
Colleges.code, COUNT(*) AS num_students
FROM
Colleges
LEFT JOIN
Students ON Colleges.code = Students.college_code
GROUP BY
Colleges.code
ORDER BY
num_students DESC, code ASC;

And this gives you the output that will impress the dean, impress all who know you and establish you as an SQL guru in the department.

+-----+---+
| KCL | 2 |
| UCL | 2 |
| CWC | 1 |
| WLC | 0 |
+-----+---+

Sampling k-combinations with linear memory and time - the bitmask/shift trick

k-combinations are like permutations but what they achieve is a string of length n made up of an alphabet of sigma (s) characters. Imagine we are dealing with RNA (alphabet = {a, u, t, g}) and we want the different combinations of triplet codons (n = 3), we know that there are sⁿ or 4³ = 64 combinations: (aaa, aau, aat, aag, uaa, uau, ..., etc). While we could generate the combinations and save them all to an array, we might come across a problem if we want a longer string. Let's say we increase n to 30, the combinations of 4³⁰ is 1152921504606846976 combinations. That's significantly more than 64, and we would struggle to generate and store that many combinations in memory. Basically, the size of our array grows exponentially with the increasing n. The longer the string you need, the more time and memory is required to produce and store all the combinations.

What I'm going to briefly explain is a way to get any combination without an expensive algorithm to generate all the combinations in order to store it in main memory or on disk. There is no need for a large array because our combinations are generated on the fly. The time complexity of getting any particular permutation is linear so it is useful for sampling or getting a small set of combinations without having to generate all of them.

Let's say we want combinations of length n:

int n = 3;

First we stick the alphabet in an array:

int s = 4;
char alphabet[s] = {'a', 'u', 'c', 'g'};

Then we work out how many bits are needed to represent the alphabet characters by rounding-up the square-root of the alphabet size:

int bits = ceil(√s); //sqrt(4) = 2

The number of combinations is of course sⁿ:

int numCombinations = pow(s, n); //64

So imagine we want to print out all the combinations, all we have to do is to loop numCombinations times and decode the index of the loop (i) into a combination using a bitmask:

unsigned long long int i, j, c;
unsigned long long int mask = (1 << bits) - 1; //00000011
for (i = 0; i < numCombinations; i++) {
    c = i;
    for (j = 0; j < n; j++) {
        printf("%c", alphabet[ (mask & c) ]);
        c = c >> bits;
    }
    printf("\n");
}

Basically, if you wanted a single random combination, just pick a random number between 0 and 63 inclusive. Create the bitmask and for each character in n shuffle the bit along a set number of bits.

Imagine we chose the number 57 at random. This is 111001 in binary. We can read two bits at a time by doing a bitwise AND using the mask (000011) and then we do a bitwise right-shift two positions:

000011 & 111001 = 000001 //index 1 (u)
shift c right 2 bits.
000011 & 001110 = 000010 //index 2 (c)
shift c right 2 bits.
000011 & 000011 = 000011 //index 3 (g)

When we print out the characters from our alphabet at index 1, 2 and 3 and together we get: ucg.

This technique is really straightforward while (bits * n) is less than or equal to the number of bits in a computer word (w). That is the immediate technical limitation. If you use a unsigned long long int to store your combinations then you can use up to 64 bits on a 64 bit computer. If you need to use more than that then it will get a little more complicated as you juggle bits across multiple computer words.

PIP installer not working (not installing Django and other packages)

I was trying to install Django using the PIP installer today and I came across this error:

~$ sudo pip install -U django
Searching for install
Reading https://pypi.python.org/simple/install/
Couldn't find index page for 'install' (maybe misspelled?)
Scanning index of all packages (this may take a while)
Reading https://pypi.python.org/simple/
No local packages or download links found for install
error: Could not find suitable distribution for Requirement.parse('install')

I found out the way to fix this problem was to reinstall pip using the easy installer and then update pip.

~$ sudo easy_install pip
~$ sudo pip install --upgrade pip

And this allowed me to install Django successfully.

An online Dynamic Programming matrix viewer

I decided to create a single page application that generates and displays the Dynamic Programming matrix of two given strings. This is for any bioinformatician or computer scientist who needs a quick and easy DP tool they can use online and be able to view the matrix that is output: DPMatrix.

It supports different algorithms - Pattern matching, Global Alignment (e.g. Needleman-Wunsch) and Local Alignment (e.g. Smith-Waterman). And it supports two DP models - Edit distance and Hamming distance. And it allows you to enter different penalty values.