Difference: JaysJournal (1 vs. 73)

09/01/2010

• Init() will now be Build(sampling rates ...), and take in paramaeters for sampling rates, with defaults set. The method will build the index from scratch.
• Load(file) will be Load(file, sampling rates ...), works the same way as Build, but it loads from a file instead of building from scratch. If sampling rates aren't provided, Load will use the sampling rates specified in the save file; otherwise, it will override the sampling rates, if it can (sampling rates should be a multiple of the ones in the save, and be larger).
• Save(file) will pretty much do what it's always done, but make sure to save the structure in such a way so that the samples are "skippable"
• I'll have SetXSamplingRate methods at some point that will allow locate and jump table parameters to be settable. Setting it to 0 will unload everything, while setting it to a different sampling rate will reload the structure. These methods will also have a file that can be passed in (optionally), which allows reloading of the structures just from a file.
• Note that the rank structures (StaticBitSequence, StaticDNASequence) won't be dynamically changeable, since those require the BWT to rebuild (I could rebuild from a file, but it's a lot of trouble!)

08/23/2010

• Have options to make building the locate table, jump table, etc. optional
• Add jumptable support
• Add text extraction support (at least be able to extract from a packed representation)
• Add a BlockArray class

Currently, I've got the jumptable support done and optional build requirements. I'm also working on some unit tests for the FM index. Finally, the bug from the last post (Rank0 sometimes fails on blocks with '$') has been fixed.

To do:

• Finish up unit tests
• Add text extraction support from a packed string (to do after block array has been made)
• Add a BlockArray class

08/18/10

The FM Index is mostly done, but without multiple-reference support. Chris says this should be an easy fix and can be updated later. Running a bunch of tests, the FM index seems quite fast when doing Locate (roughly 4x faster than bowtie). However, these results aren't the most reliable because there's no post-processing going into it, nor is there quality support. Some interesting points:

• When doing Locate, the sampling rate of the locate structure (to store positions) matters quite a bit, so this should be as small as possible.
• The StaticDNASequence sampling rate has some leeway in the beginning (in my tests, from 64-256 bases were okay), but starts affecting the time significantly after that.
• The StaticBitSequence sampling rate has pretty much no effect, even up to sampling rates of 1024! Unfortunately, changing this sampling rate also affects memory the least.
• Using SSE4.2 (hardware popcount) doesn't have too much of an effect at low StaticDNASequence sampling rates (expected); at a sampling rate of 64 bases, the SSE-enabled version was about 4% faster, but goes up to 20% higher at a sampling rate of 1024 bases.

Lately, I've been working with Daniel to integrate the new FM index into NGSA, which has taken most of the past two days. Right now, we have the BowtieMapper working (albeit with a lot of quick hacks, since multi-references isn't done yet), which maps reads with up to 3 mismatches and uses a jump table. Daniel is also going to work on the BWAMapper, which currently has some problems with unsigned integers and -1's. I'm also going to reorganize the structure of the NGSA library and make a better utilities file.

Finally, I ran some preliminary benchmarks between the new saligner and bowtie. The following were run on NC_010473.fasta, minus the 'R' and the 'Y' character in the genome, using the e_coli_2622382_1.fq readset. Tests were done on skorpios in /var/tmp/, using 64-bit versions of both aligners:

Note the times might not match so well because I used time for bowtie, but the built-in timing function for saligner, since I haven't gotten the saving/loading integrated yet. Also, I didn't use valgrind --tool=massif to profile saligner, because there are some full-text references being kept in memory somewhere, which is really raising the memory usage (I'll have to find and clear those later). The memory reported above is from the FM Index's FMIndex::MemoryUsage() function, which only reports the theoretical usage, and I kept the memory usage for saligner on the low side, to account for any extra memory that may be used for other things.

To do:

• Reorganize directory structure
• Fix a bug in the FM Index/=StaticDNASequence=, where Rank0 gives the wrong output in blocks with '$'
• Implement a jump table within FM Index
• Maybe implement multiple sequence support (Chris will get back to me on this)

08/10/10

08/11/10

To do:

• Benchmark more thoroughly
• Test more thoroughly
• Implement a feature to load/save only the BWT string, instead of the whole index, which could change based on different memory usage profiles.
• Get ready for integration?

08/11/10

I also did a few preliminary benchmarks, and the times were not great on the Locate function. I think this might be because we don't implement Locate the "proper" way, which guarantees a successful locate within sampling rate number of queries. Following Chris' suggestion, I tried graphing the number of backtracks it takes before a successful location is found on a random readset, and here are the results:

The sequence length was 65,751,813 bases long and consisted of the E. coli genome (4.6 Mbp) repeated multiple times. Reads were randomly generated 10-base sequences, and only the first 10 matches were "located". This was run at a sampling rate of 64 bases, and the average distance to locate was 341.68, which isn't very good. The x axis of the graph represents the distance it takes to locate, and the y axis is the number of aligments with that distance (it's logged). There were a total

• Reference length = 65,751,813 bp (consisting of the E. coli genome repeated multiple times)
• Randomly generated 10-base sequences, only the first 10 alignments were "located"
• Run at a sampling rate of 64 bases for the locate structure
• Average distance = 341.68, which isn't very good
• x-axis is the distance it takes to locate, y-axis is the number of alignments with that distance, semi-logged
• 855,150 locates performed in total

To do:

• Make a "proper" locate structure and compare.

08/06/10

Chris suggeted in an email that we prefetch ep ranks while doing both counts and locates, which might make it faster. I'll be discussing this with him later.

Finally, I finished the FMIndex::Locate function today. It works as expected, and I've tested it pretty thoroughly. Locate requires a bunch of LF mappings of each position in the BWT string to travel "backwards" until a position is reached whose position has been indexed. However, I can't actually get the specific character at the given position in the BWT to do the LF mapping. The FM index paper suggests doing an LF on both position i and i-1 for each character in the alphabet, since the only one that will differ in LF value is the character I'm currently on. I'm also thinking I could just implement a select function in my rank structure, which shouldn't be too hard.

To do:

• Benchmarks! I might be able to integrate the index into saligner without too much trouble, so we can get some pretty accurate comparisons?
• Implement saving/loading of the whole index. Currently, I support saving/loading of the rank structure, and the FM index is just a few more data structures so it shouldn't be too bad.

• reference = e coli genome NC_010473, minus two characters (there was one Y and one R...) (length = 4,686,136 bases)

• Finish the $ data structure implementation in the StaticDNASequence

08/04/10

Also did a few benchmarks using the new 2-bit "enhanced" version and repeated some benchmarks from last time: Test parameters:

• read length = 35
• reference = e coli genome NC_010473, minus two characters (there was one Y and one R...) (length = 4,686,136 bases)
• number of reads = 100,000, generated randomly
  • this randomly generated set was copied 1000 times, for a total of 100,000,000 calls to count.
• random seed = 1234
• times are reported in milliseconds wall time.
• Build version is 64-bit, Release mode (-O3)
• Memory usage is just the memory used by the rank structure
• There was one $ in the sequence, no N's

To do:

• Start doing a Locate method.
• Maybe do some code cleanup

07/29/10

Chris and I have also decided on a naming scheme, as follows:

• 1-bit = BitSequence
• 2-bit = DNASequence
• 3-bit = EDNASequence (for "Enhanced DNA sequence")
• Rank structures have "Static" prepended
• Packed strings don't have "Static"

I also benchmarked the various rank structures, shown below. These benchmarks were done on skorpios in the /var/tmp/ folder. Times are reported over 25,000,000 calls to rank on a 512,000,000 base sequence with varied sampling rates. The graph reports times on the y axis (milliseconds) vs memory usage on the x axis (bytes). As can be seen, the StaticEDNASequence is much slower than StaticDNASequence given the same memory usage, but one good thing is that they're not too much slower given the same sampling rate.

