#include "sam_rsem_cvt.h"
#include "utils.h"
+#include "my_assert.h"
#include "bc_aux.h"
#include "Transcript.h"
#include "Transcripts.h"
BamConverter::BamConverter(const char* inpF, const char* outF, const char* chr_list, Transcripts& transcripts)
: transcripts(transcripts)
{
- if (transcripts.getType() != 0)
- exitWithError("Genome information is not provided! RSEM cannot convert the transcript bam file!");
+ general_assert(transcripts.getType() == 0, "Genome information is not provided! RSEM cannot convert the transcript bam file!");
in = samopen(inpF, "rb", NULL);
assert(in != 0);
+ transcripts.buildMappings(in->header->n_targets, in->header->target_name);
+
bam_header_t *out_header = sam_header_read2(chr_list);
refmap.clear();
for (int i = 0; i < out_header->n_targets; i++) {
if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
// at least one segment is not properly mapped
- if ((b->core.flag & 0x0004) || isPaired && (b2->core.flag & 0x0004)) continue;
+ if ((b->core.flag & 0x0004) || (isPaired && (b2->core.flag & 0x0004))) continue;
- const Transcript& transcript = transcripts.getTranscriptAt(b->core.tid + 1);
+ const Transcript& transcript = transcripts.getTranscriptViaEid(b->core.tid + 1);
convert(b, transcript);
if (isPaired) {
int pos = b->core.pos;
int readlen = b->core.l_qseq;
- if (readlen == 0) exitWithError("One alignment line has SEQ field as *. RSEM does not support this currently!");
+ general_assert(readlen > 0, "One alignment line has SEQ field as *. RSEM does not support this currently!");
iter = refmap.find(transcript.getSeqName());
assert(iter != refmap.end());