baseComp.cc

#include <Sequence/Fasta.hpp>
#include <Sequence/stateCounter.hpp>
#include <Sequence/SeqUtilities.hpp>
#include <algorithm>
#include <iostream>
#include <functional>

using namespace std;
using namespace Sequence;

//an example of using iterators for objects of type Sequence::Seq in conjunction with
//STL algorithms

//usage: cat file | baseComp

struct outputFreq
    : public std::binary_function<std::pair<char, unsigned>, unsigned, void>
{
    inline void
    operator()(const std::pair<char, unsigned> &p, const unsigned &u) const
    {
        std::cout << p.first << " (" << double(p.second) / double(u) << ") ";
    }
};

int
main()
{
    Fasta x;
    while (!cin.fail())
        {
            cin >> x;

            //method 1--use Sequence::stateCounter
            //count base composition for the sequence using the STL algorithm for_each
            //note the use of the iterators belonging to class Sequence::Seq
            //class Sequence::stateCounter is used to keep track of the counts of each base
            stateCounter count
                = for_each(x.begin(), x.end(), stateCounter('-'));
            cout << "Base composition for sequence " << x.name << "\n";

            if (count.ndna == false)
                {
                    unsigned len
                        = x.length() - count.gap
                          - count.n; //ungapped length,exclude missing data

                    //turn the counts into percentages
                    double percA
                        = (count.a > 0) ? double(count.a) / double(len) : 0.;
                    double percG
                        = (count.g > 0) ? double(count.g) / double(len) : 0.;
                    double percC
                        = (count.c > 0) ? double(count.c) / double(len) : 0.;
                    double percT
                        = (count.t > 0) ? double(count.t) / double(len) : 0.;

                    //output
                    std::cout << "using Sequence::stateCounter: "
                              << "A (" << percA << ") "
                              << "G (" << percG << ") "
                              << "C (" << percC << ") "
                              << "T (" << percT << ")\n";
                }
            else
                {
                    cout << "non-DNA character encountered.  Skipping...\n";
                }
            //method 2--use Sequence::makeCountList (doesn't check for non-standard DNA characters)
            std::map<char, unsigned> m
                = Sequence::makeCountList(x.begin(), x.end());
            std::cout << "using Sequence::makeCountList(): ";
            std::for_each(m.begin(), m.end(),
                          std::bind2nd(outputFreq(), x.length()));
            std::cout << std::endl;
        }
}