A classy DNA nucleobase counter

Question

This question is part of a series solving the Rosalind challenges. For the previous question in this series, see Counting nucleobases in a nucleotide c. The repository with all my up-to-date solutions so far can be found here.

---

I've been wondering for a while how awful it would look if I created an OO version of the first Rosalind challenge: DNA.

Turns out it isn't totally awful, so why not put it up for review?

I've considered adding more input sanitation this time, but there's not much that can go wrong. The sequence is derived from the argument provided, so it can't be empty or ArgParse will yell at you. I've forced it to capitalize every character if for some reason the user forgets to do so and I explicitly mention it's only supposed to work for DNA sequences (a version detecting whether a nucleobase is DNA/RNA or something different altogether would be a nice class for in a library, perhaps in a later version).

I've also considered (ab)using Biopython's Seq for this solution, but it would be overkill at this moment. Again, might be added later.

As said, it looks better than I expected, but it still looks like the method behind writing it is too chaotic. It could be more logically designed and probably more idiomatic as well. Argument handling, docstrings and internal datastructures are primary concern. Would it make sense to put argument handling in it's own class?

Problem: DNA

A string is simply an ordered collection of symbols selected from some alphabet and formed into a word; the length of a string is the number of symbols that it contains. An example of a length 21 DNA string (whose alphabet contains the symbols 'A', 'C', 'G', and 'T') is "ATGCTTCAGAAAGGTCTTACG."

Given:

A DNA string s of length at most 1000 nt.

Return:

Four integers (separated by spaces) counting the respective number of times that the symbols 'A', 'C', 'G', and 'T' occur in s.

Sample Dataset:

AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC

Sample Output:

20 12 17 21

Actual Dataset:

GTAGTATATCCACATACGAATTTATGGACGCTATCGCGACGCTCTGAACGCGGATTCCGGCAAGCCATGCAGTGTCTACTTACACCCGACGTCCACTCCGCCGGATACATGAGACATATCCCGACAGCTACACGATACTCTATACCTGATCGTACTGTTCTAGCGTACAATACGCGTGTCGCCGGCGGTATCTAAGCTTATAAGTTTCCGCCCTGCGGTTGCGATATTCCGTGTGACGGGCGGCAAATACGGGCCGGGGCACTGTAAGACTCGGAGGGTAAACAGTGGTACAGTGTTGATTATCCGAAACCGCCGGGGAGAGGTGTCGATGACGCTGCAACCTGTACCGCAGTTAGGGGCCGGGATGATGTGGTCGTGCACGCCTTAATCTCCGACTAACTGCAGTTCATAGAGCTTTGATTACATTTCCCAGAGGAATGGATACCTAGTCTTAAAGACAGCACGCTACGCACCAAAATCAGCACCAAACGCAAGCCTCGTGTTCAGACAACTTATTCGCCGTGAAGCTAAGTTTAAGAGGACTCGGATACTATGGAAAATCCTTCTTACAGCAGTGCCATTCCAGTTCGTACTTGTCAGCGAGGCGTTCAGCAATTGGCACCGCAATCATGGGGCACCACGTCATCGTAACTTGTTGGGACACAGTCACGCGGGGTCGTAGTTCTACCGATGGTTAATGCCAGTACCTCAGACGGTCTACACCCTTCGCCCCGGAGAGCGGATCAATGCATCTAGGTAGGTATTCTGCATCCCCCAAGCCTGCGTAGGAACAGATGGATCAGTCCGCCCTACTATGCGTTTAGTGAAGGACGAGGCATTAATACTAGTCCCGTGCGGCCCTGATGACCTGTCTTACGCGATGTACGCCTGACGTCGTGCATTGAATGTTTCATCTTAGGGTCGTAT

Actual Output:

223 244 238 222

Of-course, plenty of other input works correctly as well.

gtagtatatccacatacgaatttatggacgctatcgcgacgctctgaacgcggattccggcaagccatgcagtgtctacttacacccgacgtccactccgccggatacatgagacatatcccgacagctacacgatactctatacctgatcgtactgttctagcgtacaatacgcgtgtcgccggcggtatctaagcttataagtttccgccctgcggttgcgatattccgtgtgacgggcggcaaatacgggccggggcactgtaagactcggagggtaaacagtggtacagtgttgattatccgaaaccgccggggagaggtgtcgatgacgctgcaacctgtaccgcagttaggggccgggatgatgtggtcgtgcacgccttaatctccgactaactgcagttcatagagctttgattacatttcccagaggaatggatacctagtcttaaagacagcacgctacgcaccaaaatcagcaccaaacgcaagcctcgtgttcagacaacttattcgccgtgaagctaagtttaagaggactcggatactatggaaaatccttcttacagcagtgccattccagttcgtacttgtcagcgaggcgttcagcaattggcaccgcaatcatggggcaccacgtcatcgtaacttgttgggacacagtcacgcggggtcgtagttctaccgatggttaatgccagtacctcagacggtctacacccttcgccccggagagcggatcaatgcatctaggtaggtattctgcatcccccaagcctgcgtaggaacagatggatcagtccgccctactatgcgtttagtgaaggacgaggcattaatactagtcccgtgcggccctgatgacctgtcttacgcgatgtacgcctgacgtcgtgcattgaatgtttcatcttagggtcgtat

