libsequence  1.9.5
baseComp.cc
1 #include <Sequence/Fasta.hpp>
4 #include <algorithm>
5 #include <iostream>
6 #include <functional>
7 
8 using namespace std;
9 using namespace Sequence;
10 
13 //an example of using iterators for objects of type Sequence::Seq in conjunction with
14 //STL algorithms
15 
16 //usage: cat file | baseComp
17 
18 struct outputFreq
19  : public std::binary_function<std::pair<char, unsigned>, unsigned, void>
20 {
21  inline void
22  operator()(const std::pair<char, unsigned> &p, const unsigned &u) const
23  {
24  std::cout << p.first << " (" << double(p.second) / double(u) << ") ";
25  }
26 };
27 
28 int
29 main()
30 {
31  Fasta x;
32  while (!cin.fail())
33  {
34  cin >> x;
35 
36  //method 1--use Sequence::stateCounter
37  //count base composition for the sequence using the STL algorithm for_each
38  //note the use of the iterators belonging to class Sequence::Seq
39  //class Sequence::stateCounter is used to keep track of the counts of each base
40  stateCounter count
41  = for_each(x.begin(), x.end(), stateCounter('-'));
42  cout << "Base composition for sequence " << x.name << "\n";
43 
44  if (count.ndna == false)
45  {
46  unsigned len
47  = x.length() - count.gap
48  - count.n; //ungapped length,exclude missing data
49 
50  //turn the counts into percentages
51  double percA
52  = (count.a > 0) ? double(count.a) / double(len) : 0.;
53  double percG
54  = (count.g > 0) ? double(count.g) / double(len) : 0.;
55  double percC
56  = (count.c > 0) ? double(count.c) / double(len) : 0.;
57  double percT
58  = (count.t > 0) ? double(count.t) / double(len) : 0.;
59 
60  //output
61  std::cout << "using Sequence::stateCounter: "
62  << "A (" << percA << ") "
63  << "G (" << percG << ") "
64  << "C (" << percC << ") "
65  << "T (" << percT << ")\n";
66  }
67  else
68  {
69  cout << "non-DNA character encountered. Skipping...\n";
70  }
71  //method 2--use Sequence::makeCountList (doesn't check for non-standard DNA characters)
72  std::map<char, unsigned> m
73  = Sequence::makeCountList(x.begin(), x.end());
74  std::cout << "using Sequence::makeCountList(): ";
75  std::for_each(m.begin(), m.end(),
76  std::bind2nd(outputFreq(), x.length()));
77  std::cout << std::endl;
78  }
79 }
bool operator()(const std::pair< key, value > &l, const std::pair< key, value > &r) const
keep track of state counts at a site in an alignment or along a sequence
FASTA sequence stream.
Definition: Fasta.hpp:49
STL namespace.
The namespace in which this library resides.
Declaration of Sequence::makeCountList (an alternative to Sequence::stateCounter), Sequence::internalGapCheck.
Declaration of Sequence::Fasta streams.
declaration of Sequence::stateCounter, a class to keep track of nucleotide counts either at a site in...
std::map< typename std::iterator_traits< Iterator >::value_type, unsigned > makeCountList(Iterator beg, Iterator end)