http://qs321.pair.com?node_id=1120876


in reply to Re: Is it possible to make reference to substrings and sort them?
in thread Is it possible to make reference to substrings and sort them?

It's for bioinformatics, so the alphabet is 4 letters, and the text can be up to 5M characters in my case.

Replies are listed 'Best First'.
Re^3: Is it possible to make reference to substrings and sort them?
by BrowserUk (Patriarch) on Mar 22, 2015 at 13:12 UTC
    the text can be up to 5M characters

    No matter how you generate them, just storing the 5 million rotations of 5 million characters (>23 Terabytes of data) is going to be a problem, long before you get to the point of trying to sort.

    If the purpose of your BWT generation is to aid compression, then there are much more effective and efficient methods of compressing genomic data.

    If the purpose is Short Read Alignment, you almost certainly need to look at the suffix trie methods described there.


    With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority". I'm with torvalds on this
    In the absence of evidence, opinion is indistinguishable from prejudice. Agile (and TDD) debunked

      Yes, this is why I wouldn't want to generate all the rotations.

      Short read alignments is what I need to do. I have a few thousand short sequences, length of ~20nt. I want to find those that are present in my text, as well as their position. BWT is the way to go.

      Suffix trie is not efficient in terms of memory, as the trie needs to store all suffixes from 1 to length(text). Thank you for the link to the paper.

        Suffix trie is not efficient in terms of memory, as the trie needs to store all suffixes from 1 to length(text)

        Hm. Maybe I'm misreading the paper, but my reading suggests that constructing a BWT the classical method -- construct all the rotations and then sort them -- takes quadratic time, which is why the normal practice is now to construct a prefix trie first; and the BWT from that (to save time).

        And further, the problem of suffux tries being inefficient of space ("nlog2n bits of working space which amounts to 12 GB for human genome") is fixed by Hon et al. (2007) which requires "only requires <1 GB memory at peak time for constructing the BWT of human genome".

        The algorithm shown in Figure 2 is quadratic in time and space. However, this is not necessary. In practice, we usually construct the suffix array first and then generate BWT. Most algorithms for constructing suffix array require at least nlog2n bits of working space, which amounts to 12 GB for human genome. Recently, Hon et al. (2007) gave a new algorithm that uses n bits of working space and only requires <1 GB memory at peak time for constructing the BWT of human genome . This algorithm is implemented in BWT-SW (Lam et al., 2008). We adapted its source code to make it work with BWA.

        Anyway. How many mismatches do you allow and still consider a site to match one of your 20nt short reads?

        I have code that can find all the 1-mismatch sites for each of 10,000 x 20nt short reads. in a 5Mnt strand, in 1.5 seconds using 1/4GB memory:

        C:\test>1120869 -BIG=5e6 -FUZZ=1 -MIN=10 -NSHORTS=10000 | wc -l Indexed search took 1.509765148 secs 2

        Allowing two mismatches per requires 6 1/2 minutes :

        C:\test>1120869 -BIG=5e6 -FUZZ=2 -MIN=6 -NSHORTS=10000 | wc -l Indexed search took 395.132736921 secs 80 C:\test>1120869 -BIG=5e6 -FUZZ=2 -MIN=6 -NSHORTS=10000 AAACTTGTTACCGGGGTACT ~ TAACTTGTTACCGGGTTACT ( AACTTGTTACCGGG TACT) [18 + of 20] at 1423486 CTCGCCATGTAAGTTATGGT ~ CCCGCCATGAAAGTTATGGT (C CGCCATG AAGTTATGGT) [18 + of 20] at 3739837 ACTCAGGAGCATATGAGCAA ~ ACTCGGGAGCATATGAGCAG (ACTC GGAGCATATGAGCA ) [18 + of 20] at 117017 GCGCACTACAACCTCTTAAA ~ GCGCACTACAACGTCTTATA (GCGCACTACAAC TCTTA A) [18 + of 20] at 2758181 CCAGTATATGATACGAACGC ~ CCAGTATCTGGTACGAACGC (CCAGTAT TG TACGAACGC) [18 + of 20] at 2681307 CTCGCACGATTACGACCGAA ~ CTCGCAGAATTACGACCGAA (CTCGCA ATTACGACCGAA) [18 + of 20] at 4658451 ...

        And allowing 3 requires 30minutes:

        C:\test>1120869 -BIG=5e6 -FUZZ=3 -MIN=5 -NSHORTS=10000 | wc -l Indexed search took 1763.693833113 secs 1455

        With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority". I'm with torvalds on this
        In the absence of evidence, opinion is indistinguishable from prejudice. Agile (and TDD) debunked