source: trunk/test/bam.cc @ 3295

// $Id: bam.cc 3295 2014-07-25 04:24:39Z peter $
//
// Copyright (C) 2013, 2014 Peter Johansson
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#include <config.h>

#include "Suite.h"

#ifdef YAT_HAVE_LIBBAM
#include "yat/omic/BamFile.h"
#include "yat/omic/BamRead.h"
#endif

#include "yat/utility/Aligner.h"

#include <cassert>
#include <string>

using namespace theplu::yat;

void test1(test::Suite& suite);

int main(int argc, char* argv[])
{
  test::Suite suite(argc, argv, true);
#ifdef YAT_HAVE_LIBBAM
  test1(suite);
#endif
  return suite.return_value();
}

#ifdef YAT_HAVE_LIBBAM
using namespace omic;

void test_aux(test::Suite& suite, const BamRead& b)
{
  suite.out() << "test aux\n";
  BamRead bam(b);

  char c = 'h';
  int size = bam.aux_size();
  bam.aux_append("ZZ", 'A', 1, (uint8_t*)(&c));
  suite.out() << "size: " << size << " " << bam.aux_size() << "\n";
  suite.add(bam.aux_size() == size+4);
  bam.aux();
  bam.aux("ZZ");
  char c1 = bam_aux2A(bam.aux("ZZ"));
  suite.out() << c << "==" << c1 << "\n";
  suite.add(c==c1);
  bam.aux_del("ZZ");
  if (!suite.add(bam.aux_size() == size))
    suite.err() << "error: incorrect size: " << bam.aux_size() << "\n";
}


void test_cigar(test::Suite& suite, const BamRead& b, const BamHeader& hdr)
{
  BamRead bam(b);

  bam.cigar();
  bam.cigar(0);
  bam.cigar_op(0);
  bam.cigar_oplen(0);

  suite.out() << "test cigar:\n";
  suite.out() << bam.cigar_str() << "\n";
  OutBamFile os("cigar_test.bam", hdr);
  utility::Aligner::Cigar cig;

  bam.cigar(cig);
  suite.out() << bam.cigar_str() << "\n";
  suite.add(bam.cigar_str()=="");
  os.write(bam);

  cig.push_back(BAM_CMATCH, bam.sequence_length());
  bam.cigar(cig);
  suite.out() << bam.cigar_str() << "\n";
  suite.add(bam.cigar_str()=="100M");
  os.write(bam);

  cig.clear();
  cig.push_back(BAM_CMATCH, 50);
  cig.push_back(BAM_CDEL, 2);
  cig.push_back(BAM_CMATCH, 50);
  bam.cigar(cig);
  suite.out() << bam.cigar_str() << "\n";
  suite.add(bam.cigar_str()=="50M2D50M");
  os.write(bam);

  // check that other elements in data* is preserved
  if (!suite.add(same_query_name(b, bam)))
    suite.err() << "error: \n'" << b.name() << "'\n'" << bam.name() << "'\n";

  if (!suite.add(b.sequence()==bam.sequence()))
    suite.err() << "error: \n'" << b.sequence() << "'\n'"
                << bam.sequence() << "'\n";

  for (int32_t i=0; i<b.core().l_qseq; ++i)
    if (bam.qual(i) != b.qual(i)) {
      suite.err() << "error: qual: " << i << " " << bam.qual(i)
                  << " not euqal " << b.qual(i) << "\n";
      suite.add(false);
    }

  os.close();
}


void test_sequence(test::Suite& suite, const BamRead& b)
{
  suite.out() << "test sequence\n";
  BamRead bam(b);
  std::string seq = b.sequence();
  std::vector<uint8_t> qual;
  qual.reserve(seq.size());
  for (size_t i=0; i<seq.size(); ++i)
    qual.push_back(bam.qual(i));
  assert(seq.size()==100);
  bam.sequence(seq, qual);

  if (bam.sequence()!=seq) {
    suite.add(false);
    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
    suite.err() << "expected: " << seq << "\n";
  }

  // just for consistency
  utility::Aligner::Cigar cig;
  bam.cigar(cig);
  bam.core().flag &= ~BAM_FUNMAP;

  // trim off one base
  seq.resize(seq.size()-1);
  qual.resize(qual.size()-1);
  bam.sequence(seq, qual);

  if (bam.sequence()!=seq) {
    suite.add(false);
    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
    suite.err() << "expected: " << seq << "\n";
  }

  // extend sequence
  seq.resize(120, 'A');
  qual.resize(120, 0);
  bam.sequence(seq, qual);
  if (bam.sequence()!=seq) {
    suite.add(false);
    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
    suite.err() << "expected: " << seq << "\n";
  }
}


