test/bam_multi_region_iterator.cc

Code
Comments
Other
Rev Date Author Line
4272 24 Jan 23 peter 1 // $Id$
4272 24 Jan 23 peter 2 //
4272 24 Jan 23 peter 3 // Copyright (C) 2023 Peter Johansson
4272 24 Jan 23 peter 4 //
4272 24 Jan 23 peter 5 // This program is free software; you can redistribute it and/or modify
4272 24 Jan 23 peter 6 // it under the terms of the GNU General Public License as published by
4272 24 Jan 23 peter 7 // the Free Software Foundation; either version 3 of the License, or
4272 24 Jan 23 peter 8 // (at your option) any later version.
4272 24 Jan 23 peter 9 //
4272 24 Jan 23 peter 10 // This program is distributed in the hope that it will be useful, but
4272 24 Jan 23 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
4272 24 Jan 23 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
4272 24 Jan 23 peter 13 // General Public License for more details.
4272 24 Jan 23 peter 14 //
4272 24 Jan 23 peter 15 // You should have received a copy of the GNU General Public License
4272 24 Jan 23 peter 16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
4272 24 Jan 23 peter 17
4272 24 Jan 23 peter 18 #include <config.h>
4272 24 Jan 23 peter 19
4272 24 Jan 23 peter 20 #include "Suite.h"
4272 24 Jan 23 peter 21
4272 24 Jan 23 peter 22 #ifdef YAT_HAVE_HTSLIB
4272 24 Jan 23 peter 23 #include "yat/omic/BamRead.h"
4272 24 Jan 23 peter 24 #include "yat/omic/BamReadIterator.h"
4272 24 Jan 23 peter 25 #include "yat/omic/BamRegion.h"
4272 24 Jan 23 peter 26 #endif
4272 24 Jan 23 peter 27
4272 24 Jan 23 peter 28 #include <yat/utility/Exception.h>
4272 24 Jan 23 peter 29
4272 24 Jan 23 peter 30 #include <algorithm>
4272 24 Jan 23 peter 31 #include <cassert>
4272 24 Jan 23 peter 32 #include <string>
4272 24 Jan 23 peter 33
4272 24 Jan 23 peter 34 using namespace theplu::yat;
4272 24 Jan 23 peter 35
4272 24 Jan 23 peter 36 void test1(test::Suite& suite);
4272 24 Jan 23 peter 37 int main(int argc, char* argv[])
4272 24 Jan 23 peter 38 {
4272 24 Jan 23 peter 39   test::Suite suite(argc, argv, true);
4272 24 Jan 23 peter 40
4272 24 Jan 23 peter 41   try {
4272 24 Jan 23 peter 42     test1(suite);
4272 24 Jan 23 peter 43   }
4272 24 Jan 23 peter 44   catch (std::runtime_error& e) {
4272 24 Jan 23 peter 45     suite.err() << "what: " << e.what() << "\n";
4272 24 Jan 23 peter 46     suite.add(false);
4272 24 Jan 23 peter 47   }
4272 24 Jan 23 peter 48
4272 24 Jan 23 peter 49   return suite.return_value();
4272 24 Jan 23 peter 50 }
4272 24 Jan 23 peter 51
4272 24 Jan 23 peter 52
4272 24 Jan 23 peter 53 void test1(test::Suite& suite)
4272 24 Jan 23 peter 54 {
4272 24 Jan 23 peter 55 #ifdef YAT_HAVE_HTSLIB
4272 24 Jan 23 peter 56
4272 24 Jan 23 peter 57   using namespace omic;
4272 24 Jan 23 peter 58   std::string fn = "../../data/foo.sorted.bam";
4272 24 Jan 23 peter 59
4272 24 Jan 23 peter 60   InBamFile file(fn);
4272 24 Jan 23 peter 61   if (!file.is_open()) {
4272 24 Jan 23 peter 62     suite.err() << "cannot read " << fn << "\n";
4272 24 Jan 23 peter 63     suite.add(false);
4272 24 Jan 23 peter 64     return;
4272 24 Jan 23 peter 65   }
4272 24 Jan 23 peter 66
4272 24 Jan 23 peter 67   std::vector<BamRegion> regions =
4272 24 Jan 23 peter 68     { BamRegion(1, 24300, 24400), BamRegion(1, 24500, 24550),
4272 24 Jan 23 peter 69       BamRegion(1, 24600, 24700), BamRegion(1, 24650, 24750)};
4272 24 Jan 23 peter 70
4272 24 Jan 23 peter 71   BamReadIterator  first(file, regions);
4272 24 Jan 23 peter 72   BamReadIterator end;
4272 24 Jan 23 peter 73
4272 24 Jan 23 peter 74   // copy reads from regions into a vector
4272 24 Jan 23 peter 75   std::vector<BamRead> reads(first, end);
4272 24 Jan 23 peter 76   if (reads.empty()) {
4272 24 Jan 23 peter 77     suite.add(false);
4272 24 Jan 23 peter 78     suite.err() << "read vector is empty\n";
4272 24 Jan 23 peter 79     return;
4272 24 Jan 23 peter 80   }
4272 24 Jan 23 peter 81   suite.out() << "number of reads: " << reads.size() << "\n";
4272 24 Jan 23 peter 82   for (const auto& bam : reads) {
4272 24 Jan 23 peter 83     bool overlap = false;
4272 24 Jan 23 peter 84     for (const auto& reg : regions) {
4272 24 Jan 23 peter 85       if (!(bam.tid() != reg.tid() ||
4272 24 Jan 23 peter 86             bam.pos() >= reg.end() ||
4272 24 Jan 23 peter 87             bam.end() <= reg.begin()))
4272 24 Jan 23 peter 88         overlap = true;
4272 24 Jan 23 peter 89     }
4272 24 Jan 23 peter 90     if (!overlap) {
4272 24 Jan 23 peter 91       suite.add(false);
4272 24 Jan 23 peter 92       suite.err() << "error: unexpected read position: " << bam.pos()
4272 24 Jan 23 peter 93                   << " : " << bam.end() << "\n";
4272 24 Jan 23 peter 94     }
4272 24 Jan 23 peter 95   }
4272 24 Jan 23 peter 96
4272 24 Jan 23 peter 97   BamLessPos less_pos;
4272 24 Jan 23 peter 98   // check that range is sorted
4272 24 Jan 23 peter 99   if (!std::is_sorted(reads.begin(), reads.end(), less_pos)) {
4272 24 Jan 23 peter 100     suite.add(false);
4272 24 Jan 23 peter 101     suite.err() << "range is not sorted\n";
4272 24 Jan 23 peter 102   }
4272 24 Jan 23 peter 103
4272 24 Jan 23 peter 104   for (auto it = reads.begin(); it!=reads.end(); ++it) {
4272 24 Jan 23 peter 105     for (auto it2 = it + 1; it2!=reads.end() && !less_pos(*it, *it2); ++it2) {
4272 24 Jan 23 peter 106       assert(it->pos() == it2->pos());
4272 24 Jan 23 peter 107       if (same_query_name(*it, *it2) && it->flag() == it2->flag()) {
4272 24 Jan 23 peter 108         suite.add(false);
4272 24 Jan 23 peter 109         suite.err() << "duplicated read: " << it->pos() << "\t"
4272 24 Jan 23 peter 110                     << it->name() << "\n";
4272 24 Jan 23 peter 111       }
4272 24 Jan 23 peter 112     }
4272 24 Jan 23 peter 113   }
4272 24 Jan 23 peter 114
4272 24 Jan 23 peter 115 #endif
4272 24 Jan 23 peter 116 }

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