From: Bo Li Date: Fri, 11 Nov 2011 00:16:05 +0000 (-0600) Subject: add optional field XS for cufflinks X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=commitdiff_plain;h=109200af109807720ab01febf9b7620c06d4378b add optional field XS for cufflinks --- diff --git a/BamWriter.h b/BamWriter.h index e29e36b..558a44a 100644 --- a/BamWriter.h +++ b/BamWriter.h @@ -240,29 +240,30 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) { push_qual(qual, readlen, transcript.getStrand(), data); free(b->data); - b->m_data = b->data_len = data.size() + 7; // 7 extra bytes for ZW tag - b->l_aux = 7; + b->m_data = b->data_len = data.size(); + b->l_aux = 0; b->data = (uint8_t*)malloc(b->m_data); for (int i = 0; i < b->data_len; i++) b->data[i] = data[i]; b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b))); + + char strand = transcript.getStrand(); + bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand); + } else { - b->m_data = b->data_len = b->data_len - b->l_aux + 7; // 7 extra bytes for ZW tag - b->l_aux = 7; + b->m_data = b->data_len = b->data_len - b->l_aux; + b->l_aux = 0; b->data = (uint8_t*)realloc(b->data, b->m_data); } - - + if (cqname != bam1_qname(b)) { if (!hmap.empty()) { for (hmapIter = hmap.begin(); hmapIter != hmap.end(); hmapIter++) { bam1_t *tmp_b = hmapIter->first.b; tmp_b->core.qual = getMAPQ(hmapIter->second); - uint8_t *p = bam1_aux(tmp_b); - *p = 'Z'; ++p; *p = 'W'; ++p; *p = 'f'; ++p; float val = (float)hmapIter->second; - memcpy(p, &val, 4); + bam_aux_append(tmp_b, "ZW", 'f', 4, (uint8_t*)&val); if (tmp_b->core.qual > 0) samwrite(out, tmp_b); // output only when MAPQ > 0 bam_destroy1(tmp_b); // now hmapIter->b makes no sense } @@ -286,10 +287,8 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) { for (hmapIter = hmap.begin(); hmapIter != hmap.end(); hmapIter++) { bam1_t *tmp_b = hmapIter->first.b; tmp_b->core.qual = getMAPQ(hmapIter->second); - uint8_t *p = bam1_aux(tmp_b); - *p = 'Z'; ++p; *p = 'W'; ++p; *p = 'f'; ++p; float val = (float)hmapIter->second; - memcpy(p, &val, 4); + bam_aux_append(tmp_b, "ZW", 'f', 4, (uint8_t*)&val); if (tmp_b->core.qual > 0) samwrite(out, tmp_b); // If MAPQ is equal to 0, do not output this alignment bam_destroy1(tmp_b); // now hmapIter->b makes no sense } @@ -367,7 +366,6 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) b->core.qual = b2->core.qual = 255; - //Do I really need this? The insert size uses transcript coordinates if (transcript.getStrand() == '-') { b->core.isize = -b->core.isize; b2->core.isize = -b2->core.isize; @@ -390,27 +388,31 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) push_qual(qual2, readlen2, transcript.getStrand(), data2); free(b->data); - b->m_data = b->data_len = data.size() + 7; // 7 extra bytes for ZW tag - b->l_aux = 7; + b->m_data = b->data_len = data.size(); + b->l_aux = 0; b->data = (uint8_t*)malloc(b->m_data); for (int i = 0; i < b->data_len; i++) b->data[i] = data[i]; free(b2->data); - b2->m_data = b2->data_len = data2.size() + 7; // 7 extra bytes for ZW tag - b2->l_aux = 7; + b2->m_data = b2->data_len = data2.size(); + b2->l_aux = 0; b2->data = (uint8_t*)malloc(b2->m_data); for (int i = 0; i < b2->data_len; i++) b2->data[i] = data2[i]; b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b))); b2->core.bin = bam_reg2bin(b2->core.pos, bam_calend(&(b2->core), bam1_cigar(b2))); + + char strand = transcript.getStrand(); + bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand); + bam_aux_append(b2, "XS", 'A', 1, (uint8_t*)&strand); } else { - b->m_data = b->data_len = b->data_len - b->l_aux + 7; // 7 extra bytes for ZW tag - b->l_aux = 7; + b->m_data = b->data_len = b->data_len - b->l_aux; + b->l_aux = 0; b->data = (uint8_t*)realloc(b->data, b->m_data); - b2->m_data = b2->data_len = b2->data_len - b2->l_aux + 7; // 7 extra bytes for ZW tag - b2->l_aux = 7; + b2->m_data = b2->data_len = b2->data_len - b2->l_aux; + b2->l_aux = 0; b2->data = (uint8_t*)realloc(b2->data, b2->m_data); } @@ -422,13 +424,9 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) tmp_b->core.qual = tmp_b2->core.qual = getMAPQ(hmapIter->second); - uint8_t *p = bam1_aux(tmp_b), *p2 = bam1_aux(tmp_b2); - *p = 'Z'; ++p; *p = 'W'; ++p; *p = 'f'; ++p; - *p2 = 'Z'; ++p2; *p2 = 'W'; ++p2; *p2 = 'f'; ++p2; - float val = (float)hmapIter->second; - memcpy(p, &val, 4); - memcpy(p2, &val, 4); + bam_aux_append(tmp_b, "ZW", 'f', 4, (uint8_t*)&val); + bam_aux_append(tmp_b2, "ZW", 'f', 4, (uint8_t*)&val); // If MAPQ is equal to 0, do not output this alignment pair if (tmp_b->core.qual > 0) { @@ -462,13 +460,9 @@ void BamWriter::work(HitWrapper wrapper, Transcripts& transcripts) tmp_b->core.qual = tmp_b2->core.qual = getMAPQ(hmapIter->second); - uint8_t *p = bam1_aux(tmp_b), *p2 = bam1_aux(tmp_b2); - *p = 'Z'; ++p; *p = 'W'; ++p; *p = 'f'; ++p; - *p2 = 'Z'; ++p2; *p2 = 'W'; ++p2; *p2 = 'f'; ++p2; - float val = (float)hmapIter->second; - memcpy(p, &val, 4); - memcpy(p2, &val, 4); + bam_aux_append(tmp_b, "ZW", 'f', 4, (uint8_t*)&val); + bam_aux_append(tmp_b2, "ZW", 'f', 4, (uint8_t*)&val); if (tmp_b->core.qual > 0) { samwrite(out, tmp_b);