void test_name(test::Suite& suite, const BamRead& b)
{
  suite.out() << "test name\n";
  BamRead bam(b);
  std::string name = "my-new-name";
  bam.name(name);
  if (bam.name() != name) {
    suite.add(false);
    suite.err() << "error: name: '" << bam.name()
                << "' expected: '" << name << "'\n";
  }
  // check that other fields are not changed
  if (bam.cigar_str()!=b.cigar_str()) {
    suite.add(false);
    suite.err() << "error: cigar str:" << bam.cigar_str() << " != "
                << b.cigar_str() << "\n";
  }
  if (bam.sequence()!=b.sequence()) {
    suite.add(false);
    suite.err() << "error: sequence:" << bam.sequence() << " != "
                << b.sequence() << "\n";
  }
}


void test_open3(test::Suite& suite, const BamHeader& header)
{
  OutBamFile out("yaya.bam", header, 0);
  out.close();
  out.open("yaya.bam", header, 9);
  out.close();
  try {
    out.open("yaya.bam", header, 10);
    suite.add(false);
    suite.err() << "open(\"yaya.bam\", header, 10): did not throw\n";
  }
  catch (std::invalid_argument& e) {
    suite.err() << "expected exception: " << e.what() << "\n";
  }
  out.close();
}


void test_write(test::Suite& suite, const BamHeader& header)
{
  OutBamFile out;
  BamRead read;
  try {
    // trying to write to a file not open
    out.write(read);
    suite.add(false);
    suite.err() << "write(BamRead): did not throw\n";
  }
  catch (utility::runtime_error& e) {
    suite.err() << "expected exception: " << e.what() << "\n";
  }
}


void test1(test::Suite& suite)
{
  std::string file = "../../data/foo.sorted.bam";

  InBamFile in;
  in.open(file);
  BamRead bam;
  in.read(bam);
  // test that functions are implemented
  bam.end();
  bam.flag();
  bam.mpos();
  bam.mtid();
  bam.pos();
  bam.tid();
  bam.swap(bam);
  std::string str = "AGCTTTGTTATTATTTGAGGGTGTGGTTCAGTTGTAAACAGTGTATG"
    "TTTTAGAATTTGTGTTATTGTGATGGCGATGACCAAGTACAACATATTTCCCA";
  std::string seq = bam.sequence();
  if (str!=seq) {
    suite.err() << "incorrect sequence: '" << bam.sequence() << "'\n";
    suite.err() << "expected:           '" << str << "'\n";
    suite.add(false);
  }
  uint8_t c = bam.sequence(24);
  if (c != 4) {
    suite.add(false);
    suite.err() << "error: sequence: " << static_cast<int>(c)
                << " expected " << 4 << "\n";
  }
  bam.sequence(24, 8);
  c = bam.sequence(24);
  if (c != 8) {
    suite.add(false);
    suite.err() << "error: sequence: " << static_cast<int>(c)
                << " expected " << 8 << "\n";
  }
  uint8_t q = bam.qual(24);
  if (q != 72) {
    suite.add(false);
    suite.err() << "error: quality: " << static_cast<int>(q)
                << " expected " << 72 << "\n";
  }
  uint8_t new_q = 60;
  bam.qual(24, new_q);
  q = bam.qual(24);
  if (q != new_q) {
    suite.add(false);
    suite.err() << "error: quality: " << static_cast<int>(q)
                << " expected " << new_q << "\n";
  }
  test_cigar(suite, bam, in.header());
  test_aux(suite, bam);
  test_sequence(suite, bam);
  test_name(suite, bam);
  test_open3(suite, in.header());
  test_write(suite, in.header());
}
#endif