6

I have a file which contains a gene sequence such as:

ATGTGGATGGTGGGTTACAATGAAGGTGGTGAGTTCAACATGGCTGATTATCCATTCAGTGGAAGGAAACTAAGGCCTCTCATTCCAAGACCAGTCCCAGTCCCTACTACTTCTCCTAACAGCACTTCAACTATAACTCCTTCCTTAAACCGCATTCATGGTGGCAATGATTTATTTTCACAATATCATCACAATCTGCAGCAGCAAGCATCAGTAGGAGATCATAGCAAGAGATCAGAGTTGAATAATAATAATAATCCATCTGCAGCAGTTGTGGTGAGTTCAAGATGGAATCCAACACCAGAACAGTTAAGAGCACTGGAAGAATTGTATAGAAGAGGAACAAGAACACCTTCTGCTGAGCAAATCCAACAAATAACTGCCCAGCTTAGAAAATTTGGAAAAATTGAAGGCAAAAATGTTTTCTATTGGTTTCAGAATCACAAAGCCAGAGAAAGGCAAAAACGACGGCGTCAAATGGAATCAGCAGCTGCTGAGTTTGATTCTGCTATTGAAAAGAAAGACTTAGGCGCAAGTAGG


ACAGTGTTTGAAGTTGAACACACTAAAAACTGGCTACCATCTACAAATTCCAGTACCAGTACTCTTCATCTTGCAGAGGAATCTGTTTCAATTCAAAGGTCAGCAGCAGCAAAAGCAGATGGATGGCTCCAATTCGATGAAGCAGAATTACAGCAAAGAAGAAACTTTATGGAAAGGAATGCCACGTGGCATATGATGCAGTTAACTTCTTCTTGTCCTACAGCTAGCATGTCCACCACAACCACAGTAACAACTAGACTTATGGACCCAAAACTCATCAAGACCCATGAACTCAACTTATTCATTTCACCTCACACATACAAAGAAAGAGAAAACGCTTTTATCCACTTAAATACTAGTAGTACTCATCAAAATGAATCTGATCAAACCCTTCAACTTTTCCCAATAAGGAATGGAGATCATGGATGCACTGATCATCATCATCATCATCATAACATTATCAAAGAGACACAGATATCAGCTTCAGCAATCAATGCACCCAACCAGTTTATTGAGTTTCTTCCCTTGAAAAACTGA

I am trying to count the number of occurrence of "ATG" substring in the above string (which is only one line without line breaks.) My file contains tens (10s) of these sequences and I want to be able to count how many "ATG" in each sequence. Each sequence is separated from others by an empty line.

I tried grep but did not know which options I should use (if at all grep can solve the problem) and I googled for any awk example but I did not find any.

  • 1
    To clarify, did you want a total for each line or one total for the whole file? – John1024 Sep 22 '15 at 4:21
  • 1
    People have answered the question as stated, but I'm interested in clarifying the question, since this is DNA sequence, and you're searching for a particular codon: do you care which reading frame the codon is in? The answers so far all assume you do not. – Tim Cutts Sep 22 '15 at 7:38
  • @john1024 I would like a total for each line. – Faiz Lotfy Sep 24 '15 at 1:54
7

Returns the number of occurrences of ATG in each line:

awk -F'ATG' 'NF{print NF-1}' testfile

This works for files with one or many lines.

Example 1

Consider this test file:

$ cat testfile
xxATGxxATG

ATGxxxATGxxx

xxATGxxxxATGxxATGxx

The code correctly counts the occurrences of ATG:

$ awk -F'ATG' 'NF{print NF-1}' testfile
2
2
3

Example 2

Using the example in the current version of the question:

$ cat >file1
ATGTGGATGGTGGGTTACAATGAAGGTGGTGAGTTCAACATGGCTGATTATCCATTCAGTGGAAGGAAACTAAGGCCTCTCATTCCAAGACCAGTCCCAGTCCCTACTACTTCTCCTAACAGCACTTCAACTATAACTCCTTCCTTAAACCGCATTCATGGTGGCAATGATTTATTTTCACAATATCATCACAATCTGCAGCAGCAAGCATCAGTAGGAGATCATAGCAAGAGATCAGAGTTGAATAATAATAATAATCCATCTGCAGCAGTTGTGGTGAGTTCAAGATGGAATCCAACACCAGAACAGTTAAGAGCACTGGAAGAATTGTATAGAAGAGGAACAAGAACACCTTCTGCTGAGCAAATCCAACAAATAACTGCCCAGCTTAGAAAATTTGGAAAAATTGAAGGCAAAAATGTTTTCTATTGGTTTCAGAATCACAAAGCCAGAGAAAGGCAAAAACGACGGCGTCAAATGGAATCAGCAGCTGCTGAGTTTGATTCTGCTATTGAAAAGAAAGACTTAGGCGCAAGTAGG


ACAGTGTTTGAAGTTGAACACACTAAAAACTGGCTACCATCTACAAATTCCAGTACCAGTACTCTTCATCTTGCAGAGGAATCTGTTTCAATTCAAAGGTCAGCAGCAGCAAAAGCAGATGGATGGCTCCAATTCGATGAAGCAGAATTACAGCAAAGAAGAAACTTTATGGAAAGGAATGCCACGTGGCATATGATGCAGTTAACTTCTTCTTGTCCTACAGCTAGCATGTCCACCACAACCACAGTAACAACTAGACTTATGGACCCAAAACTCATCAAGACCCATGAACTCAACTTATTCATTTCACCTCACACATACAAAGAAAGAGAAAACGCTTTTATCCACTTAAATACTAGTAGTACTCATCAAAATGAATCTGATCAAACCCTTCAACTTTTCCCAATAAGGAATGGAGATCATGGATGCACTGATCATCATCATCATCATCATAACATTATCAAAGAGACACAGATATCAGCTTCAGCAATCAATGCACCCAACCAGTTTATTGAGTTTCTTCCCTTGAAAAACTGA

