3

I am using the following grep script to output all the unmatched patterns:

grep -oFf patterns.txt large_strings.txt | grep -vFf - patterns.txt > unmatched_patterns.txt

patterns file contains the following 12-characters long substrings (some instances are shown below):

6b6c665d4f44
8b715a5d5f5f
26364d605243
717c8a919aa2

large_strings file contains extremely long strings of around 20-100 million characters longs (a small piece of the string is shown below):

121b1f212222212123242223252b36434f5655545351504f4e4e5056616d777d80817d7c7b7a7a7b7c7d7f8997a0a2a2a3a5a5a6a6a6a6a6a7a7babbbcbebebdbcbcbdbdbdbdbcbcbcbcc2c2c2c2c2c2c2c2c4c4c4c3c3c3c2c2c3c3c3c3c3c3c3c3c2c2c1c0bfbfbebdbebebebfbfc0c0c0bfbfbfbebebdbdbdbcbbbbbababbbbbcbdbdbdbebebfbfbfbebdbcbbbbbbbbbcbcbcbcbcbcbcbcbcb8b8b8b7b7b6b6b6b8b8b9babbbbbcbcbbbabab9b9bababbbcbcbcbbbbbababab9b8b7b6b6b6b6b7b7b7b7b7b7b7b7b7b7b6b6b5b5b6b6b7b7b7b7b8b8b9b9b9b9b9b8b7b7b6b5b5b5b5b5b4b4b3b3b3b6b5b4b4b5b7b8babdbebfc1c1c0bfbec1c2c2c2c2c1c0bfbfbebebebebfc0c1c0c0c0bfbfbebebebebebebebebebebebebebdbcbbbbbab9babbbbbcbcbdbdbdbcbcbbbbbbbbbbbabab9b7b6b5b4b4b4b4b3b1aeaca9a7a6a9a9a9aaabacaeafafafafafafafafafb1b2b2b2b2b1b0afacaaa8a7a5a19d9995939191929292919292939291908f8e8e8d8c8b8a8a8a8a878787868482807f7d7c7975716d6b6967676665646261615f5f5e5d5b5a595957575554525

How can we speed up the above script (gnu parallel, xargs, fgrep, etc.)? I tried using --pipepart and --block but it doesn't allow you to pipe two grep commands.

Btw these are all hexadecimal strings and patterns.

3 Answers 3

3

This should work:

parallel --pipepart --block -1 -a large_strings.txt grep -oFf patterns.txt |
  grep -vFf - patterns.txt > unmatched_patterns.txt

If you have ripgrep use that:

parallel --pipepart --block -1 -a large_strings.txt rg -oFf patterns.txt |
  rg -vFf - patterns.txt > unmatched_patterns.txt

If patterns.txt is huge, too, have a look at:

https://www.gnu.org/software/parallel/man.html#EXAMPLE:-Grepping-n-lines-for-m-regular-expressions

Your situation is also pretty close to the problem that BLAT solves, except BLAT is built for DNA. But I don't see that you could not use BLAT for your situation - possibly with a few changes (well, you could convert every hex value to 2 DNA letters and use it directly). BLAT is as fast as a lookup in a database, so there is no comparison to grep. http://genome.ucsc.edu/FAQ/FAQblat.html#blat3

7
  • when I try to run the script it gives me an error stating the following: Died at /home/bin/parallel line 288. Jan 25, 2021 at 0:45
  • @user3441801 The only version I can find a die in line 288 is in version 20160622. Can you upgrade to 20210122? See: git.savannah.gnu.org/cgit/parallel.git/tree/README
    – Ole Tange
    Jan 25, 2021 at 8:16
  • yes, you're right. Apparently, I was using an older version of parallel within the Miniconda environment. Switching to the right version of parallel fixed this issue. Your code works much much faster than my original code. Thank you. Jan 25, 2021 at 17:49
  • also using ripgrep is incredibly faster and accurate than the traditional grep Jan 25, 2021 at 17:50
  • Is there a way to add options --block -1 --lb to patterns.txt file as well. What I mean is that is there a way to add options --block -1 --lb to both the large_strings.txt and patterns.txt. My patterns file (19GB) is much larger than the large_strings file (9GB), @Ole Tange. Thanks. Jan 27, 2021 at 1:53
