Re: Identifying Overlapping Matches in Nucleotide Sequence

G'day FIJI42,

Biological data are typically huge. For reasons of efficiency, when dealing with this type of data, you should choose a fast solution over a slower one. Perl's string handling functions (index, substr, etc.) are measurably faster than regexes (try it yourself with Benchmark). String operators (e.g. $x eq 'some literal') will be far more efficient than an equivalent regex (e.g. $x =~ /^some literal$/).

There are times when you will need a regex. In these cases, don't give the regex more work than necessary. In "/(AA)/gi", you're asking it to pointlessly capture matches (you never use the captured data); and you're also asking for a case-insensitive match (again pointless, as all characters are uppercase). Where you are dealing with mixed-case, or unknown case, canonicalise your data once then work with that (see fc, uc, and lc).

With your specific question, index would be the best choice.

Your ultimate goal is somewhat unclear. If you want to match "AA" once in "AAA", and twice in "AAAA":

$ perl -E '
    my ($x, $y, $c, $p) = qw{AxAAxAAAxAAAA AA 0 0};
    $c++, $p += 2 while ($p = index $x, $y, $p) > -1;
    say $c;
'
4
[download]

If you want to match "AA" twice in "AAA", and thrice in "AAAA":

$ perl -E '
    my ($x, $y, $c, $p) = qw{AxAAxAAAxAAAA AA 0 0};
    $c++, $p++ while ($p = index $x, $y, $p) > -1;
    say $c;
'
6
[download]

Note that the only difference between those two is: in the first, the position ($p) is incremented by two ($p += 2); in the second, it is incremented by one ($p++).

That code is really only about quickly demonstrating the index function. It's not intended as a production-grade solution; for instance, pick better (more meaningful) variable names.

— Ken

Comment on Re: Identifying Overlapping Matches in Nucleotide Sequence Select or Download Code

