source: trunk/test/alignment.cc @ 2815

Last change on this file since 2815 was 2815, checked in by Peter, 10 years ago

closes #706. Implement a generalization of SmithWaternan? and NeedlemanWunsch? in a new class utility::Aligner. Class also provide a way to access not only alignment score but also which elements were aligned (to maximize the score). Old function SmithWaterman? did not provide this functionality and is deprecated for new member function Aligner::smith_waterman.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.4 KB
Line 
1// $Id: alignment.cc 2815 2012-08-28 02:38:36Z peter $
2
3/*
4  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen
6  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
7
8  This file is part of the yat library, http://dev.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 3 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with yat. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#include "Suite.h"
25
26#include "yat/utility/Aligner.h"
27#include "yat/utility/Alignment.h"
28
29#include "yat/utility/Matrix.h"
30
31#include <gsl/gsl_cdf.h>
32
33#include <cassert>
34#include <cmath>
35#include <fstream>
36#include <iostream>
37#include <sstream>
38#include <string>
39#include <vector>
40
41using namespace theplu::yat;
42
43void read_vector(std::vector<double>& v,std::string& str)
44{
45  std::istringstream s(str);
46  while (!s.eof()) {
47    double d;
48    s >> d;
49    if (!s.fail())
50      v.push_back(d);
51  }
52}
53
54std::istream& operator>>(std::istream& s,
55                         std::vector<std::vector<double> >& m)
56{
57  while (!s.eof()) {
58    std::string row;
59    char ch;
60    while (((ch=s.get())!='\n') && s.good())
61      row+=ch;
62    if (row.size() || s.good()) { // if stream good, read an empty peak set
63      std::vector<double> v;
64      read_vector(v,row);
65      m.push_back(v);
66    }
67  }
68  return s;
69}
70
71double match(const double x, const double y, const double s)
72{ 
73  return 2*gsl_cdf_gaussian_Q(std::abs(x-y),sqrt(2)*s); 
74}
75 
76
77double score(const std::vector<double>& l1,
78             const std::vector<double>& l2,
79             const double sigma, 
80             theplu::yat::utility::Matrix& dot_matrix,
81             std::vector<std::pair<size_t,size_t> >& path)
82{
83  dot_matrix.resize(l1.size(),l2.size());
84  for (size_t i=0; i<l1.size(); i++) 
85    for (size_t j=0; j<l2.size(); j++) {
86      assert(i<dot_matrix.rows());
87      assert(j<dot_matrix.columns());
88      dot_matrix(i,j)=match(l1[i],l2[j],sigma);
89    }
90  return theplu::yat::utility::NeedlemanWunsch(dot_matrix,path,0);
91}
92
93
94int main(int argc, char* argv[])
95{
96  test::Suite suite(argc, argv);
97
98  std::ifstream s(test::filename("data/isoform.peaks").c_str());
99  std::vector<std::vector<double> > peaksets;
100  s >> peaksets;
101
102  for (size_t i=0; i<peaksets.size()-1; i++) 
103    for (size_t j=i+1; j<peaksets.size(); j++) {
104      utility::Matrix dot_m;
105      std::vector<std::pair<size_t,size_t> > path;
106      score(peaksets[i], peaksets[j], 1.0, dot_m, path);
107    }
108
109  std::string a("AGGUUGUCCGUGGUGAGUUCGCA");
110  std::string b("GAGGUUGUCCGUGGUGAGUUCG");
111  utility::Matrix m(a.size(), b.size());
112  for (size_t j=0; j<a.size(); ++j)
113    for (size_t k=0; k<b.size(); ++k)
114      m(j,k) = (a[j]==b[k] ? 1 : -1000);
115  double score = utility::Aligner(100, 100).smith_waterman(m);
116  if (score!=21)
117    suite.add(false);
118
119  // testing ssearch
120  if (utility::ssearch("Hello", "Hll", 0.0, 1.0)!=2){
121    suite.err() << "aligning 'Hello' and 'Hll' gives score " 
122                << utility::ssearch("Hello", "Hll", 0.0, 1.0)
123                << " expected " << 2 << std::endl;
124    suite.add(false);
125  }
126  if (utility::ssearch("Hello", "Peter said you can't say 'allo", 1, 1)!=3)
127    suite.add(false);
128
129  return suite.return_value();
130}
Note: See TracBrowser for help on using the repository browser.