Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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):

enter image description here enter image description here

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):

enter image description here enter image description here

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).

enter image description here

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:

enter image description here

The code proposed by george2079 achieves:

enter image description here

@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}]
share|improve this question
1  
Can you provide complete code for a test case? –  Michael E2 Mar 17 at 13:04
    
I'm sorry for the many edits. had to fix a logical thing :) Thanks for help. –  g3kk0 Mar 17 at 13:33
    
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? –  Leonid Shifrin Mar 17 at 17:44

5 Answers 5

up vote 9 down vote accepted

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 *)
share|improve this answer
    
thanks very much! That piece of code looks perfect for my needs. I will have to test it on my data. –  g3kk0 Mar 17 at 13:53
3  
My goodness, what a difference it makes to work with lists rather than strings! This is something to keep in mind. –  David Carraher Mar 17 at 13:54
    
@DavidCarraher Oh yeah :) –  Leonid Shifrin Mar 17 at 13:55
    
@g3kk0 Thanks for the accept, but you could've waited for some more to encourage others to post - perhaps, there are still better solutions. –  Leonid Shifrin Mar 17 at 13:56
    
@Leonid Shifrin: sure, why not. I have unchecked your answer and keep it open for a while. –  g3kk0 Mar 17 at 14:06

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 *)
share|improve this answer

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.

share|improve this answer
    
Actually, on my system (MacBook Pro Core I7, MaC OS X 10.7.5), this is about 1.5x slower than my suggestion. I wonder why such a platform dependence in this case. Anyway, this is really fast, and I like that you stay in the string domain, a definite +1. –  Leonid Shifrin Mar 18 at 13:51
    
@LeonidShifrin: Including the time to get positions that yours needs/uses? Interesting - might be an amusing study re: timings vs platforms. Thx. for +1 as always. –  rasher Mar 18 at 21:11
    
My bad. When I include the that, and StringSplit in your case, then the timings are about the same (mine a tiny bit faster on my machine). –  Leonid Shifrin Mar 18 at 23:07
    
@LeonidShifrin: Thanks for checking - I'm hoping you follow-up on your lazy generation comment, always educational to see your ideas. –  rasher Mar 18 at 23:10

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..

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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