X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=BamWriter.h;h=e014f8d5915445622841b8ca6f680cdf6b9e413d;hb=237bbdf363c9e42ee24e2fd63106dccf20d9bf2f;hp=ff11397e3539d8662f31e85d6a9e8d09309e51f8;hpb=fc69cf6af24c0550e55447fc82f01cb6f90c1c42;p=rsem.git diff --git a/BamWriter.h b/BamWriter.h index ff11397..e014f8d 100644 --- a/BamWriter.h +++ b/BamWriter.h @@ -47,23 +47,13 @@ BamWriter::BamWriter(char inpType, const char* inpF, const char* fn_list, const } assert(in != 0); + //build mappings from external sid to internal sid + transcripts.buildMappings(in->header->n_targets, in->header->target_name); + //generate output's header bam_header_t *out_header = bam_header_dwt(in->header); - - if (out_header->n_targets != transcripts.getM()) { - fprintf(stderr, "Number of reference sequences recorded in the header is not correct! The header contains %d sequences while there should be %d sequences\n", out_header->n_targets, transcripts.getM()); - exit(-1); - } - for (int i = 0; i < out_header->n_targets; i++) { - const Transcript& transcript = transcripts.getTranscriptAt(i + 1); - if (out_header->target_name[i] != transcript.getTranscriptID()) { - fprintf(stderr, "Reference sequence %d's name recorded in the header is not correct! \n", i); - fprintf(stderr, "Name in the header: %s\n", out_header->target_name[i]); - fprintf(stderr, "Should be: %s\n", transcript.getTranscriptID().c_str()); - exit(-1); - } - out_header->target_len[i] = transcript.getLength(); // transcript length without poly(A) tail + out_header->target_len[i] = transcripts.getTranscriptViaEid(i + 1).getLength(); // transcript length without poly(A) tail } std::ostringstream strout; @@ -99,7 +89,7 @@ void BamWriter::work(HitWrapper wrapper) { hit = wrapper.getNextHit(); assert(hit != NULL); - assert(b->core.tid + 1 == hit->getSid()); + assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid()); convert(b, hit->getConPrb()); if (b->core.qual > 0) samwrite(out, b); // output only when MAPQ > 0 } @@ -135,8 +125,8 @@ void BamWriter::work(HitWrapper wrapper) { hit = wrapper.getNextHit(); assert(hit != NULL); - assert(b->core.tid + 1 == hit->getSid()); - assert(b2->core.tid + 1 == hit->getSid()); + assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid()); + assert(transcripts.getInternalSid(b2->core.tid + 1) == hit->getSid()); convert(b, hit->getConPrb()); convert(b2, hit->getConPrb()); @@ -159,8 +149,7 @@ void BamWriter::work(HitWrapper wrapper) { } void BamWriter::convert(bam1_t *b, double prb) { - int sid = b->core.tid + 1; - const Transcript& transcript = transcripts.getTranscriptAt(sid); + const Transcript& transcript = transcripts.getTranscriptViaEid(b->core.tid + 1); int pos = b->core.pos; int readlen = b->core.l_qseq;