2

I have a large string file seq.txt of letters, unwrapped, with over 200,000 characters. No spaces, numbers etc, just a-z.

I have a second file search.txt which has lines of 50 unique letters which will match once in seq.txt. There are 4000 patterns to match.

I want to be able to find each of the patterns (lines in file search.txt), and then get the 100 characters before and 100 characters after the pattern match.

I have a script which uses grep and works, but this runs very slowly, only does the first 100 characters, and is written out with echo. I am not knowledgeable enough in awk or perl to interpret scripts online that may be applicable, so I am hoping someone here is!

cat search.txt | while read p; do echo "grep -zoP '.{0,100}$p' seq.txt | sed G"; done > do.grep

Easier example with desired output:

>head seq.txt    
abcdefghijklmnopqrstuvwxyz

>head search.txt
fgh
pqr
uvw

>head desiredoutput.txt
cdefghijk
mnopqrstu
rstuvwxyz

Best outcome would be a tab separated file of the 100 characters before \t matched pattern \t 100 characters after. Thank you in advance!

9
  • 1
    What when there isn't a 100 chars before (or after) the pattern? (Like when a pattern is found at a position 50 in the string, or at 60 before the end.) Take whatever padding there is?
    – zdim
    Nov 9, 2020 at 8:46
  • If there are multiple identical matches, print one or all or count? Nov 9, 2020 at 8:51
  • @zdim, they are all guaranteed to be within 100 characters from the ends
    – geri
    Nov 9, 2020 at 8:58
  • @JamesBrown there are no identical matches in my case, but if there was then print all
    – geri
    Nov 9, 2020 at 8:58
  • And seq.txt is just one line? Nov 9, 2020 at 9:01

4 Answers 4

4

One way

use warnings;
use strict;
use feature 'say';

my $string;

# Read submitted files line by line (or STDIN if @ARGV is empty)
while (<>) {
    chomp;
    $string = $_;    
    last;          # just in case, as we need ONE line
}
# $string = q(abcdefghijklmnopqrstuvwxyz);   # test

my $padding = 3;  # for the given test sample

my @patterns = do { 
    my $search_file = 'search.txt';
    open my $fh, '<', $search_file or die "Can't open $search_file: $!";
    <$fh>;
};
chomp @patterns;
# my @patterns = qw(bcd fgh pqr uvw);  # test

foreach my $patt (@patterns) {
    if ( $string =~ m/(.{0,$padding}) ($patt) (.{0,$padding})/x ) {
        say "$1\t$2\t$3";
        # or
        # printf "%-3s\t%3s%3s\n", $1, $2, $3;
    }
}

Run as program.pl seq.txt, or pipe the content of seq.txt to it.

The pattern .{0,$padding} matches any character (.), up to $padding times (3 above), what I used in case the pattern $patt is found at a position closer to the beginning of the string than $padding (like the first one, bcd, that I added to the example provided in the question). The same goes for the padding after the $patt.

In your problem then replace $padding to 100. With the 100 wide "padding" before and after each pattern, when a pattern is found at a position closer to the beginning than the 100 then the desired \t alignment could break, if the position is lesser than 100 by more than the tab value (typically 8).

That's what the line with the formatted print (printf) is for, to ensure the width of each field regardless of the length of the string being printed. (It is commented out since we are told that no pattern ever gets into the first or last 100 chars.)

If there is indeed never a chance that a matched pattern breaches the first or the last 100 positions then the regex can be simplified to

/(.{$padding}) ($patt) (.{$padding})/x

Note that if a $patt is within the first/last $padding chars then this just won't match.

The program starts the regex engine for each of @patterns, what in principle may raise performance issues (not for one run with the tiny number of 4000 patterns, but such requirements tend to change and generally grow). But this is by far the simplest way to go since

  • we have no clue how the patterns may be distributed in the string, and

  • one match may be inside the 100-char buffer of another (we aren't told otherwise)

If there is a performance problem with this approach please update.


The input (and output) of the program can be organized in a better way using named command-line arguments via Getopt::Long, for an invocation like

program.pl --sequence seq.txt --search search.txt --padding 100

where each argument may be optional here, with defaults set in the file, and argument names may be shortened and/or given additional names, etc. Let me know if that is of interest

5
  • Thanks! Is there a way to input the strings from files rather than typing them? My real example has hundreds of thousands of characters for the string input
    – geri
    Nov 9, 2020 at 9:00
  • @geri Added reading from a file. In terms of robustness, there is a question of how to handle a possibility that there may be more than a line in a file; I simply break out of the loop after one read. There's other questions of how to handle possibly unexpected input, so please consider that for production use.
    – zdim
    Nov 9, 2020 at 9:13
  • @geri Forgot that you need patterns from a file as well (of course!), added
    – zdim
    Nov 9, 2020 at 9:41
  • Thank you! It does exactly what I asked without errors popping up.
    – geri
    Nov 9, 2020 at 11:41
  • @geri Great :) Let me know if questions come up
    – zdim
    Nov 9, 2020 at 19:30
2

One in awk. -v b=3 is the before context length -v a=3 is the after context length and -v n=3 is the match length which is always constant. It hashes all the substrings of seq.txt to memory so it uses it depending on the size of the seq.txt and you might want to follow the consumption with top, like: abcdefghij -> s["def"]="abcdefghi" , s["efg"]="bcdefghij" etc.

$ awk -v b=3 -v a=3 -v n=3 '
NR==FNR {
    e=length()-(n+a-1)
    for(i=1;i<=e;i++) {
        k=substr($0,(i+b),n)
        s[k]=s[k] (s[k]==""?"":ORS) substr($0,i,(b+n+a))
    }
    next
}
($0 in s) {
    print s[$0]
}' seq.txt search.txt

Output:

cdefghijk
mnopqrstu
rstuvwxyz
1
  • Any change you could test how much memory my solution uses? Nov 9, 2020 at 13:12
1

You can tell grep to search for all the patterns in one go.

sed 's/.*/.{0,100}&.{0,100}/' search.txt |
grep -zoEf - seq.txt |
sed G >do.grep

4000 patterns should be easy peasy, though if you get to hundreds of thousands, maybe you will want to optimize.

There is no Perl regex here, so I switched from the nonstandard grep -P to the POSIX-compatible and probably more efficient grep -E.

The surrounding context will consume any text it prints, so any match within 100 characters from the previous one will not be printed.

1

You can try following approach to your problem:

  • load string input data
  • load into an array patterns
  • loop through each pattern and look for it in the string
  • form an array from found matches
  • loop through matches array and print result

NOTE: the code is not tested due absence of input data

use strict;
use warnings;
use feature 'say';

my $fname_s = 'seq.txt';
my $fname_p = 'search.txt';
    
open my $fh, '<', $fname_s
    or die "Couldn't open $fname_s";
my $data = do { local $/; <$fh> };
close $fh;

open my $fh, '<', $fname_p
    or die "Couln't open $fname_p";
my @patterns = <$fh>;
close $fh;

chomp @patterns;

for ( @patterns ) {
    my @found = $data =~ s/(.{100}$_.{100})/g;
    s/(.{100})(.{50})(.{100})/$1 $2 $3/ && say for @found;
}

Test code for provided test data (added latter)

use strict;
use warnings;
use feature 'say';

my @pat  = qw/fgh pqr uvw/;
my $data = do { local $/; <DATA> }; 

for( @pat ) {
    say $1 if $data =~ /(.{3}$_.{3})/;
}

__DATA__
abcdefghijklmnopqrstuvwxyz

Output

cdefghijk
mnopqrstu
rstuvwxyz

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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