Extract preceding and trailing characters to a matched string from file in awk

Viewed 251

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!

4 Answers

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

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

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.

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
Related