source: trunk/test/bam.cc @ 2986

Last change on this file since 2986 was 2986, checked in by Peter, 10 years ago

merge 0.10-stable branch into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.0 KB
Line 
1// $Id: bam.cc 2986 2013-02-18 08:07:44Z peter $
2//
3// Copyright (C) 2013 Peter Johansson
4//
5// This program is free software; you can redistribute it and/or modify
6// it under the terms of the GNU General Public License as published by
7// the Free Software Foundation; either version 3 of the License, or
8// (at your option) any later version.
9//
10// This program is distributed in the hope that it will be useful, but
11// WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13// General Public License for more details.
14//
15// You should have received a copy of the GNU General Public License
16// along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18#include <config.h>
19
20#include "Suite.h"
21
22#ifdef HAVE_LIBBAM
23#include "yat/omic/BamFile.h"
24#include "yat/omic/BamRead.h"
25#endif
26
27#include <cassert>
28#include <string>
29
30using namespace theplu::yat;
31
32void test1(test::Suite& suite);
33
34int main(int argc, char* argv[])
35{
36  test::Suite suite(argc, argv);
37#ifndef HAVE_LIBBAM
38  suite.out() << "no libbam\n";
39  return EXIT_SKIP;
40#endif
41#ifndef HAVE_SAMTOOLS
42  suite.out() << "no samtools\n";
43  return EXIT_SKIP;
44#endif
45#ifdef HAVE_LIBBAM
46  test1(suite);
47#endif
48  return suite.return_value();
49}
50
51#ifdef HAVE_LIBBAM
52using namespace omic;
53
54void test_cigar(test::Suite& suite, const BamRead& b, const BamHeader& hdr)
55{
56  BamRead bam(b);
57
58  suite.out() << "test cigar:\n";
59  suite.out() << bam.cigar_str() << "\n";
60  OutBamFile os("cigar_test.bam", hdr);
61  std::vector<uint32_t> cig;
62
63  bam.cigar(cig);
64  suite.out() << bam.cigar_str() << "\n";
65  suite.add(bam.cigar_str()=="");
66  os.write(bam);
67
68  cig.resize(1);
69  cig[0] = bam_cigar_gen(bam.sequence_length(), BAM_CMATCH);
70  bam.cigar(cig);
71  suite.out() << bam.cigar_str() << "\n";
72  suite.add(bam.cigar_str()=="100M");
73  os.write(bam);
74
75  cig.resize(3);
76  cig[0] = bam_cigar_gen(50, BAM_CMATCH);
77  cig[1] = bam_cigar_gen(2, BAM_CDEL);
78  cig[2] = bam_cigar_gen(50, BAM_CMATCH);
79  bam.cigar(cig);
80  suite.out() << bam.cigar_str() << "\n";
81  suite.add(bam.cigar_str()=="50M2D50M");
82  os.write(bam);
83
84  // check that other elements in data* is preserved
85  if (!suite.add(same_query_name(b, bam)))
86    suite.err() << "error: \n'" << b.name() << "'\n'" << bam.name() << "'\n";
87
88  if (!suite.add(b.sequence()==bam.sequence()))
89    suite.err() << "error: \n'" << b.sequence() << "'\n'"
90                << bam.sequence() << "'\n";
91
92  for (int32_t i=0; i<b.core().l_qseq; ++i)
93    if (bam.qual(i) != b.qual(i)) {
94      suite.err() << "error: qual: " << i << " " << bam.qual(i)
95                  << " not euqal " << b.qual(i) << "\n";
96      suite.add(false);
97    }
98
99  os.close();
100}
101
102void test1(test::Suite& suite)
103{
104  std::string file = "../../data/foo.sorted.bam";
105
106  InBamFile in;
107  in.open(file);
108  BamRead bam;
109  in.read(bam);
110  std::string str = "AGCTTTGTTATTATTTGAGGGTGTGGTTCAGTTGTAAACAGTGTATG"
111    "TTTTAGAATTTGTGTTATTGTGATGGCGATGACCAAGTACAACATATTTCCCA";
112  std::string seq = bam.sequence();
113  if (str!=seq) {
114    suite.err() << "incorrect sequence: '" << bam.sequence() << "'\n";
115    suite.err() << "expected:           '" << str << "'\n";
116    suite.add(false);
117  }
118  test_cigar(suite, bam, in.header());
119}
120#endif
Note: See TracBrowser for help on using the repository browser.