File: BamWriter.cpp

package info (click to toggle)
libseqlib 1.1.2+dfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 1,508 kB
  • sloc: cpp: 7,176; sh: 805; makefile: 60
file content (148 lines) | stat: -rw-r--r-- 3,295 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#include "SeqLib/BamWalker.h"
#include "SeqLib/BamWriter.h"

#include <stdexcept>

//#define DEBUG_WALKER 1

namespace SeqLib {

  void BamWriter::SetHeader(const SeqLib::BamHeader& h) {
    hdr = h;
  }

  bool BamWriter::WriteHeader() const {
    
    if (hdr.isEmpty()) {
      std::cerr << "BamWriter::WriteHeader - No header supplied. Provide with SetWriteHeader" << std::endl;
      return false;
    }

    if (!fop) {
      std::cerr << "BamWriter::WriteHeader - Output not open for writing. Open with Open()" << std::endl;
      return false;
    }
    
    if (sam_hdr_write(fop.get(), hdr.get()) < 0) {
      std::cerr << "Cannot write header. sam_hdr_write exited with < 0" << std::endl;
      return false;
    }

    return true;
    
  }
  
  bool BamWriter::Close() {

    if (!fop)
      return false;

    fop.reset(); //tr1 compatible
    //fop = NULL; // this clears shared_ptr, calls sam_close (c++11)

    return true;
  }

bool BamWriter::BuildIndex() const {
  
  // throw an error if BAM is not already closed
  if (fop) {
    std::cerr << "Trying to index open BAM. Close first with Close()" << std::endl;
    return false;
  }

  if (m_out.empty()) {
    std::cerr << "Trying to make index, but no BAM specified" << std::endl;
    return false;
  }
  
  // call to htslib to build bai index
  if (sam_index_build(m_out.c_str(), 0) < 0) { // 0 is "min_shift", which is 0 for bai index
    std::cerr << "Failed to create index";
    return false;
  }

  return true;

}

  bool BamWriter::Open(const std::string& f) {

    // don't reopen
    if (fop)
      return false;

    m_out = f;

    // hts open the writer
    fop = SeqPointer<htsFile>(hts_open(m_out.c_str(), output_format.c_str()), htsFile_delete());

    if (!fop) {
      return false;
      //throw std::runtime_error("BamWriter::Open - Cannot open output file: " + f);
    }

    return true;
  }

  BamWriter::BamWriter(int o) {

    switch(o) {
    case BAM :  output_format = "wb"; break;
    case CRAM : output_format = "wc"; break;
    case SAM :  output_format = "w"; break;
    default : throw std::invalid_argument("Invalid writer type");
    }

  }
  

bool BamWriter::WriteRecord(const BamRecord &r)
{
  if (!fop) {
    return false;
  } else {
    if (sam_write1(fop.get(), hdr.get(), r.raw()) < 0)
      return false;
  }

  return true;
}

std::ostream& operator<<(std::ostream& out, const BamWriter& b)
{
  if (b.fop)
    out << "Write format: " << b.fop->format.format;
  out << " Write file " << b.m_out; 
  return out;
}

  //does not return false if file not found
bool BamWriter::SetCramReference(const std::string& ref) {

  if (!fop)
    return false;

  // need to open reference for CRAM writing 
  char* fn_list = samfaipath(ref.c_str()); // eg ref = my.fa  returns my.fa.fai
  if (fn_list) {

    // in theory hts_set_fai_filename should give back < 0
    // if fn_list not there, but it doesnt
    if (!read_access_test(std::string(fn_list)))
      return false;
	
    int status = hts_set_fai_filename(fop.get(), fn_list); 
    if (status < 0) {
      fprintf(stderr, "Failed to use reference \"%s\".\n", fn_list);
      return false;
    }
  } else {
    std::cerr << "Failed to get the reference for CRAM compression" << std::endl;
    return false;
  }

  return true;
}

}