Simple program to split DNA sequences based on ambiguous content

Question

The problem

I have a program that takes a file of DNA sequences as input and writes each entry to one of two output files: one file for sequences containing only canonical characters (A, C, G, and T), and another file for sequences containing non-canonical (or ambiguous) characters.

The data

DNA sequencing machines produce gobs of data, with each entry adhering to the following format.

@uniqueID free-form description and/or metadata
ACGTACGTACGTACGT
+
BBBBBBBBBBBBBBBB

The first line contains an ID and metadata, the second line contains DNA sequence (this is what the program filters on), the third line is stupid and can be ignored, and the fourth line contains ASCII-encoded quality values corresponding to each nucleotide in the DNA. Here are 3 entries from a real file.

@ERR899705.3060 HWI-D00574:77:C6A74ANXX:5:1101:2907:2432/1
GCCTTCCAGAAAAATGTTCCCAGTTCATTGATGGCAAACTCAAACAGGTTTGAAAAAAAATCTCAATATTAAATATGGATATGATGCATCAGTTCAAGGTGATTTCTGGTTAGCTGTGTGCTGACA
+
BCCCCGGGEGGGGGGGGGGGGGGGGGFGGGGGGGGGGGEGGGGGGGGE1>FCFGF1F<DDDGGGGDGGCGGGGGGGGGGGFGGGGGGEGEGG>FCEG@@0FGGGGGCEDGCG>DCGGGGG=0D=GG
@ERR899705.3061 HWI-D00574:77:C6A74ANXX:5:1101:2869:2437/1
CAAATAAATAAATAAATATTTAGTTACTGAAATAATAATGTTGATTGGTAAAATAAGAATGTATTCTCTCTAAAAACTAGTTAGTTTGTACTAAAAAATGGGCATAATTATTATCTCAAAGATTAG
+
BBBBBEGGGGGGGGGGGGGGGGGGGGGGD@GGGGGGGGEGC@FGFG>F:CGGGGGGFGGGCFCBDC11:F>FGGGGBEGGGGGGGGGC@F@FGGGGGG>000F@FGGGGGGEGGG00B0CFGGE>0
@ERR899705.3062 HWI-D00574:77:C6A74ANXX:5:1101:2786:2438/1
TTGGAAAACATTTATTACATTTTCTTCTAATTTCTTAGAATCATTTACTTGTAATATTTAGATACTCGGCGATCCTTTCAGGTATTCACTGCATATCCCCAAATATTGCATCAAACTCTCCCACGT
+
AABB0111@F111F111=1;@@;111E1?111:111=1=1111<11111<:11:C111<1EBG11:>F/////0/1:<F>C1:11=0:=0;@0>=0000==8>00000=8000<0000008000;0

When a DNA sequencing machine cannot reliably determine the molecule at a given position, it will use an N character instead of an ACGT character. Also, programs for analyzing DNA sequences will sometimes use other non-ACGT characters (IUPAC notation) to represent uncertainty about the content of DNA sequences.

The constraints

I need to handle sequences containing ambiguous characters separately from sequences containing only canonical characters. The input files are huge, so I want something as performant as possible, but I'm also trying to keep it as simple and easy to read as possible.

• The code does no sanity checking on the input data (assumes the number of lines in the input file is a multiple of 4, assumes there are no blank lines, etc). I am ok with this for my current purposes.
• It abuses stderr as a second output stream, although I'm not sure it's worth the extra clutter to parse command-line args, provide a usage statement, etc.

The code

The following minimal C++ program solves the problem as stated. Keeping the program simple, are there any ways I can improve performance or clarity of the code?

I'd welcome any feedback, but I'm specifically interested in the following.

• How should I format the if statement in the non_canon function? Everything I've tried looks hideous.
• I'm using the ternary operator to conditionally assign the output stream. My typical pattern in Python is to assign, check, and then reassign if needed, but std::ostream references cannot be reassigned in C++. Is there a clearer way of doing this?

---

#include <algorithm>
#include <iostream>
#include <string>

char non_canon(char bp)
{
    if (
        bp != 'A' && bp != 'a' && bp != 'C' && bp != 'c' &&
        bp != 'G' && bp != 'g' && bp != 'T' && bp != 't'
    ) {
        return true;
    }
    return false;
}

int main()
{
    std::string defline, seq, qual, buffer;
    while (getline(std::cin, defline)) {
        getline(std::cin, seq);
        getline(std::cin, buffer);
        getline(std::cin, qual);

        bool non_acgt = std::any_of(seq.begin(), seq.end(), non_canon);
        std::ostream& outstream = non_acgt ? std::cerr : std::cout;
        outstream << defline << '\n' << seq << "\n+\n" << qual << std::endl;
    }
    return 0;
}

