2
$\begingroup$

I'm writing the code from scratch and not using libraries for the actual indexing/search since it's a project for school. Any blackbox advice is fairly useless.

I'm planning on aligning some short reads with the D. melanogaster genome but, since I was thrown into this project with very little info/ experience I'm stuck on a few things.

Firstly, for preprocessing the fasta file:

  • What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?

  • As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?

Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?

On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?

Is there any way to validate the alignment quickly at the end? How much error should I be allowing?

This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.

Any help would be hugely appreciated!

$\endgroup$
2
  • $\begingroup$ [STAR][1] uses suffix arrays to map RNA-seq reads. It is both very fast and quite accurate. Perhaps their manuscript may help with your questions. [1]: github.com/alexdobin/STAR $\endgroup$
    – GWW
    Apr 20 at 12:35
  • $\begingroup$ AH, sorry I did not notice the "it's a homework" part. However, note that practically the mapping tools out there are open source, so you can actually browse the code and check what exactly they did. $\endgroup$
    – Kamil S Jaron
    Apr 21 at 11:20
2
$\begingroup$
  • What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?
  • As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

Good questions. I think this depends on the scope of your class project. Frankly, I think it would be best to just start with unambiguous sequences (only ATGC) and not worry about potential repeats. It's also a good software engineering exercise to solve the easy case first in a way that's flexible and allows you to expand to the additional complexities above. For example, the original Bowtie software considered references with ambiguous characters invalid. You can see for yourself here how Bowtie2 handles it.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?

Yes, the most naive implementation of a search with a suffix array will be $O(m\log n)$ time, where $m$ is the read length and $n$ is the genome length. However, with enhanced suffix arrays you can also obtain the same $O(n)$ search complexity of suffix trees. There is also the use of suffix arrays to construct Burrows-Wheeler transforms and even further compression with the FM index. There are open source implementations of all of these data structures that may be interesting to you.

Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?

Suffix arrays store the indices of the start of each suffix in order to save space, so each entry is a single integer.

On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?

Its worth trying to compute this yourself (: You can determine the number of bytes required per base pair in your input and see how that scales with your available memory. For example, if your suffix tree implementation required $20n$ bytes, you'd end up w/ roughly 2Gb needed for a fruit fly genome.

Is there any way to validate the alignment quickly at the end? How much error should I be allowing?

Again, this depends on the scope of the project. Allowing for errors in mapping reads is an active area of research. See the supplemental of the Bowtie2 paper for a dynamic programming algorithm that handles a certain number of mismatches.

This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.

All good questions, and I think many answers depend on the scope of the project. Reimplementing a read mapper is no small feat, especially when you start considering all of the different complexities such as sequencing error and ambiguous bases. I would strongly consider starting with looking at just the pure string problem of read mapping: Only consider ATCG in both reference and reads, and don't even consider reads who do not have exact matches to the reference. From there, time permitting, you can start adding on to your solution. As far as reading goes, there's a deep rabbit whole of these topics, from tools like bwa-mem, bowtie2, mini-map, to the more theory-oriented papers discussing optimized versions of the data structures. I can attach some here if you would like. Hope these answers helped!

$\endgroup$
2
  • $\begingroup$ Thank you! This was really helpful $\endgroup$
    – nibs
    Apr 22 at 7:55
  • 1
    $\begingroup$ @nibs feel free to accept the answer if it sufficiently answers your questions (: $\endgroup$ Apr 22 at 19:27
3
$\begingroup$

Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.

I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an existing solution.

By the way, the Ns in the reference are either hard-masked nucleotides (that would be repeats of some sorts) or (more often) simply unknown nucleotides. So there is no way to map your reads there anyway (and the mentioned mappers will handle the Ns for you).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.