This results in:

$ awk -F'ATG' 'NF{print NF-1}' file1
9
15

How it works

awk implicitly loops through every line of a file. Each line is divided into fields.

  • -F'ATG'

    This tells awk to use ATG as the field separator.

  • NF{print NF-1}

    For each non-empty line, this tells awk to print the number of fields minus 1.

    (On empty lines, the number of fields, NF, is zero. So, the condition NF evaluates to false on these lines, effectively skipping over them.)

  • Alternatively awk '{print gsub("ATG","&")}' – 123 Sep 22 '15 at 8:32
3

From man grep

-o, --only-matching
    Print only the matched (non-empty) parts of a matching line, with each such part on a separate output line.

So you can try

$ grep -o 'ATG' file | wc -l
3

Improving tachomi's idea of grep -o, with John1024's testfile, here is a way to include both the count as well as indicating which line of the input file has that count, in case you are generating thousands of lines, and at a later point may need to retrace exactly which line a count came from.

Firstly a sample file, based on a slightly different version of John1024's testfile,

$ cat testfile2
xxATGxxATG

ATGxxxATGxxx

xxATGxxxxATGxxATG

If you add -n to show line number in original input, you see:

$ grep -no ATG testfile2
1:ATG
1:ATG
3:ATG
3:ATG
5:ATG
5:ATG
5:ATG

Finally, use uniq -c to do the counting:

$ grep -no ATG testfile2 | uniq -c
  2 1:ATG
  2 3:ATG
  3 5:ATG

So now you have both counts of ATG, as well as what line number (of the input file) that count came from.

You can also convert it to show only counts, by making use of awk:

$ grep -no ATG testfile2 | uniq -c | awk '{print $1}'
2
2
3

Awk's $1 refers to the first field.

1

You can use nl's section -delimiter to reset counts per set.

After copying your data to my clipboard...

xsel |
sed '/./!G;s/\n/::::::&::::/;s/ATG/&\
/g' | nl -d:: -v0

...yields...

 0  ATG
 1  TGGATG
 2  GTGGGTTACAATG
 3  AAGGTGGTGAGTTCAACATG
 4  GCTGATTATCCATTCAGTGGAAGGAAACTAAGGCCTCTCATTCCAAGACCAGTCCCAGTCCCTACTACTTCTCCTAACAGCACTTCAACTATAACTCCTTCCTTAAACCGCATTCATG
 5  GTGGCAATG
 6  ATTTATTTTCACAATATCATCACAATCTGCAGCAGCAAGCATCAGTAGGAGATCATAGCAAGAGATCAGAGTTGAATAATAATAATAATCCATCTGCAGCAGTTGTGGTGAGTTCAAGATG
 7  GAATCCAACACCAGAACAGTTAAGAGCACTGGAAGAATTGTATAGAAGAGGAACAAGAACACCTTCTGCTGAGCAAATCCAACAAATAACTGCCCAGCTTAGAAAATTTGGAAAAATTGAAGGCAAAAATG
 8  TTTTCTATTGGTTTCAGAATCACAAAGCCAGAGAAAGGCAAAAACGACGGCGTCAAATG
 9  GAATCAGCAGCTGCTGAGTTTGATTCTGCTATTGAAAAGAAAGACTTAGGCGCAAGTAGG


 0  ACAGTGTTTGAAGTTGAACACACTAAAAACTGGCTACCATCTACAAATTCCAGTACCAGTACTCTTCATCTTGCAGAGGAATCTGTTTCAATTCAAAGGTCAGCAGCAGCAAAAGCAGATG
 1  GATG
 2  GCTCCAATTCGATG
 3  AAGCAGAATTACAGCAAAGAAGAAACTTTATG
 4  GAAAGGAATG
 5  CCACGTGGCATATG
 6  ATG
 7  CAGTTAACTTCTTCTTGTCCTACAGCTAGCATG
 8  TCCACCACAACCACAGTAACAACTAGACTTATG
 9  GACCCAAAACTCATCAAGACCCATG
10  AACTCAACTTATTCATTTCACCTCACACATACAAAGAAAGAGAAAACGCTTTTATCCACTTAAATACTAGTAGTACTCATCAAAATG
11  AATCTGATCAAACCCTTCAACTTTTCCCAATAAGGAATG
12  GAGATCATG
13  GATG
14  CACTGATCATCATCATCATCATCATAACATTATCAAAGAGACACAGATATCAGCTTCAGCAATCAATG
15  CACCCAACCAGTTTATTGAGTTTCTTCCCTTGAAAAACTGA

sed just appends a newline after every occurrence of the string ATG it finds in input. sed also transforms blank-lines to look like:

::::::
::::

nl divides logical pages on its section -delimiter which is here set to ::. Three section -delimiters on a line by itself indicates the end of a logical page - which nl replaces with a blank-line - and two such the start of a new one. I start each count at -v0 zero because the last sequence in any chain will not contain the ATG string.

For less descriptive output, you can simply transform every ATG sequence to a character that doesn't otherwise appear in your input, and -delete the difference.

xsel |
sed 's/ATG/./g' | tr -dc . | wc -c

...which swaps every ATG in input with a dot, deletes all bytes which are not dots, then counts bytes:

24

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.