Answer 1

Well one simple change will make it do lines in groups of four.

while (getline(std::cin, defline)
    && getline(std::cin, seq)
    && getline(std::cin, buffer)
    && getline(std::cin, qual)
      )
{
    // Only enter the loop if you correctly read all four lines.
}

We can also improve the test for a specific character. Note I am changing the return type to bool.

bool canon(char bp)
{
    bp = toupper(bp);
    return bp == 'A' || bp == 'C' || bp == 'G' || bp == 'T';
}

That will work. But its a bit slow if you get a lot of 'T'. We can optimize that if you want with a lookup table;

bool canon(unsigned char bp)
{
    static bool const is_cannon[] = {
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, // A C G
        0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // T
        0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, // a c g
        0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // t
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    return is_cannon[bp];
}

Bit of a micro optimization but if your files are huge. Now that @Edward has published an answer I like his switch statement better.

Prefer '\n' over std::endl. It (the endl) will only slow your code down. If forces a flush and this flush is unnecessary. The streams will automatically flush when they need to.

Your C++ code will run slower than C code because the C++ streams are tied to the C stream and try and keep them in sync. If you don't care about the C streams you can improve performance by unbinding them.

std::ios::sync_with_stdio(false);

Also because you are reading from the standard input this can cause slow downs (compared to reading from a file). Because the input needs to force the output to flush before getting reading input (so you see any messages that you should respond too). Since you don't want to wait for user input you can disable this behavior. In your case it probably will not make any difference but its worth adding just to be pedantic:

std::cin.tie(nullptr);

Questions

How should I format the if statement in the non_canon function? Everything I've tried looks hideous.

Don't use and if statement. My code above just returns the expression.

I'm using the ternary operator to conditionally assign the output stream. My typical pattern in Python is to assign, check, and then reassign if needed, but std::ostream references cannot be reassigned in C++. Is there a clearer way of doing this?

I think your code looks fine. But you don't need an intermediate variable.

(non_acgt ? std::cerr : std::cout) << defline << '\n' << seq << "\n+\n" << qual << "\n";

Or you can put it in a function.

 inline std::ostream& getOutputStream(bool non_acgt) {
     return non_acgt ? std::cerr : std::cout;
 }


 getOutputStream(non_acgt) << defline << '\n' << seq << "\n+\n" << qual << "\n";

Answer 2

The other review hits most of the important bits, but if you're writing software to classify DNA strings, I would bet that it's not the only processing you're ever going to need to do. With that in mind, I'd suggest that using C++ objects would not only make the code clearer, but provide a reusable object from which subsequent programs might be written to deal with the same data (e.g. searching for particular sequences within a sequence). Here's how that might look:

#include <algorithm>
#include <iostream>
#include <string>

class DNA {
public:
    friend std::istream& operator>>(std::istream &in, DNA &dna) {
        getline(in, dna.defline);
        getline(in, dna.seq);
        getline(in, dna.buffer);
        getline(in, dna.qual);
        return in;
    }
    friend std::ostream& operator<<(std::ostream &out, const DNA &dna) {
        return out << dna.defline << '\n' << dna.seq << "\n+\n" << dna.qual << '\n';
    }
    void classify(std::ostream &canon, std::ostream& noncanon) const {
        (isCanonical() ? canon : noncanon) << *this; 
    }
    bool isCanonical() const {
        return std::all_of(seq.begin(), seq.end(), is_canon);
    }
private:
    static bool is_canon(char bp) {
        switch(bp) {
            case 'A':
            case 'C':
            case 'G':
            case 'T':
            case 'a':
            case 'c':
            case 'g':
            case 't':
                return true;
        }
        return false;
    }

    std::string defline, seq, buffer, qual;
};

Now the main routine is extremely simple:

int main()
{
    std::ios::sync_with_stdio(false);
    for (DNA dna; std::cin >> dna; ) {
        // canonical sequences in the first file; others in the second
        dna.classify(std::cout, std::cerr);
    }
}

Note that by writing the character classification as a switch, we enable some compiler optimizations that aren't available in other ways (e.g. using a jump table). Also, by creating generic input and output stream operators, we can more clearly see the underlying logic.

Answer 3

In addition to the two excellent answers you have already received I would like to point out another way of simplifying non_canon:

bool non_canon(char bp){
    constexpr std::string canon_alphabet {"actgACTG"};
    return std::string::npos == canon_alphabet.find(bp);
}