0

I have to write a script for an assignment which will take the file name as a command and output every unique 3 base codon in the file and how many times it occurs in descending order. The script has to check it has an argument and if not output an error message. I'm confused as to how to begin thee and the syntax for the regular expression which I will need to have the script search for each unique codon. The point is that a file is full of nothing but repeating letters of "a" "c" "t" "g" in random orders and the purpose is to write a script that will start from the beginning of the file and display each unique 3 letter codon and its count, for example aac or acg.

For example if a file named dnafile contains the pattern aacacgaactttaacacg then the script will take the following input and output

$script dnafile              
aac 3
acg 2
ttt 1

and if the error checking shows an error it should be script cannot open file dnafile for reading.

4
  • Sample input and output would be a good start...
    – Jeff Schaller
    Oct 28, 2018 at 20:15
  • What's a "DNS codon", and do I really need to know? Please consider rewriting the question in computing terms rather than bioinformatics terms. Otherwise you're restricting yourself from the set of people who might be able to answer to a considerably smaller group of people who understand what you're asking and can answer. Oct 28, 2018 at 21:48
  • @roima DNA codon, sorry typo. The point is that a file is full of nothing but repeating letters of "a" "c" "t" "g" in random orders and the purpose is to write a script that will start from the beginning of the file and display each unique 3 letter codon, for example aac or acg. But yes I will edit my question. Oct 28, 2018 at 21:52
  • Are there newlines in the file?
    – Kusalananda
    Oct 29, 2018 at 7:01

2 Answers 2

1

To get all the complete codons from a single file:

{
    for (pos = 1; pos < length - 1; ++pos) {
        codon = substr($0, pos, 3)
        print codon
    }
}

This short awk script simply outputs every three character long substring of each line. In bioinformatic speak, it outputs the complete codons for each frame in turn.

Testing on your given data:

$ awk -f script.awk <file
aac
aca
cac
acg
cga
gaa
aac
act
ctt
ttt
tta
taa
aac
aca
cac
acg

You may then sort and count the number of times that each codon occurs:

$ awk -f script.awk <file | sort | uniq -c
   3 aac
   2 aca
   2 acg
   1 act
   2 cac
   1 cga
   1 ctt
   1 gaa
   1 taa
   1 tta
   1 ttt

If the file contains newlines, then first remove these to be able to get the codons that are broken by newlines:

$ cat file
aacacgaactttaacacg
aacacgaactttaacacg
$ tr -d '\n' <file | awk -f script.awk | sort | uniq -c
   6 aac
   4 aca
   4 acg
   2 act
   4 cac
   3 cga
   2 ctt
   3 gaa
   2 taa
   2 tta
   2 ttt

(notice how the count for cga and gaa changed to 3 instead of to 2)

If your data is large, then you would have to come up with something else to handle the transition across newlines:

{
    $0 = lastbit $0

    for (pos = 1; pos < length - 1; ++pos) {
        codon = substr($0, pos, 3)
        print codon
    }

    lastbit = substr($0, length - 1)
}

This saves the last two bases from each line in lastbit and prepends them to the next line of sequence.

Running this on the same two lines of input as just above:

$ awk -f script.awk <file | sort | uniq -c
   6 aac
   4 aca
   4 acg
   2 act
   4 cac
   3 cga
   2 ctt
   3 gaa
   2 taa
   2 tta
   2 ttt

If you only want the codons in the first frame:

{
    $0 = lastbit $0

    for (pos = 1; pos < length - 1; ++pos) {
        if ((pos + length(lastbit)) % 3 == 1) {
            codon = substr($0, pos, 3)
            print codon
        }
    }

    lastbit = substr($0, length - 1)
}

Add an extra sort -nr step to the pipelines above if you want the output sorted on the counts in decreasing order.

See also:

0

One can do this on one line without an external file like so:

tr -d '\n' <dnafile |sed -e 's/\(...\)/\1\n/g' |sort |uniq -c |sed 's/^\s*\([0-9]*\) \([acgt]*\)$/\2 \1/'
aac 3
acg 2
ttt 1

If you don't need the columns in the order you gave, you can leave off the part after uniq -c. If you're doing this on a Mac, use gsed instead of sed.

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .