libblasr-dev 0~20151014+gitbe5d1bf-2 / usr / include / blasr / format / BAMPrinterImpl.hpp

#ifndef _BLASR_BAM_PRINTER_IMPL_HPP_
#define _BLASR_BAM_PRINTER_IMPL_HPP_

#ifdef USE_PBBAM

#include <algorithm>
#include "utils/SMRTTitle.hpp"
using namespace BAMOutput;
using namespace std;
#include "pbbam/BamRecord.h"
#include "pbbam/BamFile.h"

template<typename T_Sequence>
void BAMOutput::CreateCIGARString(T_AlignmentCandidate &alignment,
        T_Sequence &read, std::string &cigarString, const bool cigarUseSeqMatch)
{
    cigarString = "";
    // All cigarString use the no clipping core
    std::vector<int> opSize;
    std::vector<char> opChar;

    SAMOutput::CreateNoClippingCigarOps(alignment, opSize, opChar, cigarUseSeqMatch);

    // Clipping needs to be added
    DNALength prefixSoftClip = alignment.QAlignStart() - read.SubreadStart();
    DNALength suffixSoftClip = read.SubreadEnd() - alignment.QAlignEnd();

    if (alignment.tStrand == 1) {
        std::swap(prefixSoftClip, suffixSoftClip);
    }
    if (prefixSoftClip > 0) {
        opSize.insert(opSize.begin(), prefixSoftClip);
        opChar.insert(opChar.begin(), 'S');
    }
    if (suffixSoftClip > 0) {
        opSize.push_back(suffixSoftClip);
        opChar.push_back('S');
    }
    SAMOutput::CigarOpsToString(opSize, opChar, cigarString);
}

template<typename T_Sequence>
void BAMOutput::SetAlignedSequence(T_AlignmentCandidate &alignment, T_Sequence &read,
        T_Sequence &alignedSeq) {
    if (alignment.tStrand == 0) {
        alignedSeq.ReferenceSubstring(read);
    }
    else {
        T_Sequence subSeq;
        subSeq.ReferenceSubstring(read);
        subSeq.MakeRC(alignedSeq);
    }
}

