What is the fastest way to extract a substring from a huge string?

Question

I develop an algorithm to find the most common substrings of a given length n up to a certain distance l in front of a pattern p in a huge sequence S. E.g. I want to find the most common hexamers (substring of length 6) that occur in front of the string pattern p="TTAA" at a distance of at most l=1000 characters in front of the pattern "TTAA" in S. The steps of the algorithm are:

1. Import the sequence S.
2. Detect the first positions i∈{1,...,k} of all occurrences of pattern p in S.
3. Extract the sites (substrings) si in front of each pattern position i from S.
4. Determine all nmers (substrings of length n) for each site si.
5. Select the nmer from all those substrings that is the most common and determine how many sites can be covered with that nmer.
6. Delete the lists that are covered.
7. Repeat 5. and 6. until all lists are covered or stop after a given number of iterations.

The result should be a list of the nmers that are the most common ones in decending order of coverage.

Consider this string for the variable sequence and the following code. I first determine the positions in S where the pattern p occurs:

sequence = StringTake[Import["http://ge.tt/api/1/files/9teGemr/1/blob?download"], 200000];
    pattern = "TTAA";
    positions = StringPosition[sequence, pattern, IgnoreCase -> True];

Now I want to extract all the substrings si of length 1000 characters in front of the pattern p that appears at position i in sequence. In principle, I already have a solution for the problem but it appears to be quite slow if S is large and the pattern matches many positions (many si have to be extracted):

substrings = Table[
    If[
     First[i] < 1001,
     StringTake[sequence, {1, First[i]-1}],
     StringTake[sequence, {First[i] - 1001, First[i] - 1}]
     ], {i, positions}
];

For the runtime I have added some tests here without (left) and with (right) parallelization (ParallelTable instead of Table):

However, as soon as the number of found positions for the pattern and the length of the sequence increases the runtime of this approach explodes. For example if I run the code for the whole sequence of 1 million characters without (left) and with (right) parallelization (ParallelTable instead of Table):

Is there a faster way to do this?

Edit 17/03/2014

The code proposed by Öskå is a little bit faster, but still I would need this code to run for much longer sequences (200 million characters) and many more hits (about 3 million).

Edit 18/03/2014

The code proposed by Leonid Shifrin seems to be a hell of a fast piece of code :) I tested it for a sequence of 20 million characters and got the following results:

The code proposed by george2079 achieves:

@LeonidShifrin comment:

By the way, for 200 million chars and a few million hits, you'll really need some lazy infrastructure, because just holding all your extracted substrings in RAM at the same time won't be possible for such large objects. Are you interested in this kind of solution?

Yes indeed, I already thought about this. Later on, each extracted si has to be splitted into overlapping substrings of length 6. I currently do this using Partition on the list of characters. My idea was to dump save each of those strings si into an MX file and then move on file by file:

ParallelTable[
 (
  cur = FromCharacterCode /@ 
    Partition[ToCharacterCode[sites[[i]]], nSize, 1];
  Export[FileNameJoin[{path, 
     "hexamers-site-" <> 
      StringJoin[ToString /@ (PadLeft[IntegerDigits[i], 9, 0])] <> 
      ".mx"}], cur]
  ), {i, 1, Length@sites}]

Answer 1

Your main inefficiency seems to be in actually extracting the found strings, not in finding the positions. Here is a version based on ToCharacterCode - FromCharacterCode, which is quite a bit faster:

ClearAll[extractAll];
extractAll[s_String,startpos_List,size_Integer?Positive]:=
    With[{codes=ToCharacterCode[s]},
        With[{len=Length[codes]},
            (FromCharacterCode[Take[codes,{Clip[#1-size-1,{1,len}],#1-1}]]&)/@
                startpos
        ]
    ]

So that

(sub = extractAll[sequence, positions[[All, 1]], 1000]); // AbsoluteTiming

(* {0.147719, Null} *)

sub === substrings

(* True *)

Answer 2

Not heavily tested, but ~1.6X faster than fastest so far (including the time in that to get positions) in my quick tests, keeps things in the string domain:

ss = StringSplit[sequence, pattern];

result = Most@ FoldList[(s = StringJoin[#, pattern, #2]; 
     If[StringLength@s >= 1000, StringTake[s, -1001], s]) &, ss[[1]], Rest@ss]

In any case, as Leonid Shifrin points out in comments, any of the answers will probably blow chunks memory-wise for an actual 200M length case. Probably something using file-backed data and streams might be applicable (just a sketch of the idea):

strm = OpenRead["c:\\users\\rasher\\downloads\\testseq.txt"];

pattern = "TTAA";

nextstr[stream_] := 
  Read[stream, Record, RecordSeparators -> {pattern}];

(* get a chunk of results - massage as desired for 1K length, etc. *)
chunk = Table[nextstr[strm], {3000}];

(* do your thing... rinse and repeat as desired *)

Close[strm];

Streams are quite quick, e.g., getting a chunk of 3K results takes a tenth of a second on my goof-off netbook.

Answer 3

Here is an other solution (even if it's much slower than Leonid's :))

pos = GatherBy[positions, First@# < 1000 &];
split = StringSplit[sequence, ""];
mysub = Flatten@{
  StringJoin@split[[1 ;; # - 1]] & /@ First /@ First@pos,
  StringJoin@split[[# - 1001 ;; # - 1]] & /@ First /@ Last@pos};)

sub === mysub
(* True *)

Answer 4

One approact is the use regular expressions:

 StringTake[ # , {1, -5}] & /@ (StringReplace[ # , "$" -> ""] & /@ 
         StringCases[StringJoin[Table["$", {1001}]] <> sequence, 
               RegularExpression[".{1001,1001}TTAA"], Overlaps -> True]);

 substrings == %

True

this takes 0.4 sec for me for the 200k example. For the full 1,000,000 character set I get 1.97 sec vs 82 using the code in the question.

Note I handled the short sequence at the beginning cases by prepending a long string of "$$$.." and then stripping it off the results. I was unable to cook up a regexp to encompass everything, something like this ought to work it seems..

RegularExpression["((.{1001,1001})|(^.{1,1001}))TTAA"]

It might get a tad faster if someone could figure out the correct incantation..

Answer 5

Build a table upon insert.

Otherwise, for any pattern << string, then find the longest distinct letter set in pattern. Then step in blocks in string for that letterset. If a letter of the set is found, then check the nearest neighbors.

You can modify the above algorithm to choose the letterset based on commonality of occurrences of a letter, E.g. If the letterset is "etions" and the string English input, then you'll have lots of potential matches where you must compare neighbors; but, if the letter set is "JXZQK" then you will be much, much faster. Refer to oxforddictionaries website for a list of most common and least common letters in English.

A third approach is the sort the incoming string into some sort of ordered tree (3-4-5) and then apply your search against that tree. Obviously, the sorting algorithm would still be based on how common the letters are in the expected sources and in the expected match patterns. Fun project.