Today, I successfully implemented a count function using the StaticEDNASequence. Also tried to work out a nice hack for getting rid of the end markers when using the StaticDNASequence, but it doesn't seem possible. So right now, I'm going to build in another data structure into the StaticDNASequence that keeps track of all the $'s and subtracts these from the A ranks. It might get complicated, but hopefully the speed difference won't be too big.

To do:

• Build the new data structure into StaticDNASequence
• Benchmark the FM index count method using both 3-bit ranks and 2-bit ranks, to get a good idea of what speed we may be able to achieve.

07/30/10

Finally, I did some benchmarks on the count function. Test parameters are:

• read length = 35
• reference = e coli genome NC_010473, minus two characters (there was one Y and one R...) (length = 4,753,080 bases)
• number of reads = 10,000,000, generated randomly
  • this randomly generated set was copied 10 times, for a total of 100,000,000 calls to count.
• random seed = 1234
• times are reported in milliseconds wall time.
• Build version is 64-bit, Release mode (-O3)
• Memory usage is just the memory used by the rank structure

For a (very) rough comparison, saligner, using readaligner's FM index, takes 64.461 seconds to do an exact mapping for 2.7 million reads (each 35 bases in length). This means it would take roughly 2,387 seconds to align 100,000,000 reads. However, this is not a very good comparison because that one's doing a locate function which is more complex, has to read from a file and write to an output file, and runs at 32 bits. Also, the rank structure might use a different sampling rate, although I believe it uses a sampling rate of 256, which is why I included that test in. So, this means we have quite a bit of room to work with!

Another note is that for the 3-bit structure, it requires a sampling rate of 2048 to be roughly comparable in memory usage. However, at this rate, it's much slower than the 2-bit structure (roughly 5-6 times slower).

To do:

• Finish the N data structure implementation in the StaticDNASequence

07/27/10

The change from iterators was actually surprisingly easy because I had the rank blocks build separately from the count blocks (even though they were part of the same array). Hopefully, this will make everything more flexible. As a bonus, I also cleaned up the rank structure classes quite a bit, removing extraneous methods and updating the documentation.

The 1-bit rank structure wasn't too bad to code, since it's so simple. I had to get rid of the interleaved count blocks because it just wasn't practical, but the memory usage is quite good because I only have to store counts for 1's (Rank0(position) = position - Rank1(position)). I haven't really done exhaustive testing for this class, so it might still have lots of bugs.

Finally, no speed tests yet, but I should have them up by tomorrow.

To do:

• Test all the new classes for speed and bugs.

07/26/10

1. Build a new rank structure that supports 3 bits (for 8 different characters)
2. "Ignore" the N's and $'s, but note their positions. After alignment, do a post-processing step where we filter out results that cross an N or $.

The first option would be easier to do, and would provide support for wildcard matching (later). The second option is much faster and memory-efficient. For the same memory footprint, my 2-bit rank structure is about twice as fast as the 3-bit structure.

Today, I started building the 3-bit rank structure, since Chris tells me that we'll need it sometime in the future anyway and the 2-bit rank structure is still fresh in my mind. It looks like the structure works pretty much as desired, though I haven't done exhaustive bug tests yet. Also, I finished the BWT string builder integration (though not for multiple chromosomes), and I've successfully built a BWT string out of chromosome 10 of the human genomes, ran the resultant string through my 2-bit rank structure and got good ranks out of it, at least in the first 10 or 15 bases.

While typing this, I realized I forgot to change the save/load methods for the 3-bit rank structure, so I'll have to get on that tomorrow. Another thing I should do is to build a 1-bit ('regular') rank structure. The problem with this structure is that I can't interleave the lower-level count blocks, since I only need one 16-bit count per sample, as opposed to four different counts. So, it would be a big waste of memory if they were interleaved. So, I'm just going to have the counts separate from the rank structure.

Finally, for the initialization functions, I'm really considering using iterators instead of requiring an operator[] to be present. char * is already an iterator in and of itself, and it's just a matter of implementing an iterator for the packed string. This way, it's much more flexible, and I might later be able to make a "streaming" class that loads from a file a little bit at a time.

To do:

• Fix the save/load function of the 3-bit rank structure
• Maybe add iterator support
• Make a one-bit rank structure

07/23/10

Initially, I had concerns over a FIXME type comment in the actual RLCSA class (which RLCSABuilder builds), which said the program crashes when you enter in a sequence > 2 GB. Looking through the code more carefully though, I think the RLCSABuilder class builds a separate RLCSA for every chromosome (each of which is < 2 GB), and then merges them.

Also, I don't think using packed strings with the RLCSABuilder is possible, because it requires that the sequence it's given has no '\0's (it uses '\0's as the end of sequence marker). Also, it uses char arrays, which is pretty much hardcoded into the code, so using anything smaller isn't really possible without some major changes.

So my plan now is to just build a BWT string, pack it, run it through my rank structure, and then see how that goes from there.

To do:

• Make sure the BWT builder works!

07/21/10

Also ran some tests on antiparos using the 32-bit and 64-bit rank structures. Since antiparos's CPU is actually 32-bit, there was a huge difference in speed between the two implementations. Of course, the 32-bit version was faster, ranging from 20% to as much as 70% faster (the difference gets wider with higher sampling rates). Given this, I will be keeping both implementations, which I've also managed to merge. The only thing I don't like is that I have to have a selective compilation in each of the rank functions, which makes the code look a little unwieldly. Fortunately, everything else was easy to merge by defining a constant.

I looked into Daniel's packed string class as well. I haven't actually started adding/modifying it yet, because I figure it'll be easier when I start merging the new index in; it's going to require a little CMake tweaking as well, because I want to give it a 64-bit implementation as well. Some changes I would like to make:

• Add a 64-bit implementation and have it selectively compile (currently, the data_ member variable is a uint32_t. This should be uint64_t in 64-bit.
• Add an operator[] method that returns the ith base; this shouldn't be too bad.
• Build the string a little better, using a hash table instead of a switch/case.

Finally, I started looking into BWT builders. The BWT builder is in the incbwt folder (kind of obvious now that I think about it), so it's just a matter of extracting that folder and trimming down on some classes. I think the rotational index is also in that folder, which is why there are so many files in it currently. Chris also suggested an implementation he had from Karkkainen (blockwise suffix sorting), which is what bowtie uses (I think?). This version is able to build the BWT string using much less memory, but the implementation may be buggy, though it may just be that it uses signed integers (too small). As well, Chris isn't sure whether the code is open-source or not, so he's looking into that, too.

To do:

• Try fixing the Karkkainen implementation of the BWT builder.
• If that fails, extract readaligner's BWT builder.

07/20/10

• Rank structure is nearly finished, just a few minor changes:
  • Should change sampling rate to reflect bases, not blocks of bases
  • Put in packed string support (use templates)
  • When taking in a regular string, assume I'm getting values of 0-3, instead of letters. Make a translate function to translate between letters and numbers.
  • Implement a static memory function that outputs how much memory a given rank structure requires (in preparation for the "memory budget")
  • Pull out some functions and hashes into a utils file
  • Try combining 32-bit and 64-bit implementations into one file to minimize redundant code
  • Look into using a "binary" bitstring instead of two bitstrings side-by-side for cleaner code.
• For the FM index, I'll be:
  • Taking a look at readaligner's FM index building function and pulling it out, since we won't be touching on that (building is quite complicated as well)
  • Also implement a memory budget feature for the index (when getting positions)
  • Start implementing the FM index itself, assuming 4 bases with no N's.

Also, I revisited the 32-bit code and was actually wondering if it was significantly faster than the 64-bit rank implementation on 32-bit systems. Chris suggested I run some tests on antiparos, which probably has a 32-bit CPU. More on that later.

To do:

• Merge the two implementations
• Implement a memory function
• Start looking at Daniel's packed string implementation. Most importantly, add an operator[] function.

07/19/10

I've started reviewing how the FM index works, in preparation for building it. More details on that tomorrow, though.

To do:

• Discuss with Chris about what to do next, get feedback, etc.

07/16/10

Also did a little reading on CUDA and using GPGPUs. I found one paper detailing a GPGPU-based Smith-Waterman implementation, which I briefly glanced through and also went through the =CUDA= developer guide.

Finally, had a long discussion with Daniel about his new BowtieK1Mapper and multithreading. Daniel's new implementations seem to modify the index a little during the mapping process (he "flips" it), which doesn't go well with multithreading, because I assume the mapping function is constant and can be run independently. Long story short, we decided to hold off on modifications to the class, since we'll be implementing a new FM index soon, which will no doubt involve a lot of changes (the least of which is probably getting rid of the NGS_FMIndex class, which is just acting as a middleman for readaligner's FMIndex). I've also decided to hold off on committing my merge to the dev trunk for multithreading, because this would introduce compile errors with BowtieK1Mapper in its current state.