template<typename T_Sequence>
void BAMOutput::AlignmentToBamRecord(T_AlignmentCandidate & alignment,
        T_Sequence & read, PacBio::BAM::BamRecord & bamRecord,
        AlignmentContext & context, SupplementalQVList & qvList,
        Clipping clipping, bool cigarUseSeqMatch) {
    // soft clipping and subread clipping are identical for BAM
    assert(clipping == SAMOutput::soft or clipping == SAMOutput::subread);

    // Build from scratch if input reads are not from pbbam files.
    // Otherwise, use API provided by pbbam to make alignments from bamRecords.
 
    bool buildFromScratch = false;
    if (dynamic_cast<SMRTSequence*>(&read) == NULL) {
        //not SMRTSequence 
        buildFromScratch = true;
    } else { //is SMRTSequence, but not copied from bam.
        if (not read.copiedFromBam) {
            buildFromScratch = true;
        }
    }

    // build cigar string.
    string cigarString;
    BAMOutput::CreateCIGARString(alignment, read, cigarString, cigarUseSeqMatch);
    PacBio::BAM::Cigar cigar = PacBio::BAM::Cigar::FromStdString(cigarString);

    T_Sequence alignedSequence;
    BAMOutput::SetAlignedSequence(alignment, read, alignedSequence);
 
    // build flag
    uint16_t flag;
    BuildFlag(alignment, context, flag);

    // Get sequence string.
    string seqString;
    seqString.assign((char*)alignedSequence.seq, alignedSequence.length);

    // Get alignment starting position on reference sequence forward strand.
    PacBio::BAM::Position pos = 0;
    PacBio::BAM::Strand strand;
    if (alignment.tStrand == 0) {
        pos = static_cast<PacBio::BAM::Position>(alignment.TAlignStart());
        strand = PacBio::BAM::Strand::FORWARD;
    } else {
        pos = static_cast<PacBio::BAM::Position>(alignment.tLength - (alignment.TAlignStart() + alignment.TEnd()));
        strand = PacBio::BAM::Strand::REVERSE;
    }

    if (buildFromScratch) {
        SMRTTitle smrtTitle(alignment.qName);
        if (smrtTitle.isSMRTTitle) {
            bamRecord.Impl().Name(smrtTitle.ToString());
        } else {
            cout << "ERROR, can not convert non-pacbio reads to pbbam record." << endl;
            exit(-1);
        }
        bamRecord.Impl().SetSequenceAndQualities(seqString, alignedSequence.qual.ToString());
        bamRecord.Impl().CigarData(cigar);
        bamRecord.Impl().Bin(0);
        bamRecord.Impl().InsertSize(0);
        bamRecord.Impl().MapQuality(static_cast<uint8_t>(alignment.mapQV));
        bamRecord.Impl().MatePosition(static_cast<PacBio::BAM::Position>(-1));
        bamRecord.Impl().MateReferenceId(static_cast<int32_t>(-1));
        bamRecord.Impl().Position(pos);
        bamRecord.Impl().ReferenceId(static_cast<int32_t>(alignment.tIndex));

        if (strand == PacBio::BAM::Strand::REVERSE) {
            bamRecord.Impl().SetReverseStrand(true);
        }

        // Add tags required for bax->bam.
        PacBio::BAM::TagCollection tags;
        tags["RG"] = context.readGroupId;
        if (dynamic_cast<CCSSequence*>(&read) == NULL) { // subread
            tags["qs"] = read.SubreadStart();
            tags["qe"] = read.SubreadEnd();
            /// Number of passes for a subread should always be 1.
            tags["np"] = 1;
        } else { // ccs read
            /// Number of passes for ccs reads.
            tags["np"] = (static_cast<CCSSequence*>(&read))->numPasses;
        }
        tags["zm"] = read.zmwData.holeNumber;

        // Build QV tags.
        // Skip tags not define in BAM specification 3.0. 
        // including XL, XT, XQ, XS, XE, YS, YE.
        //
        // Write out optional quality values.  If qvlist does not
        // have any qv's signaled to print, this is a no-op.
        //
        // First transform characters that are too large to printable ones.
        qvList.FormatQVOptionalFields(alignedSequence);

        // Add QVs to BamRecordImpl.
        string insertionQVs, deletionQVs, substitutionQVs, mergeQVs, substitutionTags, deletionTags;
        bool alnReverse = (alignment.tStrand == 1); // reverse-complement alignment
        // If this is a reverse-complement alignment, bases and QVs of 
        // alignedSequence are reverse(-complement) of the sequence read
        // from bax.h5 file. With PB BAM specification 3.0, QVs will be stored 
        // in BAM Record, therefore, reverse orders of QVs before add QV tags. 
        if (alignedSequence.GetQVs("InsertionQV", insertionQVs, alnReverse)) { 
            tags["iq"] = insertionQVs;
        }
        if (alignedSequence.GetQVs("DeletionQV", deletionQVs, alnReverse)) {
            tags["dq"] = deletionQVs;
        } 
        if (alignedSequence.GetQVs("SubstitutionQV", substitutionQVs, alnReverse)) {
            tags["sq"] = substitutionQVs;
        } 
        if (alignedSequence.GetQVs("MergeQV", mergeQVs, alnReverse)) {
            tags["mq"] = mergeQVs;
        }
        // substitutionTag is not included by default
        if (alignedSequence.GetQVs("DeletionTag", deletionTags, alnReverse)) {
            tags["dt"] = deletionTags;
        }
        bamRecord.Impl().Tags(tags);
    } else {
        // The following code can be used to hard-clip reads, if needed. 
        // PacBio::BAM::Position clipStart = read.bamRecord.QueryStart() + alignment.QAlignStart();
        // PacBio::BAM::Position clipEnd = read.bamRecord.QueryStart() + alignment.QAlignEnd();
        // bamRecord = PacBio::BAM::BamRecord::Clipped(read.bamRecord, 
        //                PacBio::BAM::ClipType::CLIP_TO_QUERY,
        //                clipStart, clipEnd).
        bamRecord = PacBio::BAM::BamRecord::Mapped(read.bamRecord,
                        static_cast<int32_t>(alignment.tIndex),
                        static_cast<PacBio::BAM::Position>(pos),
                        strand, cigar, 
                        static_cast<uint8_t>(alignment.mapQV));
    }

    // Add tags common for bax->bam and bam->bam.
    bamRecord.Impl().AddTag("AS", alignment.score);
    bamRecord.Impl().AddTag("NM", context.editDist);

    // Set Flag 
    bamRecord.Impl().Flag(static_cast<uint32_t>(flag));
}

template<typename T_Sequence>
void BAMOutput::PrintAlignment(T_AlignmentCandidate &alignment, T_Sequence &read,
        PacBio::BAM::BamWriter &bamWriter, AlignmentContext &context, 
        SupplementalQVList & qvList, Clipping clipping, bool cigarUseSeqMatch) {

    PacBio::BAM::BamRecord bamRecord;
    BAMOutput::AlignmentToBamRecord(alignment, read, bamRecord, context, qvList, clipping, cigarUseSeqMatch);
    bamWriter.Write(bamRecord);
}
#endif

#endif