1

A much more efficient answer that does not use grep:

build_k_mers() {
    k="$1"
    slot="$2"
    perl -ne 'for $n (0..(length $_)-'"$k"') {                                                                                               
       $prefix = substr($_,$n,2);                                                                                                            
       $fh{$prefix} or open $fh{$prefix}, ">>", "tmp/kmer.$prefix.'"$slot"'";                                                                
       $fh = $fh{$prefix};                                                                                                                   
       print $fh substr($_,$n,'"$k"'),"\n"                                                                                                   
    }'
}
export -f build_k_mers

rm -rf tmp
mkdir tmp
export LC_ALL=C
# search strings must be sorted for comm                                                                                                     
parsort patterns.txt | awk '{print >>"tmp/patterns."substr($1,1,2)}' &

# make shorter lines: Insert \n(last 12 char before \n) for every 32k                                                                         
# This makes it easier for --pipepart to find a newline                                                                                      
# It will not change the kmers generated                                                                                                     
perl -pe 's/(.{32000})(.{12})/$1$2\n$2/g' large_strings.txt > large_lines.txt
# Build 12-mers                                                                                                                              
parallel --pipepart --block -1 -a large_lines.txt 'build_k_mers 12 {%}'
# -j10 and 20s may be adjusted depending on hardware
parallel -j10 --delay 20s 'parsort -u tmp/kmer.{}.* > tmp/kmer.{}; rm tmp/kmer.{}.*' ::: `perl -e 'map { printf "%02x ",$_ } 0..255'`
wait
parallel comm -23 {} {=s/patterns./kmer./=} ::: tmp/patterns.??

I have tested this on a full job (patterns.txt: 9GBytes/725937231 lines, large_strings.txt: 19GBytes/184 lines) and on my 64-core machine it completes in 3 hours.

6
  • This is great, @Ole. Thanks for your kind efforts. When I try to execute the code it runs for about 7 mins and then finishes. All it does is that it sorts the patterns file successfully along with the following error: Quantifier in {,} bigger than 32766 in regex; marked by <-- HERE in m/(.{ <-- HERE 64000})(.{12})/ at -e line 1. /var/spool/slurm/slurmd/job*/slurm_script: line 23: pv: command not found Also, I can't install the PV package on the cluster as I can't seem to find the original source. Jan 28, 2021 at 16:58
  • I managed to install pv -- found the original build source (github.com/icetee/pv), @Ole. I executed the code again and this time it took around 4 mins to sort a 9GB patterns file giving the same following message: Quantifier in {,} bigger than 32766 in regex; marked by <-- HERE in m/(.{ <-- HERE 64000})(.{12})/ at -e line 1. For some reason, it doesn't pass the make shorter lines function. Can we modify this code to also output "matched" patterns too, please? Thanks. Jan 28, 2021 at 17:52
  • @user3441801 See edit.
    – Ole Tange
    Jan 28, 2021 at 19:22
  • It worked this time and it's so much faster (took 20 mins for the job to finish), @Ole Tange. I have a cluster of four nodes each with 32 cores. Changing from 64000 to 32000 worked. But the size of the output file is significantly larger (30GB) in comparison to the output file (which is around 23MB) generated using rg script: parallel --pipepart --block -1 -a large_strings.txt rg -oFf patterns.txt | rg -vFf - patterns.txt > unmatched_patterns.txt I used a smaller test patterns file of 25MB and large_strings file of 7GB so that I can compare the results with grep. Am I missing something? Jan 28, 2021 at 23:43
  • 1
    @user3441801 Yes, that should work.
    – Ole Tange
    Jan 29, 2021 at 17:40
0

No idea if this will work...

while read line; 
do 
  grep -oF "$line" large_strings.txt & 
done <patterns.txt | grep -vFf - patterns.txt > unmatched_patterns.txt
3
  • This may lead to race conditions leading to wrong results: The inner greps may spit out lines at the same time and thus you may get half a line from one grep followed by half-a-line from a different grep.
    – Ole Tange
    Jan 24, 2021 at 22:55
  • @mcede when I try to run the while loop I get the following error: -bash: fork: retry: Resource temporarily unavailable Jan 25, 2021 at 0:54
  • @ole I agree with your comment Jan 25, 2021 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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