This is a test version of Biostars. For the public version, visit https://www.biostars.org.
find and count motif in FASTA file

Hi there,

I was looking up at some options to count a specific motif in a FASTA file: say, for instance, "TGGATATATGGACCTGTTTG". I came across a few options, some I knew already some I kind of put together:

  1. a simple 
    grep -c "TGGATATATGGACCTGTTTG" <file>.fasta
  2. an option with grep using -o, as follows: 
    grep -o "TGGATATATGGACCTGTTTG" <file>.fasta | wc -l
  3. a seqtk version following this post, as follows:
    seqkit locate -i -d -p TGGATATATGGACCTGTTTG -j 40 <file>.fasta | csvtk -t freq -k -f 1 -j 40 | awk -F'\t' '{sum+=$2;} END{print sum;}'
  4. a short Python script

python
from Bio import SeqIO
from Bio.SeqUtils import nt_search
total_count = 0
pattern = "TGGATATATGGACCTGTTTG"
for record in SeqIO.parse("/path/to/<file>.fasta", "fasta"):
    results = nt_search(str(record.seq), pattern)
    count = len(results) - 1
    total_count += count
    print(f"{record.id}: {count} occurrences")
print(f"Total occurrences: {total_count}")

Now, everything would be fine if not that these approaches give different results and I cannot really figure out why... with the 1st I get 718, the 2nd gives me 5143, the 3rd yields 8166 and, finally, the 4th again 5143.

I'm mostly interested in understanding why the two grep approaches have different outcomes, but also I'm interested in something like the 3rd or 4th approach that outputs a frequency count per FASTA header/block so that I can eventually plot the distribution of these motifs somehow.

P. S. I picked-up the Python from some of my old code and suggestions online, and I'm no expert; however, the print(f"{record.id}: {count} occurrences") returns a value of 0 for each header while the print(f"Total occurrences: {total_count}") prints out a total count. Why is that?

grep python seqtk
• 1,363 views • link
updated by lieven.sterck lieven.sterck • written by Matteo Ungaro • Matteo Ungaro

as Istvan Albert also pointed out:

Be VERY careful when using grep to do (biological) pattern matching!!

If you file is block-formatted any pattern that spans multiple lines will be missed. So always make sure your fasta file is single line formatted before running grep or use more dedicated tools to avoid this issue.

• 0 views • link
written by lieven.sterck lieven.sterck

2 answers

Everything that follows is on RHEL.

to count a specific motif in a FASTA file

Is that a strict match or are you going to allow for one or more errors?

Now, everything would be fine if not that these approaches give different results

Perhaps the issue in your case is line endings (unix/windows). Have you tried dos2unix?

When using a simple test file (below, with some sequences split on multiple lines deliberately ) I get the same answer for both option 1 and 2 and seqkit.

$ more motif.fa
>test_1
AGCTAGTGGATATATGGACCTGTTTGGATGCTCGCTAGC
>test_2
AGCTAGCTAGCTAGCTACGCTACGACTACGAT
>test_3
TGGATATATGGACCTGTTTGAGCTTTAGCTAGCTGGACT
CGCTAGCTAGACTGC
>test_4
GCGCTAGCTAGCTAGAGCTACGATTAGCTAGTC
>test_5
GAGCTAGCATCGACTTACGACTACGATCGACTAGCATCGACTACGGGCGC
CGATAGCTACGCTGGATATATGGACCTGTTTGGCTAGCTAGCGAC
CGATGATTATCTCTCTAC

With option 1:

$ grep -c TGGATATATGGACCTGTTTG motif.fa
3

With option 2:

$ grep -o TGGATATATGGACCTGTTTG motif.fa | wc -l
3

Another interesting grep option that you can try :

$ grep -n -C 2 TGGATATATGGACCTGTTTG motif.fa

Finally with seqkit I get the same answer.

