test/bam_iterator.cc

Code
Comments
Other
Rev Date Author Line
2895 10 Dec 12 peter 1 // $Id$
2895 10 Dec 12 peter 2 //
3999 08 Oct 20 peter 3 // Copyright (C) 2012, 2013, 2014, 2017, 2020 Peter Johansson
2895 10 Dec 12 peter 4 //
2895 10 Dec 12 peter 5 // This program is free software; you can redistribute it and/or modify
2895 10 Dec 12 peter 6 // it under the terms of the GNU General Public License as published by
2895 10 Dec 12 peter 7 // the Free Software Foundation; either version 3 of the License, or
2895 10 Dec 12 peter 8 // (at your option) any later version.
2895 10 Dec 12 peter 9 //
2895 10 Dec 12 peter 10 // This program is distributed in the hope that it will be useful, but
2895 10 Dec 12 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
2895 10 Dec 12 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
2895 10 Dec 12 peter 13 // General Public License for more details.
2895 10 Dec 12 peter 14 //
2895 10 Dec 12 peter 15 // You should have received a copy of the GNU General Public License
2895 10 Dec 12 peter 16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
2895 10 Dec 12 peter 17
2895 10 Dec 12 peter 18 #include <config.h>
2895 10 Dec 12 peter 19
2895 10 Dec 12 peter 20 #include "Suite.h"
2895 10 Dec 12 peter 21
3883 24 Mar 20 peter 22 #ifdef YAT_HAVE_HTSLIB
2895 10 Dec 12 peter 23 #include "yat/omic/BamFile.h"
2895 10 Dec 12 peter 24 #include "yat/omic/BamRead.h"
2895 10 Dec 12 peter 25 #include "yat/omic/BamReadIterator.h"
2895 10 Dec 12 peter 26 #include "yat/omic/BamWriteIterator.h"
2943 04 Jan 13 peter 27 #endif
2895 10 Dec 12 peter 28
2895 10 Dec 12 peter 29 #include <algorithm>
2895 10 Dec 12 peter 30 #include <cassert>
2895 10 Dec 12 peter 31 #include <fstream>
2895 10 Dec 12 peter 32 #include <string>
2895 10 Dec 12 peter 33
2895 10 Dec 12 peter 34 using namespace theplu::yat;
2943 04 Jan 13 peter 35
2943 04 Jan 13 peter 36 void test1(test::Suite& suite);
2943 04 Jan 13 peter 37
2943 04 Jan 13 peter 38 int main(int argc, char* argv[])
2943 04 Jan 13 peter 39 {
3201 03 May 14 peter 40   test::Suite suite(argc, argv, true);
2983 04 Feb 13 peter 41   try {
2983 04 Feb 13 peter 42     test1(suite);
2983 04 Feb 13 peter 43   }
2983 04 Feb 13 peter 44   catch (std::runtime_error& e) {
2983 04 Feb 13 peter 45     suite.add(false);
2983 04 Feb 13 peter 46     suite.err() << "error: what: " << e.what() << "\n";
2983 04 Feb 13 peter 47   }
2943 04 Jan 13 peter 48   return suite.return_value();
2943 04 Jan 13 peter 49 }
2943 04 Jan 13 peter 50
3883 24 Mar 20 peter 51 #ifdef YAT_HAVE_HTSLIB
2895 10 Dec 12 peter 52 using namespace omic;
2895 10 Dec 12 peter 53
2895 10 Dec 12 peter 54 struct BamReadEqual
2895 10 Dec 12 peter 55 {
2895 10 Dec 12 peter 56   bool operator()(const BamRead& lhs, const BamRead& rhs) const
2895 10 Dec 12 peter 57   {
2895 10 Dec 12 peter 58     // FIXME
2895 10 Dec 12 peter 59     /*
2895 10 Dec 12 peter 60     if (lhs.l_aux() != rhs.l_aux())
2895 10 Dec 12 peter 61       return false;
2895 10 Dec 12 peter 62     if (lhs.data_len() != rhs.data_len())
2895 10 Dec 12 peter 63       return false;
2895 10 Dec 12 peter 64     if (lhs.m_data() != rhs.m_data())
2895 10 Dec 12 peter 65       return false;
2895 10 Dec 12 peter 66     for (int i=0; i<lhs.data_len(); ++i)
2895 10 Dec 12 peter 67       if (lhs.data()[i] != rhs.data()[i])
2895 10 Dec 12 peter 68         return false;
2895 10 Dec 12 peter 69     */
2895 10 Dec 12 peter 70
2895 10 Dec 12 peter 71     // comparing core
2895 10 Dec 12 peter 72     if (lhs.tid() != rhs.tid())
2895 10 Dec 12 peter 73       return false;
2895 10 Dec 12 peter 74     if (lhs.pos() != rhs.pos())
2895 10 Dec 12 peter 75       return false;
2986 18 Feb 13 peter 76     if (lhs.mtid() != rhs.mtid())
2986 18 Feb 13 peter 77       return false;
2986 18 Feb 13 peter 78     if (lhs.mpos() != rhs.mpos())
2986 18 Feb 13 peter 79       return false;
2895 10 Dec 12 peter 80     if (lhs.core().bin != rhs.core().bin)
2895 10 Dec 12 peter 81       return false;
2895 10 Dec 12 peter 82     if (lhs.core().qual != rhs.core().qual)
2895 10 Dec 12 peter 83       return false;
2895 10 Dec 12 peter 84     if (lhs.core().l_qname != rhs.core().l_qname)
2895 10 Dec 12 peter 85       return false;
2895 10 Dec 12 peter 86     if (lhs.core().flag != rhs.core().flag)
2895 10 Dec 12 peter 87       return false;
2895 10 Dec 12 peter 88     if (lhs.core().n_cigar != rhs.core().n_cigar)
2895 10 Dec 12 peter 89       return false;
2895 10 Dec 12 peter 90     if (lhs.core().l_qseq != rhs.core().l_qseq)
2895 10 Dec 12 peter 91       return false;
2895 10 Dec 12 peter 92     if (lhs.core().mtid != rhs.core().mtid)
2895 10 Dec 12 peter 93       return false;
2895 10 Dec 12 peter 94     if (lhs.core().mpos != rhs.core().mpos)
2895 10 Dec 12 peter 95       return false;
2895 10 Dec 12 peter 96     if (lhs.core().isize != rhs.core().isize)
2895 10 Dec 12 peter 97       return false;
2895 10 Dec 12 peter 98
2895 10 Dec 12 peter 99     return true;
2895 10 Dec 12 peter 100   }
2895 10 Dec 12 peter 101 };
2895 10 Dec 12 peter 102
2943 04 Jan 13 peter 103
2943 04 Jan 13 peter 104 void test1(test::Suite& suite)
2895 10 Dec 12 peter 105 {
2895 10 Dec 12 peter 106   std::string file = "../../data/foo.sorted.bam";
2895 10 Dec 12 peter 107
2895 10 Dec 12 peter 108   InBamFile in;
2895 10 Dec 12 peter 109   in.open(file);
2895 10 Dec 12 peter 110   std::string outname("bam_iterator.bam");
2895 10 Dec 12 peter 111   OutBamFile out(outname, in.header());
3698 25 Sep 17 peter 112   BamWriteIterator default_constructed;
3698 25 Sep 17 peter 113   test::avoid_compiler_warning(default_constructed);
3150 23 Dec 13 peter 114   BamWriteIterator bam_out(out);
3150 23 Dec 13 peter 115   std::copy(BamReadIterator(in), BamReadIterator(), bam_out);
3150 23 Dec 13 peter 116   if (false) // avoid running; only compile
3150 23 Dec 13 peter 117     suite.test_output_iterator(bam_out);
2895 10 Dec 12 peter 118   out.close();
2895 10 Dec 12 peter 119
2895 10 Dec 12 peter 120   // compare that copy has equally many elements
2895 10 Dec 12 peter 121   in.close();
2895 10 Dec 12 peter 122   in.open(file);
2895 10 Dec 12 peter 123   size_t n1 = std::distance(BamReadIterator(in), BamReadIterator());
2943 04 Jan 13 peter 124   suite.out() << "entries in " << file << ": " << n1 << "\n";
2895 10 Dec 12 peter 125   in.close();
2895 10 Dec 12 peter 126
2895 10 Dec 12 peter 127   InBamFile copy(outname);
2895 10 Dec 12 peter 128   size_t n2 = std::distance(BamReadIterator(copy), BamReadIterator());
2895 10 Dec 12 peter 129   copy.close();
2943 04 Jan 13 peter 130   suite.out() << "entries in " << outname << ": " << n2 << "\n";
2895 10 Dec 12 peter 131   if (n1!=n2) {
2943 04 Jan 13 peter 132     suite.out() << "error: not same number of entries in copy\n";
2943 04 Jan 13 peter 133     suite.add(false);
2895 10 Dec 12 peter 134   }
2895 10 Dec 12 peter 135
2895 10 Dec 12 peter 136   copy.open(outname);
2895 10 Dec 12 peter 137   in.open(file);
2895 10 Dec 12 peter 138   if (!std::equal(BamReadIterator(in), BamReadIterator(),
2895 10 Dec 12 peter 139                   BamReadIterator(copy), BamReadEqual())) {
2943 04 Jan 13 peter 140     suite.out() << "error: " << file << " and " << outname << " not equal\n";
2943 04 Jan 13 peter 141     suite.add(false);
2895 10 Dec 12 peter 142   }
2895 10 Dec 12 peter 143   copy.close();
2895 10 Dec 12 peter 144   in.close();
2895 10 Dec 12 peter 145   unlink(outname.c_str());
3150 23 Dec 13 peter 146
2895 10 Dec 12 peter 147 }
3209 04 May 14 peter 148 #else
3209 04 May 14 peter 149 void test1(test::Suite& suite) {}
2943 04 Jan 13 peter 150 #endif