To do:

• More reading.
• Talk to Chris about the FM index.

07/15/10

Yesterday, I managed to actually finish up my two-level implementation (turns out there were a few bugs in certain cases, due to +1/-1 errors). I also implemented save/load functionality, so that the structure itself can be saved. To use, just call RankStructure::Save on a RankStructure that has already been initialized (by calling RankStructure::Init). Loading works the same way, just call RankStructure::Load with a filename supplied. Note that when loading a RankStructure, RankStructure::Init doesn't have to be called first. I'm not sure whether I should be taking in a file name as a string or whether I should just take in a stream pointer; readaligner's rank structure seems to just take in a pointer, so I might switch to that. Only problem is, I have to assume that the file's been opened as binary.

Today, I managed to finish my 32-bit implementation, which runs much faster on my machine (operating system is 32 bits). I also did a bunch of rigorous tests on both the 32-bit and 64-bit versions of the code, trying to fix all the bugs I could find. Both implementations are now pretty much valgrind-error free, and I'm pretty sure they both work correctly, too. One very annoying bug I had was during memory allocation, where I had an expression, size << 3, where size was a 32-bit unsigned integer. I later corrected this to a more generalized size * BITS/8, which ended up blowing up at sizes of around 2,048,000,000 characters. Well, it turns out the multiplication step is done first, which, when done, makes it too big for a 32-bit integer and rolls over. Sigh...

I also did some tests on both implementations. Initially, I found the 32-bit version to be faster on skorpios, which was very strange. However, I realized that this was because of the way I chose to represent the sampling rate (I chose to represent it in blocks, so a "sampling factor" of 1 is actually a sampling rate of either 64 or 32, depending on the implementation). So, the 64-bit implementation actually had twice the sampling rate given the same "sampling factor", which was what made it slower. Accounting for this, the 64-bit implementation is (rightly) faster.

Finally, I used the save functionality to test the memory size, just to make sure everything was going right. Using my two-level 64-bit implementation, I put in a 1,024,000,000 character sequence with a sampling factor of 1 (so sampling rate = 64) and saved it to a file. The file turned out to be 366.7 MB, which should be roughly the memory footprint of the structure. At 32 bits, with the same parameters (note sampling factor = 1, so sampling rate = 32), I get a file size of 488.8 MB. At sampling factor = 2 (sampling rate = 64), I get 366.7 MB again, which is exactly the same as the 64-bit implementation.

To do:

• Not really sure, actually. I'm going to finalize the rank structure, then either...
• Do some work on the Local Aligner, like optimize the FullAlignment methods
• Read some papers/articles on GPGPUs and CUDA and also RNA folding pathways to decide a bit better on what I would do for my thesis.

07/13/10

So now, any sample factor that is inputted is rounded down to the next smallest power of two. This way, I can use bit operations instead of divisions and multiplications, which are slightly slower. Running a few tests, I see about a 10% performance gain.

To do:

• Start on 32-bit implementation

07/12/10

• There is one "block" per 64 characters.
• The top level will be 32-bit unsigned integers.
• The bottom level will be made of 16-bit unsigned integers.
• There will be one bottom level count every (sampling rate) blocks.
• There will be one top level count every 1023 (the maximum number of blocks that can be counted with 16 bits) blocks

The implementation shouldn't be too bad. I've decided to stick with the interleaved counts implementation for now, where all the counts are interleaved with the sequence blocks. Adding a higher level would just mean adding more blocks, this time with a constant sampling rate of 1023. Ranks would just be the sum of the nearest low-level block and the nearest top-level block. I don't think I will interleave the top level blocks, since they're much too far apart to benefit anyway.

Another optimization I've been thinking about is to limit the sampling rate to powers of 2. That way, when finding the nearest blocks, I can just do bit shifts instead of costly divisions. I'm also thinking about using a sampling rate of 512 for the top-level blocks, so that I can use bit shifts as well. The difference in space usage isn't very large, less than a megabyte (in theory) for a 4-billion base sequence (calculation here, multiplied by 4 for all four bases).

Finally, here is the graph I promised on Friday, compared to sampling rate. Tests were done on skorpios with sequence size of 512,000,000 over 25,000,000 times each and a constant seed:

It's a little hard to tell, but it looks like the count interleaved implementation shows minor improvements in speed (about 5%) up until sampling rate=128, where it starts to get much slower.

To do:

• Finish implementation of 2-level rank structure
• Start on 32-bit implementation

07/10/10

07/09/10

Times are in milliseconds, tested on skorpios with a constant seed, and ranked 75,000,000 times per test. Squares are interleaved, diamonds are 3-base, and triangles are interleaved count implementations. The three different versions all use different popcount functions: blue is the function provided in readaligner, which uses a table; yellow is one provided by Wikipedia, which uses a bit-mangling method; and orange uses the hardware instruction.

In previous tests, I noted that the table version was faster than the bit-mangling version; it's interesting to note that in these implementations, the bit-mangling popcount is actually better. My thoughts are that once everything isn't able to be fit nicely into the cache, then the popcount table may encounter cache misses, which slows it down dramatically. It's also interesting to see that when using the tabular popcount, interleaving the counts is actually the slowest implementation. Again, this might be due to excessive cache misses from having to cache the table, and it's actually what led me to try the bit-mangling popcount one more time. It's nice to see that interleaving the counts has some effect, especially in the hardware popcount case, though.

Finally, I talked with Anne and Chris about a new project, possibly for my honours thesis, today. Chris proposed a GPGPU version of the Smith-Waterman aligner I vectorized earlier, or even GPGPU-capable succinct index. Anne also suggested that I look into the RNA folding pathway research currently going on as well, and I'll probably start reading some papers on that soon.

To do:

• Look into both the GPGPU project and get more familiarized with the RNA folding pathway project
• Start implementing a second level in the rank structure
• Look into a 32-bit implementation

Edit: I realized I forgot to benchmark everything WRT sampling rate, so that's another item on the to-do list.

07/08/10

I also managed to finish an interleaved count implementation, which was much harder. In fact, I pretty much had to do a complete recode of the core build/rank functions. I need to run a few more tests on the implementation, but it should be running fine.

Since skorpios is back now, I'll probably do a bunch of experiments tomorrow. I should also make up a 3-base sampling rate implementation, just to test the cost; again, however, if we end up going with the count-interleaved version, the 3-base version really isn't very practical.

To do:

• Benchmark!
• 3-base implementation (for sampling-rate enabled, interleaved)
• 32-bit implementation, look into two-level rank structure.

07/06/10

• Use larger sequences (>=16,000,000 bases)
• Use random access. More specifically, build an array of random numbers, then use those numbers to call random ranks. This should prevent streaming.
• Try both hardware and software popcount.

Also, I finished implementing count blocks; currently, each 64-bit block has a count block (which is 32 bits), with 4 count blocks per rank block, one for each base.

So here are the results, run on skorpios at -O3, with 40,000,000 rank calls each. Times are wall time:

To do:

• Try out using a different count structure (only store 3/4 bases, and calculate the last base) and see the performance hit.
• Try out interleaving the count structure as well.
• Implement adjustable sampling rates.