Replies are listed 'Best First'.
Re^2: Identifying Overlapping Matches in Nucleotide Sequence by roboticus (Chancellor) on Oct 27, 2017 at 17:13 UTC
kcott: When you mentioned the overhead of capturing and case-insensitivity, it piqued my interest, so I put together a quickie, and it shows: $ time perl ../ex_pattern_matching.pl <GGCGATGGCGTTCTAGCGCGTAAAA> x_cap_i: 7396 x_cap: 7396 x_i: 7396 x: 7396 poz: 9202 idx: 9202 ovr_cap_i: 9202 ovr_cap: 9202 ovr_i: 9202 ovr +: 9202 Rate ovr_cap ovr poz ovr_cap_i ovr_i x_cap x_cap_i x_i + idx x ovr_cap 11.1/s -- -11% -59% -60% -62% -71% -73% -73% + -91% -92% ovr 12.4/s 12% -- -54% -55% -57% -68% -69% -70% + -90% -91% poz 26.9/s 143% 117% -- -2% -7% -31% -34% -35% + -78% -79% ovr_cap_i 27.5/s 148% 122% 2% -- -5% -29% -32% -34% + -78% -79% ovr_i 28.8/s 160% 133% 7% 5% -- -26% -29% -31% + -77% -78% x_cap 38.8/s 250% 213% 44% 41% 35% -- -4% -7% + -69% -70% x_cap_i 40.5/s 266% 227% 51% 48% 41% 4% -- -3% + -67% -69% x_i 41.6/s 275% 236% 55% 51% 44% 7% 3% -- + -66% -68% idx 123/s 1012% 895% 358% 348% 327% 217% 204% 196% + -- -6% x 131/s 1079% 955% 386% 375% 353% 237% 222% 214% + 6% -- real 0m37.721s user 0m37.656s sys 0m0.015s [download] The entries starting with "x" have the same problem that the OP had: the count was too low since they skip overlapping matches. The idx solution you provided was the clear winner (so long as you want correct results). I was surprised at the interaction between case insenstivity and capturing: turning on either makes it slow, but the combination isn't any slower. Since the simple (but incorrect) match was as quick as the index solution, I tried to "fix" the overlap by using the pos function on the simple match (the poz entry), but mucking about with the pos function put poz into the realm of the other overlapping match functions. The code: Read more... (3 kB) ...roboticus When your only tool is a hammer, all problems look like your thumb.	[reply] [d/l] [select]
Re^3: Identifying Overlapping Matches in Nucleotide Sequence by kcott (Archbishop) on Oct 27, 2017 at 21:12 UTC
G'day roboticus, [You should be safe to drink coffee while reading this one. :-) ] Thanks for writing that benchmark. I have previously posted benchmarks (as have others) comparing string functions with regexes: I believe the results are fairly well known. I left writing another one as an optional exercise for the OP. It was, however, interesting that the various combinations of capturing and case-sensitivity didn't make that much difference (as you noted). I got the same sort of results as you (using Perl v5.26.0): <ACACGAAGCGCTCGTGTGATTATCT> x_cap_i: 7521 x_cap: 7521 x_i: 7521 x: 7521 poz: 9496 idx: 9496 ovr_cap_i: 9496 ovr_cap: 9496 ovr_i: 9496 ovr +: 9496 Rate ovr_cap ovr ovr_cap_i poz ovr_i x_cap x_cap_i x_i + idx x ovr_cap 5.77/s -- -12% -69% -70% -71% -80% -81% -83% + -92% -92% ovr 6.60/s 14% -- -64% -66% -67% -77% -79% -80% + -91% -91% ovr_cap_i 18.4/s 219% 179% -- -4% -8% -36% -41% -45% + -75% -76% poz 19.1/s 231% 190% 4% -- -5% -34% -39% -43% + -74% -75% ovr_i 20.1/s 248% 204% 9% 5% -- -30% -36% -40% + -73% -74% x_cap 28.8/s 399% 337% 56% 51% 44% -- -7% -14% + -61% -62% x_cap_i 31.2/s 440% 372% 69% 63% 55% 8% -- -7% + -58% -59% x_i 33.3/s 477% 405% 81% 74% 66% 16% 7% -- + -55% -56% idx 74.1/s 1183% 1023% 302% 287% 269% 157% 138% 122% + -- -3% x 76.3/s 1222% 1057% 315% 299% 280% 165% 145% 129% + 3% -- [download] By the way, I thought this was rather clever: `$cnt++ while $pos = 1 + index $str, 'AA', $pos;` [download] Out of interest, I added these two subroutines: `sub lc_x_cap_i { my ($str, $cnt) = @_; ++$cnt while $str =~ /(aa)/gi; return $cnt; } sub lc_x_i { my ($str, $cnt) = @_; ++$cnt while $str =~ /aa/gi; return $cnt; }` [download] Then ran your benchmark code again with only the "`lc_`" and "`x`" routines. I ran it a few times; this seemed to be a fairly representative result: `<AGGTATGATGGTGTAGAGTAACTAG> x_cap_i: 7506 lc_x_cap_i: 7506 x_cap: 7506 x_i: 7506 lc_x_i: 7506 + x: 7506 Rate lc_x_cap_i x_cap_i x_cap lc_x_i x_i + x lc_x_cap_i 27.7/s -- -1% -12% -13% -14% + -66% x_cap_i 28.0/s 1% -- -11% -12% -13% + -66% x_cap 31.6/s 14% 13% -- -0% -2% + -61% lc_x_i 31.7/s 15% 13% 0% -- -2% + -61% x_i 32.4/s 17% 16% 2% 2% -- + -60% x 81.3/s 193% 190% 157% 156% 151% + --` [download] So clearly "x" (`/AA/g`) was faster than the rest (being ~80/s whereas the others were ~30/s, in all runs). I'd say the others were too close to call: although there may appear to be a trend, in two runs "x_cap" was the slowest. — Ken	[reply] [d/l] [select]


XP is just a number
	PerlMonks