Count the number of times a substring appears in a file and place it in a new column

Question

Question:

I have 2 files, file 1 is a TSV (BED) file that has 23 base-pair sequences in column 7, for example:

1 779692 779715 Sample_3 + 1  ATGGTGCTTTGTTATGGCAGCTC
1 783462 783485 Sample_4 - 1  ATGAATAAGTCAAGTAAATGGAC

File 2 is a FASTA file (hg19.fasta) that looks like this. Although it breaks across the lines, this continous string of A,C,G, and T's reprsents a continous string (i.e. a chromsosome). This file is the entire human reference genome build 19, so the two > headers followed by sequences essentially occurs 23 times for each of the 23 chromosomes:

>1 dna:chromosome chromosome:GRCh37:1:1:249250621:1
AATTTGACCAGAAGTTATGGGCATCCCTCCCCTGGGAAGGAGGCAGGCAGAAAAGTTTGGAATCTATGTAGTAAAATATG
TTACTCTTTTATATATATGAATAAGTCAAGTAAATGGACATACATATATGTGTGTATATGTGTATATATATATACACACA
TATATACATACATACATACATACATATTATCTGAATTAGGCCATGGTGCTTTGTTATGGCAGCTCTCTGGGATACATGTG
CAGAATGTACAGGTTTGTTACACAGGTATACACCTGCCATGGTTGTTTGCTGCACCCATCAACTCACCATCTACATTAGG
TATTTCTCCTAACGTTATCCCTCATGAATAAGTCAAGTAAATGGAC


>2 dna:chromosome chromosome:GRCh37:1:1:2492300:1

AATTTGACCAGAAGTTATGGGCATCCCTCCCCTGGGAAGGAGGCAGGCAGAAAAGTTTGGAATCTATGTAGTAAAATATG
TTACTCTTTTATATATATGAATAAGTCAAGTAAATGGACATACATATATGTGTGTATATGTGTATATATATATACACACA
TATATACATACATACATACATACATATTATCTGAATTAGGCCATGGTGCTTTGTTATGGCAGCTCTCTGGGATACATGTG

I want to 1) Find out how many times each 23bp sequence appears in the second file without overlapping any others including sequences that break across the lines and 2) append this number to a new column next to the sequence so the new file looks like this:

Desired Output:

1 779692 779715 Sample_3  + 1  ATGGTGCTTTGTTATGGCAGCTC 1 
1 783462 783485 Sample_4 - 1  ATGAATAAGTCAAGTAAATGGAC 2 

My attempt:

I imagine solving the first part will be some variation on grep, so far I've managed:

grep -o ATGGTGCTTTGTTATGGCAGCTC "$file_2" | grep -c ""

which gets the count of a specific sequence, but not each sequence in the column. I think appending the grep results will require awk and paste but I haven't gotten that far!

Any help is appreciated as always! =)

Updates and Edits:

• The actual size of these files is relatively massive (30mb or ~500,000 lines for each tsv/BED file) and the FASTA file is the entire human reference genome build 19, which is ~60,000,000 lines. The perl solution proposed by @choroba works, but doesn't scale well to these sizes.

• Unfortunately, because of the need to identify matches across the lines, the awk and bash/grep solutions memtnioned below won't work.

• I want multiple non-overlapping hits in the same chromosome to count as the actual number of hits. I.e. If you search for a sequence and get 2 hits in a single chromosome and 1 in another chromosome, the total count should be 3.

• Ted Lyngmo is very kindly helping me develop a solution in C++ that allows this to be run in a realistic timeframe, there's more detail on his post in this thread. And link to the Github for this is here =)

Answer 1

I'd reach for a programming language like Perl.

#!/usr/bin/perl
use warnings;
use strict;

my ($fasta_file, $bed_file) = @ARGV;
open my $fasta, '<', $fasta_file or die "$fasta_file: $!";
open my $bed,   '<', $bed_file   or die "$bed_file: $!";


my $seq;
while (<$fasta>) {
    $seq .= "\n", next if /^>/;
    chomp;
    $seq .= $_;
}

while (<$bed>) {
    chomp;
    my $short_seq = (split /\t/, $_)[-1];
    my $count = () = $seq =~ /\Q$short_seq\E/g;
    print "$_\t$count\n";
}

To count overlapping sequences, change the regex to a lookahead.

my $count = () = $seq =~ /(?=\Q$short_seq\E)/g;

Answer 2

If the second file is significantly bigger than the first one, I would try this awk script:

awk 'v==1 {a[$7];next}          # Get the pattern from first file into the array a
     v==2 {                     # For each line of the second file
       for(i in a){             # Loop though all patterns
         a[i]=split($0,b,i)-1   # Get the number of pattern match in the line
       }
     }
     v==3 {print $0,a[$7]}      # Re-read first file to add the number of pattern matches
 ' v=1 file1 v=2 file2 v=3 file1

Answer 3

Since grep -c seems to give you the correct count (matching lines, not counting multiple occurances on the same line) you could read the 7 fields from the TSV (BED) file and just print them again with the grep output added to the end:

#!/bin/bash

# read the fields into the array `v`:
while read -ra v
do
    # print the 7 first elements in the array + the output from `grep -c`:
    echo "${v[@]:0:7}" "$(grep -Fc "${v[6]}" hg19.fasta)"
done < tsv.bed > outfile

outfile will now contain

1 779692 779715 Sample_3 + 1 ATGGTGCTTTGTTATGGCAGCTC 1
1 783462 783485 Sample_4 - 1 ATGAATAAGTCAAGTAAATGGAC 2

Benchmarks

This table is a comparison of the three different solutions presented as answers here, with timings to finish different amount of tsv/bed records with the full hg19.fa file (excluding the records containing only N:s). hg19.fa contains 57'946'726 such records. As a baseline I've used two versions of a C++ program (called hgsearch/hgsearchmm). hgsearch reads the whole hg19.fa file into memory and then searches it in parallel. hgsearchmm uses a memory mapped file instead and then searches that (also in parallel).

search \ beds	1	2	100	1000	10000
awk	1m0.606s	17m19.899s	-	-	-
perl	13.263s	15.618s	4m48.751s	48m27.267s	-
bash/grep	2.088s	3.670s	3m27.378s	34m41.129s	-
hgsearch	8.776s	9.425s	30.619s	3m56.508s	38m43.984s
hgsearchmm	1.942s	2.146s	21.715s	3m28.265s	34m56.783s

The tests were run on an Intel Core i9 12 cores/24 HT:s in WSL/Ubuntu 20.04 (SSD disk). The sources for the scripts and baseline programs used can be found here