test/cigar.cc

Code
Comments
Other
Rev Date Author Line
3200 03 May 14 peter 1 // $Id$
3200 03 May 14 peter 2
3200 03 May 14 peter 3 /*
4207 26 Aug 22 peter 4   Copyright (C) 2014, 2019, 2020, 2022 Peter Johansson
3200 03 May 14 peter 5
3200 03 May 14 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3200 03 May 14 peter 7
3200 03 May 14 peter 8   The yat library is free software; you can redistribute it and/or
3200 03 May 14 peter 9   modify it under the terms of the GNU General Public License as
3200 03 May 14 peter 10   published by the Free Software Foundation; either version 3 of the
3200 03 May 14 peter 11   License, or (at your option) any later version.
3200 03 May 14 peter 12
3200 03 May 14 peter 13   The yat library is distributed in the hope that it will be useful,
3200 03 May 14 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3200 03 May 14 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3200 03 May 14 peter 16   General Public License for more details.
3200 03 May 14 peter 17
3200 03 May 14 peter 18   You should have received a copy of the GNU General Public License
3200 03 May 14 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3200 03 May 14 peter 20 */
3200 03 May 14 peter 21
3200 03 May 14 peter 22 #include <config.h>
3200 03 May 14 peter 23
3200 03 May 14 peter 24 #include "Suite.h"
3200 03 May 14 peter 25
3200 03 May 14 peter 26 #include "yat/utility/Aligner.h"
4003 15 Oct 20 peter 27 #include "yat/utility/Exception.h"
3200 03 May 14 peter 28 #include "yat/utility/Matrix.h"
4003 15 Oct 20 peter 29 #include "yat/utility/utility.h"
3200 03 May 14 peter 30
4127 14 Jan 22 peter 31 #include <cassert>
3200 03 May 14 peter 32 #include <sstream>
3200 03 May 14 peter 33 #include <string>
3200 03 May 14 peter 34 #include <vector>
3200 03 May 14 peter 35
3200 03 May 14 peter 36 using namespace theplu::yat;
3200 03 May 14 peter 37
3200 03 May 14 peter 38 int main(int argc, char* argv[])
3200 03 May 14 peter 39 {
3200 03 May 14 peter 40   test::Suite suite(argc, argv);
3200 03 May 14 peter 41
3200 03 May 14 peter 42   std::string ref = "XXXXThis is a test string, oh yes!!!!!!!!!!!!!!";
3200 03 May 14 peter 43   std::string seq = ref.substr(3,17);
3200 03 May 14 peter 44   seq[4] = 'G';
3200 03 May 14 peter 45   seq.erase(10,1);
3200 03 May 14 peter 46
3200 03 May 14 peter 47   suite.out() << "ref: " << ref << "\n";
3200 03 May 14 peter 48   suite.out() << "seq: " << seq << "\n";
3200 03 May 14 peter 49
3200 03 May 14 peter 50   utility::Matrix m(ref.size(), seq.size(), -3);
3200 03 May 14 peter 51   for (size_t i=0; i<m.rows(); ++i)
3200 03 May 14 peter 52     for (size_t j=0; j<m.columns(); ++j)
3200 03 May 14 peter 53       if (ref[i] == seq[j])
3200 03 May 14 peter 54         m(i,j) = 1;
3200 03 May 14 peter 55
3200 03 May 14 peter 56   utility::Matrix alignment_score(m.rows()+1, m.columns()+1, 0);
3200 03 May 14 peter 57   utility::Aligner aligner(2,3);
3200 03 May 14 peter 58   aligner(m, alignment_score);
3200 03 May 14 peter 59
3200 03 May 14 peter 60   size_t besti=0;
3200 03 May 14 peter 61   size_t bestj=0;
3200 03 May 14 peter 62   for (size_t i=0; i<alignment_score.rows(); ++i)
3200 03 May 14 peter 63     for (size_t j=0; j<alignment_score.columns(); ++j)
3200 03 May 14 peter 64       if (alignment_score(i,j) > alignment_score(besti, bestj)) {
3200 03 May 14 peter 65         besti = i;
3200 03 May 14 peter 66         bestj = j;
3200 03 May 14 peter 67       }
3200 03 May 14 peter 68
3200 03 May 14 peter 69
3200 03 May 14 peter 70   suite.out() << "best i: " << besti << "\n";
3200 03 May 14 peter 71   suite.out() << "best j: " << bestj << "\n";
3200 03 May 14 peter 72
3200 03 May 14 peter 73   utility::Aligner::Cigar cig = aligner.cigar(besti, bestj);
3200 03 May 14 peter 74   std::ostringstream os;
3200 03 May 14 peter 75   os << cig;
3200 03 May 14 peter 76   suite.out() << os.str() << "\n";
3200 03 May 14 peter 77   if (os.str()!="10M1D6M") {
3200 03 May 14 peter 78     suite.add(false);
3200 03 May 14 peter 79     suite.err() << "error: incorrect string\n";
3200 03 May 14 peter 80   }
3840 13 Aug 19 peter 81   for (utility::Aligner::Cigar::const_iterator c=cig.begin();c!=cig.end();++c){
3840 13 Aug 19 peter 82     std::string str(BAM_CIGAR_STR);
3840 13 Aug 19 peter 83     assert(*c < str.size());
3840 13 Aug 19 peter 84     suite.out() << str[*c];
3840 13 Aug 19 peter 85   }
3840 13 Aug 19 peter 86   suite.out() << "\n";
4003 15 Oct 20 peter 87   utility::Aligner::Cigar cig2(os.str());
4003 15 Oct 20 peter 88   if (std::distance(cig.begin(), cig.end())!=
4003 15 Oct 20 peter 89       std::distance(cig2.begin(), cig2.end()) ||
4003 15 Oct 20 peter 90       !std::equal(cig.begin(), cig.end(), cig2.begin())) {
4003 15 Oct 20 peter 91     suite.add(false);
4003 15 Oct 20 peter 92     suite.err() << "error: incorrect construction from Cigar(string)\n";
4003 15 Oct 20 peter 93     suite.err() << ": " << cig2 << "\nexpected:" << cig << "\n";
4003 15 Oct 20 peter 94   }
4003 15 Oct 20 peter 95   try {
4003 15 Oct 20 peter 96     utility::Aligner::Cigar cig3("BANANA");
4003 15 Oct 20 peter 97     suite.add(false);
4003 15 Oct 20 peter 98     suite.err() << "error: Cigar(\"BANANA\") did not throw\n";
4003 15 Oct 20 peter 99     suite.err() << "cigar: " << cig3 << "\n";
4003 15 Oct 20 peter 100   }
4003 15 Oct 20 peter 101   catch (utility::invalid_argument& e) {
4003 15 Oct 20 peter 102     suite.out() << "expected exception: ";
4003 15 Oct 20 peter 103     utility::print_what(e, suite.out());
4003 15 Oct 20 peter 104     suite.out() << "\n";
4003 15 Oct 20 peter 105   }
3200 03 May 14 peter 106
3213 05 May 14 peter 107   if (!suite.add(cig.length()==17))
3213 05 May 14 peter 108     suite.err() << "error: length: " << cig.length() << " expected 17\n";
3213 05 May 14 peter 109   if (!suite.add(cig.query_length()==16))
3213 05 May 14 peter 110     suite.err() << "error: query length: " << cig.query_length()
3213 05 May 14 peter 111                 << " expected 16\n";
3213 05 May 14 peter 112   if (!suite.add(cig.reference_length()==17))
3213 05 May 14 peter 113     suite.err() << "error: length: " << cig.reference_length()
3213 05 May 14 peter 114                 << " expected 17\n";
3213 05 May 14 peter 115
3201 03 May 14 peter 116   if (!suite.add(cig.op(0)==0))
3201 03 May 14 peter 117     suite.err() << "error: cig.op(0): " << cig.op(0) << " expected 0\n";
3201 03 May 14 peter 118   suite.out() << cig << "\n";
3201 03 May 14 peter 119   if (!suite.add(cig.opchr(0)=='M'))
3201 03 May 14 peter 120     suite.err() << "error: cig.opchr(0): " << cig.opchr(0) << " expected M\n";
3201 03 May 14 peter 121   if (!suite.add(cig.oplen(0)==10))
3201 03 May 14 peter 122     suite.err() << "error: cig.oplen(0): " << cig.oplen(0) << " expected 10\n";
3200 03 May 14 peter 123   cig.pop_front(2);
3201 03 May 14 peter 124   suite.out() << cig << "\n";
3201 03 May 14 peter 125   if (!suite.add(cig.oplen(0)==8))
3201 03 May 14 peter 126     suite.err() << "error: cig.oplen(0): " << cig.oplen(0) << " expected 8\n";
3201 03 May 14 peter 127
3200 03 May 14 peter 128   cig.pop_front(9);
3201 03 May 14 peter 129   suite.out() << cig << "\n";
3201 03 May 14 peter 130   if (!suite.add(cig.oplen(0)==6))
3201 03 May 14 peter 131     suite.err() << "error: cig.oplen(1): " << cig.oplen(1) << " expected 6\n";
3201 03 May 14 peter 132
3200 03 May 14 peter 133   cig.push_back(3);
3201 03 May 14 peter 134   suite.out() << cig << "\n";
3200 03 May 14 peter 135   cig.pop_back(1);
3201 03 May 14 peter 136   suite.out() << cig << "\n";
3201 03 May 14 peter 137
3201 03 May 14 peter 138   if (!suite.add(cig.oplen(0)==6))
3201 03 May 14 peter 139     suite.err() << "error: cig.oplen(1): " << cig.oplen(1) << " expected 6\n";
3201 03 May 14 peter 140
3200 03 May 14 peter 141   cig[0];
3200 03 May 14 peter 142   suite.out() << "size: " << cig.size() << "\n";
3200 03 May 14 peter 143
3200 03 May 14 peter 144   cig.clear();
3200 03 May 14 peter 145   return suite.return_value();
3200 03 May 14 peter 146 }

Generated on Sat Oct 14 02:20:45 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)