07/07/10

Again, times are in wall time, compiled at -O3, single threaded and run with 50,000,000 ranks each.

The 1024M times are again a bit dubious. I fixed the segfault bug, but I think my computer was running out of memory, and it looks like a big jump in time. Regardless, I'm going to test these again on skorpios when it gets back. It doesn't look like there's too much of a performance hit with the 3-base implementation, but if we stick with interleaved counts, there's really no point in going that route, since the 4 32-bit counts fit nicely into 2 64-bit "rank blocks". Using only 3 counts wouldn't do much for memory, since there will still be those 32 bits available.

I also experimented with some different popcount functions, found on the Wikipedia entry. So far, it seems like the current popcount function (which I stole from readaligner) is the faster than the one proposed by Wikipedia. However, the current popcount function relies on a table of values, which might make it slower if there's a lot of cache misses on it, though it outperformed Wikipedia's function even at 1024M bases (don't have the exact values, another thing to benchmark).

Finally, I started implementing the sampling rates; there's currently a bug somewhere, but it shouldn't be too bad to fix.

PS: Before I forget, I realized I did actually lock the whole single end output processing section in my multithreaded driver (I originally thought I just locked the part where it writes to output, whereas now I realized that it locks before that, when it's processing the SamEntry into a string). I should experiment with locking it at a different place (though at this point, it's not that easy, unless I want to put locks in the sam file class itself.)

To do:

• Benchmark above on skorpios.
• Finish sampling rate implementation.
• Get started on 32-bit implementation?

07/05/10

I also did some benchmarks: On skorpios at O3, using a test sequence of 1280000 characters, ranking from position 1230000-1280000:

It's interesting to note that on the parallel implementation, ranking C and G is actually faster than T! This is strange because the calculation for T involves an & operation on the two blocks, whereas C/G involves an & and a ~ operation. The times are more as expected under the sequential implementation, although T ranking should still be faster (it's got fewer operations!). What's more surprising is that the speeds for C/G are faster on parallel!

Using cachegrind, it seems that there are far fewer (about 75%) "D refs" when ranking C/G than A/T in the parallel implementation (turns out "D refs" are main memory references). Specifically, it seems that there are much fewer D writes occurring in the C/G cases (reads are fewer as well, but not as significantly as writes). Furthermore, under the sequential implementation, the D references for all bases are more on par with the A/T rank functions in the parallel implementation. It seems there might be some specific optimization in the parallel C/G case that the compiler is making. I kind of want to investigate this a bit further; ideally, I can somehow end up with the parallel C/G times with the sequential A/T times.

To do:

• Investigate above behaviour a little more.
• Start implementing precomputed blocks.

07/02/10

I also did a few tests with the new calculation method ( ~second & first) & (0xFFFF...FF >> BITS-i ). I tried making an array of length BITS to store the precalculated masks for every i. That way might lead to fewer operations (not really sure on that point). Preliminary tests didn't show much difference, especially on -O3. The strange thing is, I actually compiled a minimal rank function using either of the two calculations as assembly, and at -O3, there wasn't even any difference in the assembly code! Strange, since the two aren't really numerically equivalent, just "popcount equivalent".

To do:

• Explore sequential vs parallel blocks
• Add precalculated blocks.

07/01/10

Also worked on the rank structure a little more. I can now build sequences greater than 64 bases in length (so as big as needed, with memory restrictions of course). I figured out how to align to 16 bytes with posix_memalign. I've almost completed the rank function as well, except I'm just getting a small bug with a boundary condition (I'm pretty sure I'm just missing a +1 somewhere).

