source: trunk/test/pileup.cc @ 3369

// $Id: pileup.cc 3369 2015-02-10 07:15:55Z peter $

/*
  Copyright (C) 2014 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library 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.

  The yat library 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Suite.h"

#ifdef HAVE_LIBBAM
#include <yat/omic/BamFile.h>
#include <yat/omic/BamReadIterator.h>
#include "yat/omic/Pileup.h"
#endif

#include <yat/utility/Aligner.h>

#include <algorithm>
#include <cassert>
#include <string>
#include <vector>

using namespace theplu::yat;

// The simplest of genotypers; reporting the first and second "thirdiles"
std::string genotype(std::string str)
{
  std::string res("  ");
  std::sort(str.begin(), str.end());
  res[0] = str[str.size()/3];
  res[1] = str[2*str.size()/3];
  return res;
}


#ifdef YAT_HAVE_LIBBAM
char seq_nt16(size_t x)
{
#if YAT_HAVE_HTSLIB
  return seq_nt16_str[x];
#elif HAVE_BAM_NT16_REV_TABLE
  return bam_nt16_rev_table[x];
#else
  std::string table = "=ACMGRSVTWYHKDBN";
  return table[c];
#endif
}

void test1(test::Suite& suite, const std::string& fn)
{
  omic::InBamFile in;
  in.open(fn);

  using omic::BamReadIterator;
  BamReadIterator first(in);
  BamReadIterator last;

  size_t count = 0;
  std::vector<std::string> correct_gt;
  correct_gt.push_back("AA");
  correct_gt.push_back("CC");
  correct_gt.push_back("CC");
  correct_gt.push_back("TT");
  correct_gt.push_back("=A");
  correct_gt.push_back("AA");
  correct_gt.push_back("AA");
  correct_gt.push_back("GG");
  correct_gt.push_back("AA");
  correct_gt.push_back("AA");
  correct_gt.push_back("AA");

  correct_gt.push_back("GG");
  correct_gt.push_back("GG");
  correct_gt.push_back("CC");
  correct_gt.push_back("CC");
  correct_gt.push_back("CC");

  correct_gt.push_back("TT");
  correct_gt.push_back("TT");
  correct_gt.push_back("=G");
  correct_gt.push_back("=C");

  using omic::Pileup;
  Pileup<BamReadIterator> pileup(first, last);
  while (pileup.good()) {
    if (pileup.pos()+1 < 24341) {
      pileup.increment();
      continue;
    }
    Pileup<BamReadIterator>::const_iterator iter = pileup.begin();
    std::string str;
    for (; iter!=pileup.end(); ++iter) {
      if (iter->bam().end() <= pileup.pos()) {
        suite.err() << "error: " << pileup.pos() << " "
                    << iter->bam().pos() << " "
                    << iter->bam().end() << " "
                    << iter->bam().cigar_str()
                    << "\n";
        suite.add(false);
        continue;
      }
      str += seq_nt16(iter->sequence());
    }

    if (pileup.pos()+1 == 24341 && count!=0) {
      suite.err() << "incorrect count\n";
      suite.add(false);
    }

    if (count < correct_gt.size())
      if (genotype(str) != correct_gt[count]) {
        suite.err() << "incorrect genotype: " << genotype(str) << "\n";
        suite.err() << "expected: " << correct_gt[count] << "\n";
        suite.add(false);
      }

    if (count <= 3) {
      if (pileup.pos()+1 != static_cast<int>(count+24341)) {
        suite.err() << "incorrect pos or count\n";
        suite.err() << "pos: " << pileup.pos()+1 << "\n";
        suite.err() << "count: " << count << "\n";
        suite.add(false);
      }
    }
    else if (count == 4 && !pileup.skip_ref() ) {
      suite.out() << "expected skip ref\n";
      suite.add(false);
    }
    else if (count==5) {
      if (pileup.pos()+1 != static_cast<int>(count+24340)) {
        suite.err() << "incorrect pos or count\n";
        suite.err() << "pos: " << pileup.pos()+1 << "\n";
        suite.err() << "count: " << count << "\n";
        suite.add(false);
      }
    }
    pileup.increment();
    ++count;
  }
}


void test2(test::Suite& suite, const std::string& fn)
{
  std::ostringstream error;
  theplu::yat::omic::InBamFile in(fn);
  using omic::BamRead;
  using omic::BamReadIterator;
  using omic::Pileup;
  BamReadIterator first(in, 1, 24150, 25000);
  BamReadIterator last;

  std::vector<BamRead> reads(2);
  reads[0] = *first;
  ++first;
  reads[1] = *first;
  in.close();
  std::string ref = reads[0].sequence();

  suite.out() << "no modifications\n";
  // no modification
  typedef std::vector<BamRead>::iterator iterator;
  typedef Pileup<iterator>::const_iterator plp_iterator;
  for (Pileup<iterator> plp(reads.begin(), reads.end());plp.pos()<24070;
       plp.increment()) {
    for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
      if (i->bam().pos()!=reads[1].pos())
        continue;
      char nt = seq_nt16(i->sequence());
      if (nt != ref[plp.pos()-reads[0].pos()]) {
        error << "'" << nt << "' not equal '"
              << ref[plp.pos()-reads[0].pos()]
              << "' at pos " << reads[0].pos() << " for bam with start pos: "
              << i->bam().pos();
        throw std::runtime_error(error.str());
      }
    }
  }

  // insertion
  suite.out() << "insertion\n";
  assert(reads[1].sequence_length()==100);
  utility::Aligner::Cigar cig;
  cig.push_back(BAM_CMATCH, 10);
  cig.push_back(BAM_CINS, 1);
  cig.push_back(BAM_CMATCH, 89);
  reads[1].cigar(cig);
  int count = 0;
  for (Pileup<iterator> plp(reads.begin(), reads.end());plp.good();
       plp.increment()) {
    for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
      if (!same_query_name(i->bam(), reads[1]))
        continue;
      ++count;
      suite.err() << count << "\n";
      char nt = seq_nt16(i->sequence());
      if (nt != reads[1].sequence()[count-1])
        error << "nt: " << nt << " not as expected. count: " << count << "\n";
      // first 10 matches
      if (count <= 10) {
        if (plp.skip_ref())
          error << count << " unexpected skip_ref\n";
        if (plp.pos() != reads[1].pos()+count-1)
          error << "count: " << count << " with pos: " << plp.pos() << "\n";
      }
      else if (count==11) { // insertion
        if (!plp.skip_ref())
          error << count << " unexpected skip_ref\n";
      }
      else {
        if (plp.skip_ref())
          error << count << " unexpected skip_ref\n";
        if (plp.pos() != reads[1].pos()+count-2)
          error << "count: " << count << " with pos: " << plp.pos() << "\n";
      }

      if (error.str().size())
        throw std::runtime_error(error.str());
    }
  }


  // deletion
  suite.out() << "deletion\n";
  assert(reads[1].sequence_length()==100);
  cig.clear();
  cig.push_back(BAM_CMATCH, 10);
  cig.push_back(BAM_CDEL, 1);
  cig.push_back(BAM_CMATCH, 90);
  assert(cig.size()==3);
  reads[1].cigar(cig);
  count = 0;
  for (Pileup<iterator> plp(reads.begin(), reads.end()); plp.good();
       plp.increment()) {
    for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
      if (!same_query_name(i->bam(), reads[1]))
        continue;
      ++count;
      char nt = seq_nt16(i->sequence());
      if (plp.pos() != reads[1].pos()+count-1)
        error << "count: " << count << " with pos: " << plp.pos() << "\n";
      // first 10 matches
      if (count <= 10) {
        if (plp.skip_ref())
          error << count << " unexpected skip_ref\n";
        if (nt != reads[1].sequence()[count-1])
          error << "nt: " << nt << " not as expected\n";
      }
      else if (count==11) {
        if (plp.skip_ref())
          error << count << " unexpected skip_ref\n";
        if (!i->is_del())
          error << count << " expected deletion\n";
      }
      else {
        if (plp.skip_ref())
          error << count << " unexpected skip_ref\n";
        if (nt != reads[1].sequence()[count-2])
          error << "nt: " << nt << " not as expected\n";
      }

      if (error.str().size())
        throw std::runtime_error(error.str());
    }
  }

  // soft clipped
  suite.out() << "soft clipped\n";
  assert(reads[1].sequence_length()==100);
  cig.clear();
  cig.push_back(BAM_CSOFT_CLIP, 10);
  cig.push_back(BAM_CMATCH, 90);
  reads[1].cigar(cig);
  count = 0;
  for (Pileup<iterator> plp(reads.begin(), reads.end()); plp.good();
       plp.increment()) {
    for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
      if (!same_query_name(i->bam(), reads[1]))
        continue;
      ++count;
      char nt = seq_nt16(i->sequence());
      if (plp.skip_ref())
        error << count << " unexpected skip_ref\n";
      if (plp.pos() != reads[1].pos()+count-1)
        error << "count: " << count << " with pos: " << plp.pos() << "\n";
      if (nt != reads[1].sequence()[count+9])
        error << "count: " << count << " nt: " << nt << " not as expected\n";

      if (error.str().size())
        throw std::runtime_error(error.str());
    }
  }
}
#else
// tests if libbam is not available, should never be run
void test1(test::Suite& suite, const std::string& fn) {assert(0);}
void test2(test::Suite& suite, const std::string& fn) {assert(0);}
#endif

int main(int argc, char* argv[])
{
  test::Suite suite(argc, argv, true);
  std::string fn = "../../data/foo.sorted.bam";
  try {
    test1(suite, fn);
    test2(suite, fn);
  }
  catch (std::runtime_error& e) {
    suite.err() << "test failed: " << e.what() << "\n";
    suite.add(false);
  }
  return suite.return_value();
}