test/bam.cc

Code
Comments
Other
Rev Date Author Line
2972 26 Jan 13 peter 1 // $Id$
2972 26 Jan 13 peter 2 //
3999 08 Oct 20 peter 3 // Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2020 Peter Johansson
2972 26 Jan 13 peter 4 //
2972 26 Jan 13 peter 5 // This program is free software; you can redistribute it and/or modify
2972 26 Jan 13 peter 6 // it under the terms of the GNU General Public License as published by
2972 26 Jan 13 peter 7 // the Free Software Foundation; either version 3 of the License, or
2972 26 Jan 13 peter 8 // (at your option) any later version.
2972 26 Jan 13 peter 9 //
2972 26 Jan 13 peter 10 // This program is distributed in the hope that it will be useful, but
2972 26 Jan 13 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
2972 26 Jan 13 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
2972 26 Jan 13 peter 13 // General Public License for more details.
2972 26 Jan 13 peter 14 //
2972 26 Jan 13 peter 15 // You should have received a copy of the GNU General Public License
2972 26 Jan 13 peter 16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
2972 26 Jan 13 peter 17
2972 26 Jan 13 peter 18 #include <config.h>
2972 26 Jan 13 peter 19
2972 26 Jan 13 peter 20 #include "Suite.h"
2972 26 Jan 13 peter 21
3883 24 Mar 20 peter 22 #ifdef YAT_HAVE_HTSLIB
2972 26 Jan 13 peter 23 #include "yat/omic/BamFile.h"
2972 26 Jan 13 peter 24 #include "yat/omic/BamRead.h"
2972 26 Jan 13 peter 25 #endif
2972 26 Jan 13 peter 26
3295 25 Jul 14 peter 27 #include "yat/utility/Aligner.h"
3674 01 Aug 17 peter 28 #include "yat/utility/FileUtil.h"
3674 01 Aug 17 peter 29 #include "yat/utility/utility.h"
3295 25 Jul 14 peter 30
2972 26 Jan 13 peter 31 #include <cassert>
2972 26 Jan 13 peter 32 #include <string>
2972 26 Jan 13 peter 33
2972 26 Jan 13 peter 34 using namespace theplu::yat;
2972 26 Jan 13 peter 35
2972 26 Jan 13 peter 36 void test1(test::Suite& suite);
2972 26 Jan 13 peter 37
2972 26 Jan 13 peter 38 int main(int argc, char* argv[])
2972 26 Jan 13 peter 39 {
3201 03 May 14 peter 40   test::Suite suite(argc, argv, true);
3883 24 Mar 20 peter 41 #ifdef YAT_HAVE_HTSLIB
2972 26 Jan 13 peter 42   test1(suite);
2972 26 Jan 13 peter 43 #endif
2972 26 Jan 13 peter 44   return suite.return_value();
2972 26 Jan 13 peter 45 }
2972 26 Jan 13 peter 46
3883 24 Mar 20 peter 47 #ifdef YAT_HAVE_HTSLIB
2972 26 Jan 13 peter 48 using namespace omic;
2972 26 Jan 13 peter 49
3395 23 Mar 15 peter 50 void test_alignment_score(test::Suite& suite, const BamRead& b,
3395 23 Mar 15 peter 51                           const std::string& ref, double correct)
3395 23 Mar 15 peter 52 {
3395 23 Mar 15 peter 53   int score = alignment_score(b, ref.begin(), 3, 2, 5);
3395 23 Mar 15 peter 54   if (!suite.equal(score, correct)) {
3395 23 Mar 15 peter 55     suite.add(false);
3395 23 Mar 15 peter 56     suite.err() << "error: alignment_score: " << score
3395 23 Mar 15 peter 57                 << " expected " << correct << "\n";
3395 23 Mar 15 peter 58     suite.err() << "query: " << b.sequence() << "\n";
3395 23 Mar 15 peter 59     suite.err() << "cigar: " << b.cigar_str() << "\n";
3395 23 Mar 15 peter 60     suite.err() << "ref:   " << ref << "\n";
3395 23 Mar 15 peter 61   }
3395 23 Mar 15 peter 62 }
3395 23 Mar 15 peter 63
3395 23 Mar 15 peter 64
3395 23 Mar 15 peter 65 void test_alignment_score(test::Suite& suite, const BamRead& b)
3395 23 Mar 15 peter 66 {
3395 23 Mar 15 peter 67   std::string reference = b.sequence();
3395 23 Mar 15 peter 68   test_alignment_score(suite, b, reference, 100);
3395 23 Mar 15 peter 69   reference[48] = 'N';
3395 23 Mar 15 peter 70   test_alignment_score(suite, b, reference, 96);
3395 23 Mar 15 peter 71
3395 23 Mar 15 peter 72   utility::Aligner::Cigar cig;
3395 23 Mar 15 peter 73   cig.push_back(BAM_CEQUAL, 48);
3395 23 Mar 15 peter 74   cig.push_back(BAM_CDIFF, 1);
3395 23 Mar 15 peter 75   cig.push_back(BAM_CEQUAL, 51);
3395 23 Mar 15 peter 76   BamRead bam(b);
3395 23 Mar 15 peter 77   bam.cigar(cig);
3395 23 Mar 15 peter 78   test_alignment_score(suite, bam, reference, 96);
3395 23 Mar 15 peter 79
3395 23 Mar 15 peter 80   cig.clear();
3395 23 Mar 15 peter 81   cig.push_back(BAM_CSOFT_CLIP, 10);
3395 23 Mar 15 peter 82   cig.push_back(BAM_CEQUAL, 38);
3395 23 Mar 15 peter 83   cig.push_back(BAM_CDIFF, 1);
3395 23 Mar 15 peter 84   cig.push_back(BAM_CEQUAL, 51);
3395 23 Mar 15 peter 85   bam.cigar(cig);
3395 23 Mar 15 peter 86   test_alignment_score(suite, bam, reference, 86);
3395 23 Mar 15 peter 87
3395 23 Mar 15 peter 88   cig.clear();
3395 23 Mar 15 peter 89   cig.push_back(BAM_CSOFT_CLIP, 10);
3395 23 Mar 15 peter 90   cig.push_back(BAM_CEQUAL, 38);
3395 23 Mar 15 peter 91   cig.push_back(BAM_CDEL, 2);
3395 23 Mar 15 peter 92   cig.push_back(BAM_CEQUAL, 10);
3395 23 Mar 15 peter 93   cig.push_back(BAM_CINS, 2);
3395 23 Mar 15 peter 94   cig.push_back(BAM_CEQUAL, 10);
3395 23 Mar 15 peter 95   cig.push_back(BAM_CSOFT_CLIP, 2);
3395 23 Mar 15 peter 96   bam.cigar(cig);
3395 23 Mar 15 peter 97   test_alignment_score(suite, bam, reference, 58 - (5+2*2) - (5+2*2));
3395 23 Mar 15 peter 98 }
3395 23 Mar 15 peter 99
3395 23 Mar 15 peter 100
3028 21 Apr 13 peter 101 void test_aux(test::Suite& suite, const BamRead& b)
3028 21 Apr 13 peter 102 {
3028 21 Apr 13 peter 103   suite.out() << "test aux\n";
3028 21 Apr 13 peter 104   BamRead bam(b);
3028 21 Apr 13 peter 105
3028 21 Apr 13 peter 106   char c = 'h';
3028 21 Apr 13 peter 107   int size = bam.aux_size();
3028 21 Apr 13 peter 108   bam.aux_append("ZZ", 'A', 1, (uint8_t*)(&c));
3028 21 Apr 13 peter 109   suite.out() << "size: " << size << " " << bam.aux_size() << "\n";
3028 21 Apr 13 peter 110   suite.add(bam.aux_size() == size+4);
3295 25 Jul 14 peter 111   bam.aux();
3028 21 Apr 13 peter 112   bam.aux("ZZ");
3028 21 Apr 13 peter 113   char c1 = bam_aux2A(bam.aux("ZZ"));
3028 21 Apr 13 peter 114   suite.out() << c << "==" << c1 << "\n";
3028 21 Apr 13 peter 115   suite.add(c==c1);
3028 21 Apr 13 peter 116   bam.aux_del("ZZ");
3028 21 Apr 13 peter 117   if (!suite.add(bam.aux_size() == size))
3028 21 Apr 13 peter 118     suite.err() << "error: incorrect size: " << bam.aux_size() << "\n";
3028 21 Apr 13 peter 119 }
3028 21 Apr 13 peter 120
3028 21 Apr 13 peter 121
2979 31 Jan 13 peter 122 void test_cigar(test::Suite& suite, const BamRead& b, const BamHeader& hdr)
2979 31 Jan 13 peter 123 {
2979 31 Jan 13 peter 124   BamRead bam(b);
2979 31 Jan 13 peter 125
3295 25 Jul 14 peter 126   bam.cigar();
3295 25 Jul 14 peter 127   bam.cigar(0);
3295 25 Jul 14 peter 128   bam.cigar_op(0);
3295 25 Jul 14 peter 129   bam.cigar_oplen(0);
3295 25 Jul 14 peter 130
2979 31 Jan 13 peter 131   suite.out() << "test cigar:\n";
2979 31 Jan 13 peter 132   suite.out() << bam.cigar_str() << "\n";
2979 31 Jan 13 peter 133   OutBamFile os("cigar_test.bam", hdr);
3295 25 Jul 14 peter 134   utility::Aligner::Cigar cig;
2979 31 Jan 13 peter 135
2979 31 Jan 13 peter 136   bam.cigar(cig);
2979 31 Jan 13 peter 137   suite.out() << bam.cigar_str() << "\n";
3359 23 Nov 14 peter 138   if (!suite.add(bam.cigar_str()==""))
3359 23 Nov 14 peter 139     suite.err() << "error: expected cigar to be empty\n";
2979 31 Jan 13 peter 140   os.write(bam);
2979 31 Jan 13 peter 141
3295 25 Jul 14 peter 142   cig.push_back(BAM_CMATCH, bam.sequence_length());
2979 31 Jan 13 peter 143   bam.cigar(cig);
2979 31 Jan 13 peter 144   suite.out() << bam.cigar_str() << "\n";
2979 31 Jan 13 peter 145   suite.add(bam.cigar_str()=="100M");
2979 31 Jan 13 peter 146   os.write(bam);
2979 31 Jan 13 peter 147
3295 25 Jul 14 peter 148   cig.clear();
3295 25 Jul 14 peter 149   cig.push_back(BAM_CMATCH, 50);
3295 25 Jul 14 peter 150   cig.push_back(BAM_CDEL, 2);
3295 25 Jul 14 peter 151   cig.push_back(BAM_CMATCH, 50);
2979 31 Jan 13 peter 152   bam.cigar(cig);
2979 31 Jan 13 peter 153   suite.out() << bam.cigar_str() << "\n";
2979 31 Jan 13 peter 154   suite.add(bam.cigar_str()=="50M2D50M");
2979 31 Jan 13 peter 155   os.write(bam);
2979 31 Jan 13 peter 156
2979 31 Jan 13 peter 157   // check that other elements in data* is preserved
2979 31 Jan 13 peter 158   if (!suite.add(same_query_name(b, bam)))
2979 31 Jan 13 peter 159     suite.err() << "error: \n'" << b.name() << "'\n'" << bam.name() << "'\n";
3384 11 Mar 15 peter 160   if (!suite.add(!less_query_name(b, bam)))
3384 11 Mar 15 peter 161     suite.err() << "error: less_query_name:\n'"
3384 11 Mar 15 peter 162                 << b.name() << "'\n'" << bam.name() << "'\n";
2979 31 Jan 13 peter 163
3384 11 Mar 15 peter 164   BamLessName less_name;
3384 11 Mar 15 peter 165   if (!suite.add(!less_name(b, bam)))
3384 11 Mar 15 peter 166     suite.err() << "error: BamLessName:\n'"
3384 11 Mar 15 peter 167                 << b.name() << "'\n'" << bam.name() << "'\n";
3384 11 Mar 15 peter 168
2979 31 Jan 13 peter 169   if (!suite.add(b.sequence()==bam.sequence()))
2979 31 Jan 13 peter 170     suite.err() << "error: \n'" << b.sequence() << "'\n'"
2979 31 Jan 13 peter 171                 << bam.sequence() << "'\n";
2979 31 Jan 13 peter 172
2979 31 Jan 13 peter 173   for (int32_t i=0; i<b.core().l_qseq; ++i)
2979 31 Jan 13 peter 174     if (bam.qual(i) != b.qual(i)) {
2979 31 Jan 13 peter 175       suite.err() << "error: qual: " << i << " " << bam.qual(i)
2979 31 Jan 13 peter 176                   << " not euqal " << b.qual(i) << "\n";
2979 31 Jan 13 peter 177       suite.add(false);
2979 31 Jan 13 peter 178     }
2979 31 Jan 13 peter 179
2979 31 Jan 13 peter 180   os.close();
2979 31 Jan 13 peter 181 }
2979 31 Jan 13 peter 182
3031 27 Apr 13 peter 183
3306 21 Aug 14 peter 184 void test_end(test::Suite& suite, BamRead bam)
3306 21 Aug 14 peter 185 {
3306 21 Aug 14 peter 186   suite.out() << "test end:\n";
3306 21 Aug 14 peter 187   suite.out() << bam.pos() << "-" << bam.end() << " " << bam.cigar_str()
3306 21 Aug 14 peter 188               << "\n";
3306 21 Aug 14 peter 189   uint32_t end = bam.end();
3306 21 Aug 14 peter 190   if (!suite.add(end==24093)) {
3306 21 Aug 14 peter 191     suite.err() << "error: incorrect end: " << end << "\n";
3306 21 Aug 14 peter 192   }
3306 21 Aug 14 peter 193
3306 21 Aug 14 peter 194   // change the cigar to include a DIFF and EQUAL
3306 21 Aug 14 peter 195   utility::Aligner::Cigar cig;
3306 21 Aug 14 peter 196   cig.push_back(BAM_CEQUAL, 55);
3306 21 Aug 14 peter 197   cig.push_back(BAM_CDIFF, 1);
3306 21 Aug 14 peter 198   cig.push_back(BAM_CEQUAL, 44);
3306 21 Aug 14 peter 199   bam.cigar(cig);
3306 21 Aug 14 peter 200   suite.out() << bam.pos() << "-" << bam.end() << " " << bam.cigar_str()
3306 21 Aug 14 peter 201               << "\n";
3306 21 Aug 14 peter 202   end = bam.end();
3306 21 Aug 14 peter 203   if (!suite.add(end==24093)) {
3306 21 Aug 14 peter 204     suite.err() << "error: incorrect end: " << end << "\n";
3306 21 Aug 14 peter 205   }
3306 21 Aug 14 peter 206 }
3306 21 Aug 14 peter 207
3306 21 Aug 14 peter 208
3031 27 Apr 13 peter 209 void test_sequence(test::Suite& suite, const BamRead& b)
3031 27 Apr 13 peter 210 {
3031 27 Apr 13 peter 211   suite.out() << "test sequence\n";
3031 27 Apr 13 peter 212   BamRead bam(b);
3031 27 Apr 13 peter 213   std::string seq = b.sequence();
3031 27 Apr 13 peter 214   std::vector<uint8_t> qual;
3031 27 Apr 13 peter 215   qual.reserve(seq.size());
3031 27 Apr 13 peter 216   for (size_t i=0; i<seq.size(); ++i)
3031 27 Apr 13 peter 217     qual.push_back(bam.qual(i));
3031 27 Apr 13 peter 218   assert(seq.size()==100);
3031 27 Apr 13 peter 219   bam.sequence(seq, qual);
3031 27 Apr 13 peter 220
3031 27 Apr 13 peter 221   if (bam.sequence()!=seq) {
3031 27 Apr 13 peter 222     suite.add(false);
3031 27 Apr 13 peter 223     suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
3031 27 Apr 13 peter 224     suite.err() << "expected: " << seq << "\n";
3031 27 Apr 13 peter 225   }
3031 27 Apr 13 peter 226
3889 25 Mar 20 peter 227
3031 27 Apr 13 peter 228   // just for consistency
3295 25 Jul 14 peter 229   utility::Aligner::Cigar cig;
3031 27 Apr 13 peter 230   bam.cigar(cig);
3031 27 Apr 13 peter 231   bam.core().flag &= ~BAM_FUNMAP;
3031 27 Apr 13 peter 232
3031 27 Apr 13 peter 233   // trim off one base
3031 27 Apr 13 peter 234   seq.resize(seq.size()-1);
3031 27 Apr 13 peter 235   qual.resize(qual.size()-1);
3031 27 Apr 13 peter 236   bam.sequence(seq, qual);
3031 27 Apr 13 peter 237
3031 27 Apr 13 peter 238   if (bam.sequence()!=seq) {
3031 27 Apr 13 peter 239     suite.add(false);
3031 27 Apr 13 peter 240     suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
3031 27 Apr 13 peter 241     suite.err() << "expected: " << seq << "\n";
3031 27 Apr 13 peter 242   }
3031 27 Apr 13 peter 243
3031 27 Apr 13 peter 244   // extend sequence
3031 27 Apr 13 peter 245   seq.resize(120, 'A');
3031 27 Apr 13 peter 246   qual.resize(120, 0);
3031 27 Apr 13 peter 247   bam.sequence(seq, qual);
3031 27 Apr 13 peter 248   if (bam.sequence()!=seq) {
3031 27 Apr 13 peter 249     suite.add(false);
3031 27 Apr 13 peter 250     suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
3031 27 Apr 13 peter 251     suite.err() << "expected: " << seq << "\n";
3031 27 Apr 13 peter 252   }
3889 25 Mar 20 peter 253
3889 25 Mar 20 peter 254
3889 25 Mar 20 peter 255   // set arbitrary sequence
3889 25 Mar 20 peter 256   seq = "CCCCGGTGCGAAAAAAATATATATATATT";
3889 25 Mar 20 peter 257   qual.resize(seq.size(), 60);
3889 25 Mar 20 peter 258   bam.sequence(seq, qual);
3889 25 Mar 20 peter 259   if (bam.sequence()!=seq) {
3889 25 Mar 20 peter 260     suite.add(false);
3889 25 Mar 20 peter 261     suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
3889 25 Mar 20 peter 262     suite.err() << "expected: " << seq << "\n";
3889 25 Mar 20 peter 263   }
3031 27 Apr 13 peter 264 }
3031 27 Apr 13 peter 265
3031 27 Apr 13 peter 266
3055 16 Jun 13 peter 267 void test_name(test::Suite& suite, const BamRead& b)
3055 16 Jun 13 peter 268 {
3055 16 Jun 13 peter 269   suite.out() << "test name\n";
3055 16 Jun 13 peter 270   BamRead bam(b);
3055 16 Jun 13 peter 271   std::string name = "my-new-name";
3055 16 Jun 13 peter 272   bam.name(name);
3055 16 Jun 13 peter 273   if (bam.name() != name) {
3055 16 Jun 13 peter 274     suite.add(false);
3055 16 Jun 13 peter 275     suite.err() << "error: name: '" << bam.name()
3055 16 Jun 13 peter 276                 << "' expected: '" << name << "'\n";
3055 16 Jun 13 peter 277   }
3055 16 Jun 13 peter 278   // check that other fields are not changed
3055 16 Jun 13 peter 279   if (bam.cigar_str()!=b.cigar_str()) {
3055 16 Jun 13 peter 280     suite.add(false);
3055 16 Jun 13 peter 281     suite.err() << "error: cigar str:" << bam.cigar_str() << " != "
3055 16 Jun 13 peter 282                 << b.cigar_str() << "\n";
3055 16 Jun 13 peter 283   }
3055 16 Jun 13 peter 284   if (bam.sequence()!=b.sequence()) {
3055 16 Jun 13 peter 285     suite.add(false);
3055 16 Jun 13 peter 286     suite.err() << "error: sequence:" << bam.sequence() << " != "
3055 16 Jun 13 peter 287                 << b.sequence() << "\n";
3055 16 Jun 13 peter 288   }
3055 16 Jun 13 peter 289 }
3055 16 Jun 13 peter 290
3055 16 Jun 13 peter 291
3160 13 Jan 14 peter 292 void test_open3(test::Suite& suite, const BamHeader& header)
3160 13 Jan 14 peter 293 {
3160 13 Jan 14 peter 294   OutBamFile out("yaya.bam", header, 0);
3160 13 Jan 14 peter 295   out.close();
3160 13 Jan 14 peter 296   out.open("yaya.bam", header, 9);
3160 13 Jan 14 peter 297   out.close();
3160 13 Jan 14 peter 298   try {
3160 13 Jan 14 peter 299     out.open("yaya.bam", header, 10);
3160 13 Jan 14 peter 300     suite.add(false);
3160 13 Jan 14 peter 301     suite.err() << "open(\"yaya.bam\", header, 10): did not throw\n";
3160 13 Jan 14 peter 302   }
3160 13 Jan 14 peter 303   catch (std::invalid_argument& e) {
3160 13 Jan 14 peter 304     suite.err() << "expected exception: " << e.what() << "\n";
3160 13 Jan 14 peter 305   }
3160 13 Jan 14 peter 306   out.close();
3160 13 Jan 14 peter 307 }
3160 13 Jan 14 peter 308
3160 13 Jan 14 peter 309
3162 15 Jan 14 peter 310 void test_write(test::Suite& suite, const BamHeader& header)
3162 15 Jan 14 peter 311 {
3162 15 Jan 14 peter 312   OutBamFile out;
3162 15 Jan 14 peter 313   BamRead read;
3162 15 Jan 14 peter 314   try {
3162 15 Jan 14 peter 315     // trying to write to a file not open
3162 15 Jan 14 peter 316     out.write(read);
3162 15 Jan 14 peter 317     suite.add(false);
3162 15 Jan 14 peter 318     suite.err() << "write(BamRead): did not throw\n";
3162 15 Jan 14 peter 319   }
3162 15 Jan 14 peter 320   catch (utility::runtime_error& e) {
3162 15 Jan 14 peter 321     suite.err() << "expected exception: " << e.what() << "\n";
3162 15 Jan 14 peter 322   }
3463 03 Feb 16 peter 323   // test ticket #855, passing empty header when opening OutBamFile
3463 03 Feb 16 peter 324   try {
3463 03 Feb 16 peter 325     BamHeader hdr;
3463 03 Feb 16 peter 326     out.open("aha", hdr);
3463 03 Feb 16 peter 327     suite.add(false);
3463 03 Feb 16 peter 328     suite.err() << "OutBamFile.open did not throw\n";
3463 03 Feb 16 peter 329   }
3463 03 Feb 16 peter 330   catch (utility::runtime_error& e) {
3463 03 Feb 16 peter 331     suite.err() << "expected exception: " << e.what() << "\n";
3463 03 Feb 16 peter 332   }
3162 15 Jan 14 peter 333 }
3162 15 Jan 14 peter 334
3162 15 Jan 14 peter 335
3674 01 Aug 17 peter 336 void remove_if(const std::string& fn)
3674 01 Aug 17 peter 337 {
3674 01 Aug 17 peter 338   utility::FileUtil file(fn);
3674 01 Aug 17 peter 339   if (file.exists())
3674 01 Aug 17 peter 340     utility::remove(fn);
3674 01 Aug 17 peter 341 }
3674 01 Aug 17 peter 342
3674 01 Aug 17 peter 343
3674 01 Aug 17 peter 344 void test_index(test::Suite& suite, const std::string& fn)
3674 01 Aug 17 peter 345 {
3674 01 Aug 17 peter 346   try {
3674 01 Aug 17 peter 347     suite.out() << "test_index\n";
3674 01 Aug 17 peter 348     std::string foo = "foo.bam";
3674 01 Aug 17 peter 349     std::string bai = foo + ".bai";
3674 01 Aug 17 peter 350     remove_if(foo);
3674 01 Aug 17 peter 351     remove_if(bai);
3674 01 Aug 17 peter 352
3674 01 Aug 17 peter 353     omic::InBamFile in(fn);
3674 01 Aug 17 peter 354     omic::OutBamFile out(foo, in.header());
3674 01 Aug 17 peter 355     omic::BamRead bam;
3674 01 Aug 17 peter 356     while (in.read(bam))
3674 01 Aug 17 peter 357       out.write(bam);
3674 01 Aug 17 peter 358     in.close();
3674 01 Aug 17 peter 359     out.close();
3674 01 Aug 17 peter 360     suite.out() << "try creating index from OutBamFile '" << foo << "'\n";
3674 01 Aug 17 peter 361     out.build_index();
3674 01 Aug 17 peter 362
3674 01 Aug 17 peter 363     suite.out() << "try reading " << foo << ".bai\n";
3674 01 Aug 17 peter 364
3674 01 Aug 17 peter 365     in.open(foo);
3674 01 Aug 17 peter 366     in.index();
3674 01 Aug 17 peter 367     in.close();
3674 01 Aug 17 peter 368
3674 01 Aug 17 peter 369     utility::remove(bai);
3674 01 Aug 17 peter 370
3674 01 Aug 17 peter 371     suite.out() << "try creating index from InBamFile '" << foo << "'\n";
3674 01 Aug 17 peter 372
3674 01 Aug 17 peter 373     in.open(foo);
3674 01 Aug 17 peter 374     if (!in.have_index())
3674 01 Aug 17 peter 375       in.build_index();
3674 01 Aug 17 peter 376     suite.out() << "try reading " << foo << ".bai\n";
3674 01 Aug 17 peter 377     in.index();
3674 01 Aug 17 peter 378     in.close();
3674 01 Aug 17 peter 379   }
3674 01 Aug 17 peter 380   catch (std::exception& e) {
3674 01 Aug 17 peter 381     suite.add(false);
3674 01 Aug 17 peter 382     suite.err() << "error in test_index: what(): " << e.what() << "\n";
3674 01 Aug 17 peter 383   }
3674 01 Aug 17 peter 384 }
3674 01 Aug 17 peter 385
3674 01 Aug 17 peter 386
2972 26 Jan 13 peter 387 void test1(test::Suite& suite)
2972 26 Jan 13 peter 388 {
2972 26 Jan 13 peter 389   std::string file = "../../data/foo.sorted.bam";
2972 26 Jan 13 peter 390
2972 26 Jan 13 peter 391   InBamFile in;
2972 26 Jan 13 peter 392   in.open(file);
3464 09 Feb 16 peter 393   // test idxstats
3464 09 Feb 16 peter 394   if (in.n_mapped(1)!=795) {
3464 09 Feb 16 peter 395     suite.add(false);
3464 09 Feb 16 peter 396     suite.err() << "error: InBamFile::n_mapped: "
3464 09 Feb 16 peter 397                 << in.n_mapped(1) << " expected: 795\n";
3464 09 Feb 16 peter 398   }
3464 09 Feb 16 peter 399   if (in.n_unmapped(1)!=5) {
3464 09 Feb 16 peter 400     suite.add(false);
3464 09 Feb 16 peter 401     suite.err() << "error: InBamFile::n_unmapped: "
3464 09 Feb 16 peter 402                 << in.n_unmapped(1) << " expected: 5\n";
3464 09 Feb 16 peter 403   }
3464 09 Feb 16 peter 404   if (in.n_no_coordinate()!=0) {
3464 09 Feb 16 peter 405     suite.add(false);
3464 09 Feb 16 peter 406     suite.err() << "error: InBamFile::n_no_coordinate: "
3464 09 Feb 16 peter 407                 << in.n_no_coordinate() << " expected: 0\n";
3464 09 Feb 16 peter 408   }
3698 25 Sep 17 peter 409   // test that idxstats works also when tid is empty (no reads, bug #895)
3698 25 Sep 17 peter 410   try {
3698 25 Sep 17 peter 411     if (in.n_mapped(0)!=0) {
3698 25 Sep 17 peter 412       suite.add(false);
3698 25 Sep 17 peter 413       suite.err() << "error: InBamFile::n_mapped: "
3698 25 Sep 17 peter 414                   << in.n_mapped(0) << " expected: 0\n";
3698 25 Sep 17 peter 415     }
3698 25 Sep 17 peter 416     if (in.n_unmapped(0)!=0) {
3698 25 Sep 17 peter 417       suite.add(false);
3698 25 Sep 17 peter 418       suite.err() << "error: InBamFile::n_unmapped: "
3698 25 Sep 17 peter 419                   << in.n_unmapped(0) << " expected: 0\n";
3698 25 Sep 17 peter 420     }
3698 25 Sep 17 peter 421   }
3698 25 Sep 17 peter 422   catch (utility::runtime_error& e) {
3698 25 Sep 17 peter 423     suite.add(false);
3698 25 Sep 17 peter 424     suite.err() << "error: exception: " << e.what() << "\n";
3698 25 Sep 17 peter 425   }
3464 09 Feb 16 peter 426
2972 26 Jan 13 peter 427   BamRead bam;
2972 26 Jan 13 peter 428   in.read(bam);
3295 25 Jul 14 peter 429   // test that functions are implemented
3306 21 Aug 14 peter 430   test_end(suite, bam);
3295 25 Jul 14 peter 431   bam.flag();
3732 11 Apr 18 peter 432   bam.flag_bit_is_set(48);
3731 11 Apr 18 peter 433   bam.get_flag_bit(32);
3731 11 Apr 18 peter 434   bam.set_flag_bit(32);
3731 11 Apr 18 peter 435   bam.unset_flag_bit(32);
3731 11 Apr 18 peter 436   bam.flip_flag_bit(32);
3733 11 Apr 18 peter 437   bam.strand();
3733 11 Apr 18 peter 438   bam.mstrand();
3295 25 Jul 14 peter 439   bam.mpos();
3295 25 Jul 14 peter 440   bam.mtid();
3295 25 Jul 14 peter 441   bam.pos();
3295 25 Jul 14 peter 442   bam.tid();
3295 25 Jul 14 peter 443   bam.swap(bam);
2972 26 Jan 13 peter 444   std::string str = "AGCTTTGTTATTATTTGAGGGTGTGGTTCAGTTGTAAACAGTGTATG"
2972 26 Jan 13 peter 445     "TTTTAGAATTTGTGTTATTGTGATGGCGATGACCAAGTACAACATATTTCCCA";
2972 26 Jan 13 peter 446   std::string seq = bam.sequence();
2972 26 Jan 13 peter 447   if (str!=seq) {
2972 26 Jan 13 peter 448     suite.err() << "incorrect sequence: '" << bam.sequence() << "'\n";
2972 26 Jan 13 peter 449     suite.err() << "expected:           '" << str << "'\n";
2972 26 Jan 13 peter 450     suite.add(false);
2972 26 Jan 13 peter 451   }
2987 18 Feb 13 peter 452   uint8_t c = bam.sequence(24);
2987 18 Feb 13 peter 453   if (c != 4) {
2987 18 Feb 13 peter 454     suite.add(false);
2987 18 Feb 13 peter 455     suite.err() << "error: sequence: " << static_cast<int>(c)
2987 18 Feb 13 peter 456                 << " expected " << 4 << "\n";
2987 18 Feb 13 peter 457   }
2987 18 Feb 13 peter 458   bam.sequence(24, 8);
2987 18 Feb 13 peter 459   c = bam.sequence(24);
2987 18 Feb 13 peter 460   if (c != 8) {
2987 18 Feb 13 peter 461     suite.add(false);
2987 18 Feb 13 peter 462     suite.err() << "error: sequence: " << static_cast<int>(c)
2987 18 Feb 13 peter 463                 << " expected " << 8 << "\n";
2987 18 Feb 13 peter 464   }
3026 06 Apr 13 peter 465   uint8_t q = bam.qual(24);
3026 06 Apr 13 peter 466   if (q != 72) {
3026 06 Apr 13 peter 467     suite.add(false);
3026 06 Apr 13 peter 468     suite.err() << "error: quality: " << static_cast<int>(q)
3026 06 Apr 13 peter 469                 << " expected " << 72 << "\n";
3026 06 Apr 13 peter 470   }
3026 06 Apr 13 peter 471   uint8_t new_q = 60;
3026 06 Apr 13 peter 472   bam.qual(24, new_q);
3026 06 Apr 13 peter 473   q = bam.qual(24);
3026 06 Apr 13 peter 474   if (q != new_q) {
3026 06 Apr 13 peter 475     suite.add(false);
3026 06 Apr 13 peter 476     suite.err() << "error: quality: " << static_cast<int>(q)
3026 06 Apr 13 peter 477                 << " expected " << new_q << "\n";
3026 06 Apr 13 peter 478   }
3395 23 Mar 15 peter 479   test_alignment_score(suite, bam);
2979 31 Jan 13 peter 480   test_cigar(suite, bam, in.header());
3028 21 Apr 13 peter 481   test_aux(suite, bam);
3031 27 Apr 13 peter 482   test_sequence(suite, bam);
3055 16 Jun 13 peter 483   test_name(suite, bam);
3160 13 Jan 14 peter 484   test_open3(suite, in.header());
3162 15 Jan 14 peter 485   test_write(suite, in.header());
3674 01 Aug 17 peter 486   test_index(suite, file);
2972 26 Jan 13 peter 487 }
2972 26 Jan 13 peter 488 #endif

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