I was thinking about the calculation needed for the rank structure some more. Currently, the calculation goes something like this (for example, for a C): (~second & first) >> (BITS-i-1). Chris told me that this takes something like three operations to do, since some operations that are independent can be done in parallel. I figure the following equation would also work for calculation: (~second & first) & (0xFFFF... >> BITS-i) (these two formulas aren't numerically equivalent, but they should give the same popcount result). I'm not sure if this will do anything, but it looks like I might be able to drop the operations to 2, if I can somehow calculate the mask fast. I'm going to have to play around with that a bit more.

To do:

• Fix that bug with the rank.
• Play around with the rank calculation.
• Find out what's going on with multithreading.
• Start implementing precalculated blocks.

06/30/2010

Tested on commit 68bab30017bd18dd03f0689760111f229d595238, running everything from the /var/tmp/ directory without nice. Tested with the single NC_010473.fasta genome, with the 2622382 readset, using regular exact matching (no jump-start):

Same as above, but using OneMismatchMapper:

Notes:

• Didn't go past 5 threads, because apparently skorpios is actually an i5 processor (not an i7), without hyperthreading, so there's not that much point.
• Ratio was calculated comparing user times with the first entry (saligner compiled with multithreading set to OFF)
• There seems to be a slight delay at the end of the alignment, when saligner finishes (all console output is done), but doesn't actually end. This wait comprises 2-3 seconds regardless of thread implementation and might be one factor of the poor ratios (though not much; even when this is taken into account, the difference is minimal, less than 0.05)
• User time seems to make a big jump between 2 threads and 3 threads for both matchers. Maybe this is a clue?

Shifting gears, I've also started to implement the new rank structure, as discussed yesterday. I've completed step 1 and can now build one 64-bit block from a 64-character sequence of A/T/G/C. I'm using Daniel's representation for now, where A=0, C=1, G=2, T=3. It's interesting to note one rank operation is slightly faster than the other 3 rank operations, so it may be advantageous to change the alphabet assignment up to whichever base is the most common? Determining the most common is probably too slow and/or the speedup probably wouldn't be worth it though.

To do:

• Rank structures!
• Figure out why the multi-threaded implementation isn't working as well as hoped on skorpios.

06/29/10

I realized that saligner wasn't running at high priority when I was doing tests yesterday, so doing that made it utilize both CPU cores on my machine. Preliminary tests on my machine shows that running with two threads drops the runtime to something like 51-52% of original speed, so there's not too much overhead from thread switching. Adding more threads does slow it down a little.

Finally, I got a new project from Chris. I'm to build a new rank structure he proposed, which should utilize the cache more efficiently and does away with the wavelet tree. Basically, we will be using two arrays for the popcount, so that we can accurately represent the four bases (as opposed to one array, which would only give us two values). The implementation will involve:

1. Implementing the structure for one block of 64-bits
2. Implementing the structure for multiple blocks of 64-bits (there's two ways to store the blocks, parallel or serial, so I should experiment in this step or step 3.)
3. Have multiple blocks of 64-bits, with counts for every block (there's multiple ways to store counts, so I should try experimenting here)
4. Adjusting sampling rate for every 32-bits, 64, 128, 256, 512.

To do:

• Finish up the multi-threaded implementation: do some code inlining, check everything over, test on skorpios and merge.
• Start on doing the rank structure!

06/26/10

I noticed something strange while running the threaded driver: it seems that when I specify 2 threads, it doesn't use both cores of my CPU, but when I specify >=3 threads, it seems to start utilizing both cores.

The threaded driver does work though; at least, it aligns all the 125k reads fine. I still have to lock up the driver files (single end driver and pair end driver), because I'm still running into a minor race condition when approaching the end of the file. Also, I still have to add the alignment status reporting into the threaded driver (e.g. report how many reads were aligned, not aligned, failed, etc). Finally, I need to verify the output of my threaded aligner, which I'm not really sure how to do; it seems to me that the threaded aligner might print reads out of order, so a simple diff won't suffice...

To do:

• Add the "bells and whistles" to threaded driver.
• Lock up driver classes to avoid race conditions.
• Verify output.
• Benchmark.

06/25/10

Google guidelines

Anyways, I've added options in the CMakeLists.txt file for configuring multithreading, and also SSE support. For threading, I've made a new driver called MultithreadedDriver (I'm going to change this to ThreadedDriver), which takes in another driver. It will then create x number of threads, each calling the driver's Run method. If the user wants no multithreading, then ThreadedDriver just calls the Run method of the driver it gets without making any new threads.

To do:

• Finish up threading.

06/24/10

Also began looking at multi-threading. It doesn't look too hard, and I've verified that the only places requiring locks should be with FastqFile and SamFile. Also, the two mappers, MismatchMapper and GapMapper, which use readaligner's query classes (MismatchQuery and IndelQuery respectively), might need locks when querying, because these seem to have shared variables. This might not be such a good plan because the bulk of the calculation is in the query segment, so it might as well be single threaded. To make matters worse, none of the mapper classes are copyable, so making copies of the mappers for each thread isn't an option. Nonetheless, this shouldn't be too big of a problem, especially if we're building a new index, since we probably won't be using the query classes much anyway.

On another note, I tried integrating the Boost.Thread libraries in, but I really don't know how with CMake. I tried extracting the Thread library using boost's bcp tool, and then adding that library into our build but I got a lot of compile errors; I think the thread libraries are compiled specially with bjam. So I'll have to ask Chris about that tomorrow.

To do:

• Integrate the Boost.Thread library.
• Finish up multi-threading.

06/23/2010

Anyways, some outstanding bugs:

• MutantMapper output isn't quite the same as the previous aligner. I think there are two different bugs:
  • CIGAR strings on some reads add up to 36 on 35 nucleotide reads, which is strange...for example, one string will show up as "33M3S", where the old aligner shows "33M2S".
  • For reverse strand reads, some positions seem to be 1 base less than the old aligner. I'm really not sure why...
• KmerMapper output also isn't quite the same:
  • CIGAR strings for some reads also add up to 36, just like MutantMapper
  • Regardless of forward/reverse strand, some alignment positions are off, but by more than 1. This also screws up the CIGAR string.

The similarities of the above two cases seem to point to something that they have in common going wrong. They both use a LocalAligner as the final step of their alignment, so I'm thinking there might be something wonky going on there. I haven't really touched the LocalAligner yet, because I want to make some changes to it later on anyway, so I might as well put that off. Also, these two mappers won't be used too much for the time being (I don't think), so I think I"ll just put that off for now...

Chris also suggested I take a look at threading the aligner in the meantime, while he's still thinking up something for me to do. He suggested looking into Boost threads, which I did. I don't know...from the arguments I saw, the arguments seem to be more along the lines of "boost threads are easier to use" and "why not use boost threads?", so I think I'll just stick with those. Now I just have to figure out how to add the library to CMake for building. The actual threading looks kind of simple, since I can just have each thread work on a different read. The only shared object the two aligner threads would have are the input/output classes (SamFile and FastqFile) as far as I can see, and I think locking them up won't be too hard.

To do:

• Push my latest changes after doing some final robustness tests (I'm kind of nervous...don't want to break anything!)
• Add multithreading support!

Chris suggested three things from the top of his head:

• Work on building the index (which currently requires something like 12 GB for the human genome) in a memory-efficient fashion
• Build an FM index from scratch, which would require things like a rank structure, etc.
• Work on RNA Seq stuff (intron support, splice graph, etc)

So far, option 2 seems the most interesting to me, but all of them are pretty good!

Anyways, managed to finish merging everything. Mismatches and Exact matching work fine so far. The jumpstart mapper doesn't work (segfaults somewhere), but I think that might be something to do with the jump table not saving to the right place. Daniel also made some bug fixes, so I think I'll have to do some final (minor) merges, and then I'll be done the conversion!

To do:

• Think about some bigger parts of the project I might like to do.
• Finish up the C-string conversion (hopefully, done by tomorrow).

06/22/10

Worked on moving over from C strings still. I managed to fix all the segmentation faults, even with the index-building step. Also fixed most of the bugs that arose from the conversion. Now, the exact aligner returns exactly the same output. The OneMismatchMapper doesn't quite do it, but I'm hoping the merge will take care of that.

Also worked on merging the branches, which might take a bit of time. There's so many new changes!

06/18/10

Finally got my code to compile for the C string conversion. Now it's off to fix bugs and segmentation faults! I managed to clear up all the segmentation faults during the actual alignment process, but the output's still not quite right. My new build seems miss most of the reads that are alignable from an old build. Also, I get a segmentation fault during the index building process, which is really strange because I haven't even touched the code for that. No new code is even called until after the index has been built! Maybe I'll have to run the build through the debugger...

I also spent a large portion of the day experimenting with an idea I had last night. When calculating the rank value for a position i, instead of going to the next smallest pre-calculated position and doing popcounts up until i is reached, I modified the function so it goes to the nearest pre-calculated position, and either goes up or down from there. According to callgrind, popcount only gets called about 60% as much on the 125000 read dataset (170,862,051 vs the original 254,181,103 calls). Unfortunately, this only results in a very small speed increase, something like 3-5%, which is strange because I expected the gain to be a lot more.

Edit: Just noticed Daniel's started to implement a different rank structure, so my change might not even be relevant :(.

Finally, I noticed Daniel's made a lot of changes...not looking forward to the huge merge later...

To do:

• Fix bugs in the new C string version.

06/11/10

I also did a few benchmarks on a bunch of string vs char array functions, and I've developed the following guidelines for any future compatability issues or conversions:

• Converting a string to char array using std::string::c_str() or std::string::data() are negligibly fast, so don't worry about doing it!
• There's no difference as far as I can see between std::string::data() and std::string::c_str() in terms of speed.
• Converting a char array to string (using std::string(const char *) or std::string::operator=(const char *)) is very slow (O(n)), avoid as much as possible.
• Using strncpy with some pointer math and appending a NULL is much faster than string::substr (something like 4-5x faster).
• std::string::operator[] is a bit slower than char array's [] operators as well, probably because the string does bounds checking (not sure on this). Even doing std::string::c_str(), then using the [] operator is faster (though not by much)
• Also, be wary of the string's operator=(const char *) assignment, which sometimes masks a char array --> string conversion. Try being as explicit as possible when setting a string equal to a char array (by using the constructor), to avoid confusion

Finally, I think I might take Chris up on his offer for taking a break to study for my MCAT. I'll see how far I get on the weekend, but I might take Monday off for an extra study day. I guess we'll see after the weekend.

To do:

• Finish the conversion by getting stuff to compile! (I really hope saligner ends up being faster).

06/10/10

#define NGSA_READLEN 128
typedef char SequenceBase;
typedef char QualityBase;
typedef SequenceBase SequenceString[NGSA_READLEN+1];
typedef QualityBase QualityString[NGSA_READLEN+1];

To do:

• More conversion...

Finally, spent the rest of my time conerting strings to char arrays, though I'm still not sure if it's doing anything. Regardless, it should be more future-proof.

There's also one point I'd like to look at in the ExactMapper, though I doubt it has too much of a performance impact. When the total results of an alignment is greater than max_results, it seems we do a Locate on all the results, not only max_results. This should be a pretty easy fix.

To do:

• More converting
• Look at fix mentioned above.

06/09/10

I also benchmarked my changes from yesterday (swapping out the Fastq read function). It seems saligner is now faster than readaligner (at least in exact matching). Anyways, here are the times:

The mismatch:exact ratio is still a bit off from bowtie's and readaligner's, but that might also be because mismatches still unecessarily go through the CIGAR generator (which I should fix sometime).

After the meeting, I did a bit of tinkering with the FASTQ reader. After some discussion with Daniel, we decided it might be better to trade robustness for speed in the reader, as most FASTQ output these days are more standardized (for example, no wrap-around for sequence or quality strings). So by taking all those checks out, I was able to get another roughly 10-15% increase in speed. I also started doing some work on converting our C++-style strings to C-strings. I started with the FASTQ reader again, but it's going to end up being a huge change that'll effect almost all of our classes. At the same time, I think I should also start looking into using regular arrays instead of vectors where applicable.

To do:

• Pretty much refactoring (converting C++ strings to C strings).

06/08/10

It's a lot closer to readaligner now. I ended up changing all the osstreams and ofstreams to sprintf and fprintf, which made it quite a bit faster. I'm debating whether I should change the reading code to cstdio as well, but it might get a bit more complicated, since I don't know the actual size of each line and the C input functions all require a char array buffer. Our saligner also does a bit more than readaligner. For example, our FASTQ file processing is much more robust and follows the FASTQ spec better.

Chris also mentioned that something else I can do is to replace the rank structure. He mentioned Daniel was looking at a 64-bit implementation, so that might be something to look at? I'll need to get more details later.

For tomorrow:

• A few more minor optimizations.
• Take a look at bowtie source code, maybe?

06/07/10

Run on skorpios with saligner commit 5f7f77021a321eab0019825bec36969209e707b6 and using the time command.

06/04/10

FM indexes

With my new-found knowledge, I took another look at the exact matching implementation in readaligner and now I feel like I get it a bit more. I tried to look up some papers on wavelet trees (I think that's what wv_tree stands for) and the rank functions, since that's pretty much the bottleneck but I couldn't really find much stuff on it (not that I tried that hard).

So for later:

• Find some literature on wavelet trees and take a look at the rank function some more
• Do some more reading...

06/03/10

Up until now, I've really only been testing exact matches with saligner over the 125,000 reads file. I was doing one final test on the 2M read file (to verify output and to do a preliminary speed test), I noticed a lot of disk activity after about 40k reads! I'm guessing this is what cache misses look like (it could just be it reading from the reads file, I'm not really sure on this, have to confirm), so I'm going to have to do some more profiling with Callgrind. Granted, my system doesn't have as much memory as skorpios, but this might be one of the areas we can improve on.

Other than that, I think it might be a good idea to start reading up on FM indexes, because I really can't do much else without knowing how the aligner works. So tomorrow, I might start on my reading (probably following in Daniel's footsteps).

To do:

• Check if we're really getting lots of cache misses.
• Start on my reading!

64-bit
Built on sycamore.umiacs.umd.edu
Sat Apr 24 15:56:29 EDT 2010
Compiler: gcc version 4.1.2 20080704 (Red Hat 4.1.2-46)
Options: -O3 -m64  -Wl,--hash-style=both
Sizeof {int, long, long long, void*, size_t, off_t}: {4, 8, 8, 8, 8, 8}

• In Exact and MismatchMappers, just make a cigar on the fly instead of passing stuff to SamEntry.

06/02/10

Another optimization I did today was make a convenience method in SamEntry to generate an entry representing an exact match (as per Chris' suggestion yesterday). Because exact matches are so predictable, it saves us having to generate the CIGAR, and a bunch of the tags. We also don't need to extract the original reference from the index, since it's not really needed. Mismatches should also be similar, since mismatches aren't represented differently from matches in CIGARs. However, we do have an MD:Z tag that supports mismatches (like BWA), but I'm not sure whether this is needed at all.

Finally, I added support for arguments in saligner, where we can specify reference, reads and output for convenience, so we don't have to recompile every time we need to use a different read file.

Anyways, I decided to run some tests on my newest changes using the 2622382 reads file just for completeness and for more accuracy, since I didn't do any tests on that file yesterday, either.

Under commit c6f55e33aa732c8952d1f56fa4c1fe6aa3875677, with the 2822382 reads file and Release build:

And with Bowtie (same version info as above):

To do:

• Do some profiling again...I haven't done that in a while and there's been a lot of changes to the code.
• Start reading up on FM index.

To do:

• In Exact and MismatchMappers, just make a cigar on the fly instead of passing stuff to SamEntry.
• More low level optimizations.
• Make a better saligner wrapper.
• Fix a bug in cigar generation (ticket 53df76)

06/01/10

Using commit 5a82392eeeb5653d1b75b7a8dd7ccdd499605cfa and saligner, with the 125000 readset (real time)

Using BWA with flags -l 99 -k 2 outputting to SAM:

Using readaligner with flags -k 2 --sam --fastq:

Bowtie with version info:

64-bit Built on sycamore.umiacs.umd.edu Sat Apr 24 15:56:29 EDT 2010 Compiler: gcc version 4.1.2 20080704 (Red Hat 4.1.2-46) Options: -O3 -m64  -Wl,--hash-style=both Sizeof {int, long, long long, void*, size_t, off_t}: {4, 8, 8, 8, 8, 8}

Note: I'm not really sure BWA is actually working right, because the output file looks way too small. More on that later.

Finally, I revisited the LocalAligner and tried to make a more efficient FullAlignment function by using arrays instead of vectors and speeding up the recursive Traceback function (it had a bunch of string-passing by value). I ended up failing, getting segmentation faults everywhere, but I think it was just because of me getting my columns and rows mixed up in my matrices. Anyways, more of that tomorrow, when I have a fresher mind.

To do:

• I'm going to table the SequentialMapper thing for now, since it's not very high priority. Chris also thought it might be a good idea to look into why MismatchMapper and GapMapper are so slow at k=0, instead of just replacing them with ExactMapper (ideally, ExactMapper shouldn't even be needed).
• Unabstracting Query - it's not very necessary on second look, since I usually use the specific Query I need instead of the general Query object anyway.
• Find more areas where reserving space will benefit (still not touched)
• Improve the AlignerInterface::FullAlignment functions.
• Benchmark BWA on skorpios.

05/31/10

Tests were done using the provided E. coli genome (NC_010473.fasta) on 375,000 reads (e_coli_125000_1.fq copied 3 times). Tested using the single-end aligner (saligner) release build, times taken from saligner output.

So there were slight speed improvements to exact and mismatch mappers and a negligible speed different in the gap mapper. The speed improvements are probably from some minor fixes I made during the cleanup, but I haven't really touched the readaligner Query classes, which takes up much more time in GapMapper than the other two mappers. Also, the more calls to the readaligner align functions in Exact/Mismatch mappers, the less speed improvement we see.

05/28/10

Preliminary testing on my own machine shows that saligner, using the GapMapper, is slightly faster than the old version (about 5-10% faster). Further tests need to be done on Skorpios to confirm. Another interesting point is that the GapMapper actually uses the LocalAligner to output to SAM format. However, it uses LocalAligner::FullAlignment functions to do so, which have not been optimized at all. I can also see GapMapper benefiting from using the BandedLocalAligner instead, since gaps tend to be small (BandedLocalAligner supports up to 15 gaps). It might also be worth vectorizing the FullAlignment function; currently, only ScoreOnlyAlignment is vectorized.

I also came across an interesting result, today. I noticed that at k=0, GapMapper, MismatchMapper and ExactMapper all have the same output (no difference whatsoever using diff), even though ExactMapper is faster than the other two. In fact, ExactMapper is something like 16x faster than GapMapper at k=0 on my machine! So, I think an interesting optimization would be to have ExactMapper do an initial alignment and using the more complex mappers only after it has returned negative. This follows through with existing behaviour anyway (both GapMapper and MismatchMapper will try the lowest k value, k=0, first before moving on to progressively higher k values). With data that is mostly exact, this could lead to a serious improvement, since most reads won't even need to touch the more complex alignment functions. This does lead to a question of whether or not the more complex mappers should even start at k=0 at all, since in most cases, ExactMapper will be there to cover it. The question is, I'm not sure if I should have Gap/MismatchMapper call ExactMapper directly within their Map functions, or whether I should let the user decide to call ExactMapper first. I'm currently leaning towards the latter option, since, as implied by their names, Mismatch/GapMapper are there to find alignments with gaps or mismatches, not to find alignments with exact matches. So, I'm proposing to start both mappers at k=1 instead of k=0, and let ExactMapper handle the k=0 condition (but the user has to do that explicitly). Finally, one more point to note is that even if those mappers start at k=1, they still find exact results, which leads me to think that Query::align checks for k=0 and k=1 when I call it for k=1 (and etc. for higher k). This might be an area of optimization worth looking into, since there would be a lot of overlap as we go higher in k value. However, it might not actually be possible to optimize this; I'm not familiar enough with the Query::align functions to know.

Finally, Daniel and I finished our discussion on class hierarchy and I think we're both pretty happy with the design. I've sent an email off to Chris for the finalized version (available on the discussion page of the Twiki).

To do:

• Add in a new mapper, SequentialMapper, which allows the user to pass in two mappers. SequentialMapper will then map using the first mapper, and use the second mapper only if the first mapper fails. This should introduce more flexibility in the aligner.
• Look into unabstracting the Query class and seeing how it fares (still haven't touched this yet).
• Benchmark on Skorpios
• Find more areas where reserving space will benefit (haven't touched this either)
• Make a vectorized LocalAligner::FullAlignment function.
• Make a BandedLocalAligner::FullAlignment function (both vectorized and non-vectorized).

05/27/10

In ExactMapper, I didn't really get a speed improvement. This was pretty much expected because the alignment function it uses is part of TextCollection (as opposed to the other two mappers, which use the Query class). The method, TextCollection::locate, pretty much searches the entire index for all occurrences of the sequence, not for the first x occurences. So, adding a max results parameter pretty much just prevents overflow matches from being outputted, which wasn't that slow to begin with.

MismatchMapper gave me a better result, though the news is bittersweet. With multimaps turned off (setting max_results to 1), I get about a 10-20% speed improvement over the old mapper (all numbers are preliminary, I haven't actually tested on Skorpios yet). However, with any sort of multimapping, even with max_results = 2, I get pretty significant slowdown WRT the old mapper (set to the same max_results value). This isn't really a regression, though; the old MismatchMapper code actually had different behaviour for multimaps. In cases where alignments are found at a lower mismatch value, the old mapper would not continue searching at higher mismatch values, even if there are fewer alignments than max_results. The new MismatchMapper actually continues searching until it's reached its maximum mismatch value or it's found the maximum number of results. This behaviour is more in line with readaligner, although I'm considering pushing it back to the old behaviour because it makes a little more sense to me (and is much faster at high max_results). I should discuss with Daniel/Chris about this.

While testing the new MismatchMapper, I noticed something interesting: if I increase the max_results value (it's report in readaligner) for both saligner and readaligner, I get about the same 2-4 second difference in runtime, no matter what value I use. Since increasing max_results effectively increases the number of alignment steps, the fact that the time difference is pretty much constant supports my previous conjecture that it's the pre/post-processing (parsing of FASTQ/FASTA files, outputting to SAM format, etc) of our aligner that's leading to the time difference. So, again, that should be an area worth investigating.

To do:

• Discuss class design with Daniel and come up with a finalized design we're both happy with. This may or may not mean changing the current design (which is good, except for the bit with Index and Mapper (haven't really done it today...)
• Look into unabstracting the Query class and see how that fares. If it's a big improvement, look into policy classes (also hasn't been visited at all)
• Add a max_results option to the GapMapper and move NGS_FMIndex::LocateWithGaps into GapMapper::Map.
• Find more areas where reserving space for strings and vectors would be useful.
• Benchmark the new implementations vs the old on Skorpios.
• Find other areas to optimize.

05/26/10

Daniel and I have also decided to put an option in the Mapper classes for the user to choose the maximum number of results they want. This is a minor optimization that would allow us to terminate the alignment function once a set amount of results have been reached, rather than going through with the whole alignment. This might also allow us to use fixed-length arrays instead of vectors (or at least reserve the space in the vector beforehand to prevent reallocs).

Another area I've investigated today was the cost of abstraction in our Index classes. I tried unabstracting the NGS_FMIndex and comparing runtimes. On 125k reads, I get about 9.8-10.1 seconds for both implementations, so the abstraction cost isn't really an issue. I just realized that Query objects are also abstracted those methods are called much more often. This might also be a good area to investigate.

Finally, I just realized that string 's can reserve space, too. I think this might be a good area to improve on, as I'm seeing a lot of calls to malloc from both string and vector. Considering both of these could get quite large in some cases, reserving memory beforehand may be a good idea so we don't get frequent realloc 's.

To do:

• Discuss class design with Daniel and come up with a finalized design we're both happy with. This may or may not mean changing the current design (which is good, except for the bit with Index and Mapper.
• Look into unabstracting the Query class and see how that fares. If it's a big improvement, look into policy classes.
• Add an option to limit number of results.
• Move instantiation of Query objects into the Mapper 's (if Daniel agrees)

05/25/10

In readaligner, they instantiate a single Pattern object for each read. In our saligner, we instantiate one Pattern per read per max mismatch value. When the Pattern is checked, it's first checked "as-is", and then reversed and checked again. Since readaligner's Pattern object is kept between k values, it will check the forward and reverse strands in different order than saligner. If a match is found, the alignment function terminates.

Unfortunately, this "discrepancy" isn't really a slowdown, since it's dependent on the given read set. For example, if all reads given match on the forward strand, then our saligner would be "faster". But, since I've finally found the error (or lack thereof), we can finally start to optimize.

To do:

• Optimize!
• Read up on policy classes (Chris mentioned this briefly but said it was a low priority, so I just wanted to jot it down before I forgot about it).

05/21/10

• On 125,000 reads, readaligner runs in about 9-10 seconds, while baligner runs in between 15-16 seconds.
• If I comment out the lines that actually do the alignment (the call to Query::align()), readaligner takes <1 second, while baligner takes about 3 seconds. This means baligner has a lot more overhead when starting and finishing (i.e. not actually "aligning"), but this still doesn't account for the 6-7 second difference.
• On 1000 reads, baligner and readaligner call Query::align() the same amount of times (if I use the Query class's own reverse functionality and comment out Daniel's). However, it seems calls after that start to diverge, but only slightly. For example, the call after Query::align(), MismatchQuery::firstStep(), gets called 3634 times on baligner, but only 3630 times on readaligner. It looks like this small difference is propagated downward through the calls, and ends in a difference of about 1000 more calls in baligner in a low-level function. This small difference in the beginning might be the source of the problem, since the increase in function calls could potentially get huge in 125,000 reads. However, I haven't profiled that yet, since it would take so long. Perhaps this is a good next step?
• baligner is much more complicated than readaligner in terms of class structure, so this might be giving us some overhead. Also, I've noticed vectors are used a lot when passing information around; this might be a source of overhead, especially since just the vectors seem to be returned, instead of a pointer to a vector. This could be a future point of optimization, but I doubt it would save on very much time.

I think the most promising area to look into right now is still the Query::align() area, and why it's calling MismatchQuery::firstStep() slightly more than the readaligner. I think I'm going to run a profile comparing 125,000 reads on both aligners, just to confirm if it is that portion of the code that's giving us the problem. Also, it might just be a bunch of small optimization problems and overhead from using vectors, etc. that could be resulting in the performance loss, although I kind of doubt it.

Todo:

• More profiling!

Afternoon update:

I realized an alignment size of 500nt isn't that realistic, since the vectorized banded version can only handle 255nt before it starts suffering from integer rollovers.

Anyways, I've also been profiling the baligner some more and trying to figure out why it's roughly 2x slower than the readaligner. I thought I had a breakthrough when I noticed baligner called Query::align() twice as many times as readaligner (for aligning reverse sequences), when Query already dealt with reverse alignments, but it turned out Daniel had already accounted for that. So, I'm pretty much back to square one!

To do:

• More profiling to figure out what's wrong; maybe play with report, k, and iteration limit values.
• Try making CMake compile different implementations of the local aligner based on CPU features.

05/20/10

So in the end, there's quite a big difference between the regular and vectorized versions. It's interesting to see the improvement in the normal vectorized aligner is greater than the banded aligner, even though the banded aligner does twice as many calculations per vector. It's probably due to me having to code in my own max function in the banded aligner.

05/19/10

Note these were not done on a Release build because the Release build messes up the vectorized banded aligner for some reason (outputs all 0). I'll have to look into why building as "Release" build is wrong but a regular build is not. I think it might be going wrong in my "custom" signed max function.

I also started profiling the baligner, after making sure all the output was correct. Since aligning the 125000 reads Daniel usually uses would have been way too slow, I used the e_coli_1000_1.fq reads file as a test. I ran memcheck first and fixed all the memory leaks that arose. I also ran callgrind over baligner and I found a significant bottleneck in the gap aligner, where a constant 256*256 matrix was being initialized and reinitialized multiple times over all the reads, when it only needed to be initialized once. After the fix, the speed increase was significant:

For tomorrow:

• Find out why the vectorized banded local aligner isn't working on a Release build and fix it.
• Benchmark the local aligner using Release builds.
• Do more profiling of the baligner.

05/18/10

05/17/10

I also cleaned up my implementations and optimized the existing local aligner implementation (only the ScoreOnlyAlignment method) to be more memory-efficient and faster. All public interfaces for both vectorized/non-vectorized versions of the aligners now match in both input and output and should be more-or-less equivalent in usage (with the exception of the free gap bits in the regular local aligner).

I'm not really sure what else to do for now, so I'll probably have to ask Chris/Daniel tomorrow.

05/14/10

For Monday:

• Edit the CMakeList.txt files so that it compiles vector/non-vector versions based whether SSE2 is available or not.
• Polish up all the code I've written and merge.
• Add in my improvements to the original local aligner.

05/13/10

I came in today and realized my method of finding the bounding lines on the banded non-vectorized version was pretty complicated; instead of finding the equations for the two bounding lines given some k distance of separation from the diagonal (which involves finding the slope and y-intercept and has square roots in it), I instead take the diagonal line and calculate all cells that are a horizontal distance k from it. This is much easier, and works out to pretty much the same thing. It has the added advantage of letting me take in a horizontal k value for my band radius, which will let me guarantee a correct score given the number of gaps <= k.

I managed to finish up the non-vectorized version as well as make up a (I think) pretty comprehensive test suite for it. In the end, I passed all my own tests, though I did miss a "-1" in one line of code that had me banging my head for a while. I also managed to improve the memory efficiency of the score-only local alignment function in general based on a few tricks I learned while combing through the vectorized Smith-Waterman code, though I've only implemented them in my banded aligner for now. Improvements include:

• Keeping only one read-sized row of scores instead of 2.
• Keeping only one read-sized row to keep track of gaps in the x-axis instead of 2 rows.
• Having just one int to keep track of the last match, instead of 2 read-sized rows.
• Having just one int to keep track of the last gap in y, instead of 2 read-sized rows.

To do:

• Start tackling the vectorized banded algorithm! That'll be pretty challenging and probably keep me busy all day.
• (Much later) Add in my improvements to the original local aligner.
• (Still much later) Polish up the different aligners and have CMake compile either the vectorized or non-vectorized versions of the aligners based on CPU support.

I do remember doing a "flip" in CS320, but I think it was with a square, whereas this matrix could be rectangular (shouldn't be that much different). I should look up my old notes later on, too.

For tomorrow, I will:

• Finish up my banded non-vectorized algorithm and test it
• Try coming up with a vectorized version (may involve having to "flip" the matrix).

05/12/10

I did a few preliminary tests, using a randomly generated genome/read sequences of length 1000 each. Over 1000 iterations, my results were:

• 38.836576365 seconds using affine gap penalty
• 28.907083597 seconds using linear gap penalty

Now it's on to doing the "banded" algorithm. I'm going to try making a non-vectorized version first, then move on to vectorized. I'm not really sure how I'm going to test it though, since the results are going to be a little bit different...

• The SHRiMP code required a call to initialization, which was troublesome and not really required given the way the NGSA LocalAligner works.

05/11/10

• I'll have to use that code again, but modified when I do the full alignment (with traceback).
• The SHRiMP code required a call to initialization, which was troublesome and not really required given the way the NGSA LocalAligner works.
• Since SHRiMP was in C, I couldn't wrap their functions in namespaces.

The good news is, I got all the tests to pass on the score-only alignment. Now all I have to do is get a full alignment working. I also want to try to optimize for the linear gap penalty method, since that case is much simpler than the affine gap penalty. Although the linear case is just affine with a gap open penalty set to 0, I think I can do away with some of the vectors in the calculation, and so make it a little faster. That'll be a low priority, though. So for tomorrow, I think I'll be able to finish:

• The full Smith-Waterman alignment (vectorized)
• Hopefully get that to pass all the tests (i.e. make the behaviour the same as the non-vectorized version)
• Take a look at optimizing for linear gap penalty if I have time.

05/10/10

I'm about a quarter complete on a score-only vectorized version for the local aligner, and I've got a basic game plan laid out for completing it; I've also found that the best way to understand some code is to try to code it yourself! I now have a much firmer grasp on the whole SHRiMP aligner and on why it does some of the things it does (for example, one question I had was why it stored values for b_gap, but no a_gap array). I think by tomorrow, I can:

• Get the score-only vectorized Smith-Waterman finished and hopefully tested out.
• Time permitting, get started on the full alignment with traceback.
• Read up on CMake a bit more; it's still giving me a strange amount of trouble...

05/07/10

Chris also gave me the source code for a vectorized Smith-Waterman function (from the SHRiMP library) and I spent the rest of the afternoon looking at that and trying to understand it. I found the documentation for the SSE2 functions on the Intel website (link to come later, not on the right now). From my understanding (and Chris' explanation), the function goes through the matrix using a diagonal vector ("positive slope") operating on 8*16bit values, and calculates the highest score from that vector. It uses the affine gap penalty method to calculate scores. From my discussion with Chris, I should now be working on:

• Integrating the vectorized Smith-Waterman function into the NGSA library (I can either try using the source code given or make my own).
• Configuring CMake so that it detects support for SSE2 (and maybe SSE3) so that it will compile either the vectorized Smith-Waterman or non-vectorized Smith-Waterman, depending on the machine.
• Take a look at the banded Smith-Waterman function and try integrating that, too, perhaps. This function only checks a "strip" of cells going down the diagonal of the score matrix, since in practice, alignments aren't as useful when there are large variations from a full match.
• Finally, in the distant future, start thinking about the free gap algorithms and vectorizing them.

End-of-the-day post:

I've managed to modify Daniel's test/correctness_test.cc so that it works on the existing code, and all the test cases pass. The test cases Daniel wrote were pretty complete but only tested the ScoreOnlyAlignment method in the local aligner. So next up will be writing some test cases for the FullAlignment, which may be a little harder since it doesn't just output numbers. I've written a few tests so far, and I think I may have found one or two bugs (but I'm going to have to look into that some more to see if it was my test case or the code that was failing). So tomorrow, I will mostly be doing...

• Finishing up the test cases for the FullAlignment method.
• Fixing any bugs that might come up along the way.
• If I get a chance...maybe even start looking into the vectorized methods!

05/06/10

• Read through a bunch of Daniel's journal entries to get myself a little more familiar with the project so far.
• Had a meeting with Daniel and Chris where we discussed my first task:
  • First, I should make or beef up the existing test suite.
  • Then, I'm to make a vectorized version of the Smith-Waterman aligner that Daniel has working.
  • I should also familiarize myself with existing tools (e.g. Bowtie, BWA, Mosaik, and SamTools).
• I ended up spending a lot of time reading through the various guides Daniel posted, as well as the Git guide and the Google C++ style guide to get myself more familiarized with what I should be doing and how I should be doing it.
• I also spent some time reading up on the Smith-Waterman algorithm and reading up on sequence alignment in general to get myself a little more familiar with the project.
• Checked out Daniel's existing code, especially local_aligner.cc bit (which I'm assuming my vectorized SW will be based a lot off of). I can understand the majority of the code, although I'm a little iffy on some of the gap penalty calculations.

So after all that, I think I'm finally ready to do a little coding. Today, I'll mostly be working on the test suite and trying to get that out of the way.