$ seqkit locate -i -d -p TGGATATATGGACCTGTTTG -j 40 motif.fa 
seqID   patternName     pattern strand  start   end     matched
test_1  TGGATATATGGACCTGTTTG    TGGATATATGGACCTGTTTG    +       7       26      TGGATATATGGACCTGTTTG
test_3  TGGATATATGGACCTGTTTG    TGGATATATGGACCTGTTTG    +       1       20      TGGATATATGGACCTGTTTG
test_5  TGGATATATGGACCTGTTTG    TGGATATATGGACCTGTTTG    +       63      82      TGGATATATGGACCTGTTTG

With bbduk.sh from BBMap suite for fun (converting the motif into lower case characters):

$ bbduk.sh -Xmx4g in=motif.fa out=stdout.fa literal=TGGATATATGGACCTGTTTG k=7 kmask=lc

generates

Allocating kmer table:  0.023 seconds.
Initial:
Memory: max=4116m, total=4116m, free=3922m, used=194m

Added 14 kmers; time:   0.002 seconds.
Memory: max=4116m, total=4116m, free=3837m, used=279m

Input is being processed as unpaired
Started output streams: 0.013 seconds.
>test_1
AGCTAGtggatatatggacctgtttgGATGCTCGCTAGC
>test_2
AGCTAGCTAGCTAGCTACGCTACGACTACGAT
>test_3
tggatatatggacctgtttgAGCTTTAGCTAGCTGGACTCGCTAGCTAGACTGC
>test_4
GCGCTAGCTAGCTAGAGCTACGATTAGCTAGTC
>test_5
GAGCTAGCATCGACTTACGACTACGATCGACTAGCATCGACTACGGGCGCCGATAGCTACGCtggatata
tggacctgtttgGCTAGCTAGCGACCGATGATTATCTCTCTAC
Processing time:                0.012 seconds.

Input:                          5 reads                 271 bases.
KMasked:                        3 reads (60.00%)        60 bases (22.14%)
Total Removed:                  0 reads (0.00%)         0 bases (0.00%)
Result:                         5 reads (100.00%)       271 bases (100.00%)
• 0 views • link
written by GenoMax GenoMax

GenoMax indeed I'm looking for exact matches. Potentially, my issues is this one

Perhaps the issue in your case is line endings (unix/windows)?

I just overlooked at how the different commands could have behaved in respect to that. How can I fix it? Should I simply set all line endings to \n, or maybe attempt to remove spaces? Could this also explain why Python is returning 0s for all FASTA headers, and if not is there a way to understand what's happening with it?

Thanks again, much appreciated!

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

How can I fix it?

You can simply use dos2unix utility, which should be available on most unix distros. If not you could use (do not use -i if you want to preserve original files, redirect to new files).

For bash:

sed -i 's/\r$//' file.txt

or for macOS

sed -i '' $'s/\r$//' file.txt
• 0 views • link
written by GenoMax GenoMax

GenoMax thanks a lot, it worked. I'm still having some troubles with seqtk for some reason, but I think I can figure it out! One thing, about the Python what can be a way to have the actual number of occurrences printed correctly instead of having header: 0 occurrences from the:

print(f"{record.id}: {count} occurrences")
• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

If I understand what you are asking right, you could insert the following if statement as a check

if count > 0:
    print(f"{record.id}: {count} occurrences")

so nothing is printed if the count is 0.

So using your code should produce the following with the test file above.

$ python count.py 
test_1: 1 occurrences
test_3: 1 occurrences
test_5: 1 occurrences
Total occurrences: 3
• 0 views • link
written by GenoMax GenoMax

GenoMax once more, thanks again! This also fixed the issue with the Python attempt I made :)

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

First on the discrepancy. There is a big difference between counting the lines that match or the matches themselves:

echo AAA | grep  A
AAA

echo AAA | grep -o A
A
A
A

But in general don't use line oriented tools like grep to match motifs that may span over multiple lines. Use dedicated tools like the one GenoMax showed or as an alternative seqkit grep or seqkit locate

https://bioinf.shenwei.me/seqkit/usage/#grep

• 0 views • link
written by Istvan Albert Istvan Albert

Log in to answer this question.

More posts like this