The above provides the same output and the following (invalid) input will fail as expected:

AACTAIN


Invalid DNA nucleobase encountered. Aborting.

Usage: python DNA.py <sequence>

DNA.py:

from argparse import ArgumentParser


class DnaNucleobaseSequence:
    """
    Contains a DNA nucleobase sequence.
    Counts the occurence of the available bases.
    """
    def __init__(self, sequence):
        self.sequence = sequence.upper()
        self.counted_nucleobases = self.count_nucleobases()

    def count_nucleobases(self):
        """
        Return the count of nucleobases in sequence.

        Since this is explicitly a DNA sequence, only adenine, cytosine,
        guanine and thymine are allowed to be present.

        Exit on KeyError (invalid nucleobase(s)).
        """
        count = {
            "A": 0,
            "C": 0,
            "G": 0,
            "T": 0,
        }
        try:
            for nucleobase in self.sequence:
                count[nucleobase] = count[nucleobase] + 1
            return count
        except KeyError:
            print("Invalid DNA nucleobase encountered. Aborting.")
            raise SystemExit

    def space_separated_nucleobases(self):
        """
        Return the count of A, C, G and T separated by a single space.
        """
        return "{0} {1} {2} {3}".format(
            self.counted_nucleobases["A"],
            self.counted_nucleobases["C"],
            self.counted_nucleobases["G"],
            self.counted_nucleobases["T"]
        )


def main():
    parser = ArgumentParser(description='DNA nucleobase sequence')
    parser.add_argument(
        'sequence',
        type=str,
        help='Should contain a DNA nucleobase sequence'
        )

    args = parser.parse_args()

    seq = DnaNucleobaseSequence(args.sequence)
    print(seq.space_separated_nucleobases())

if __name__ == "__main__":
    main()

Answer 1

• I don't really like self.sequence, it's only used in count_nucleobases, and doesn't seem to be needed as an attribute to the class.

• I'm not a fan of raise SystemExit. First if I'm reading things correctly, it causes the program to have an exit status of 0. Which means that it ran without errors, which is not true here. Secondly, I think of it more as something main should raise. Say you have a main that wraps all your code in a try except, which may display an error message if you set debugging on, but otherwise causes the program to have an exit status of 1. Such as most Linux programs.

  And so I don't think it's right to put it in the middle of a class.

• I'd rather use count[nucleobase] += 1, it's just a short hand, but I prefer it.

• You may prefer to use collections.Counter. Personally I think it makes the class as it is obsolete.

  However, it does change your code to consume the entire sequence before erroring.

• You can index objects when using format strings. And so you can simplify space_separated_nucleobases to:

  "{0[A]} {0[C]} {0[G]} {0[T]}".format(self.counted_nucleobases)
  

• If I were to do this, I'd remove the KeyError code, and move all the code to the main function.

And so I'd change your class to:

class DnaNucleobaseSequence:
    """
    Contains a DNA nucleobase sequence.
    Counts the occurence of the available bases.
    """
    def __init__(self, sequence):
        self.counted_nucleobases = self.count_nucleobases(sequence.upper())

    def count_nucleobases(self, sequence):
        """
        Return the count of nucleobases in sequence.

        Since this is explicitly a DNA sequence, only adenine, cytosine,
        guanine and thymine are allowed to be present.

        Exit on KeyError (invalid nucleobase(s)).
        """
        count = Counter(sequence)
        if set(count) - set('ACGT'):
            raise KeyError("Invalid DNA nucleobase encountered. Aborting.")
        return count

    def space_separated_nucleobases(self):
        """
        Return the count of A, C, G and T separated by a single space.
        """
        return "{0[A]} {0[C]} {0[G]} {0[T]}".format(self.counted_nucleobases)

Answer 2

@Peilonrayz already touched a bit about that but I'd like to comment on the use of that SystemExit.

Instead of raising it manually, you should consider using sys.exit which is both clearer (as it directly remind the C call) and more versatile: feed it an integer and it's the return code of the program. Feed it an object and it's __str__ gets printed to stderr and the return code is 1.