// BamAlignment.cpp (c) 2009 Derek Barnett
// Marth Lab, Department of Biology, Boston College
// ---------------------------------------------------------------------------
-// Last modified: 10 October 2011 (DB)
+// Last modified: 4 December 2012 (DB)
// ---------------------------------------------------------------------------
// Provides the BamAlignment data structure
// ***************************************************************************
-#include <api/BamAlignment.h>
-#include <api/BamConstants.h>
+#include "api/BamAlignment.h"
+#include "api/BamConstants.h"
using namespace BamTools;
using namespace std;
*/
/*! \var BamAlignment::QueryBases
\brief 'original' sequence (as reported from sequencing machine)
+
+ \note Setting this field to "*" indicates that the sequence is not to be stored on output.
+ In this case, the contents of the Qualities field should be invalidated as well (cleared or marked as "*").
*/
/*! \var BamAlignment::AlignedBases
\brief 'aligned' sequence (includes any indels, padding, clipping)
+
+ This field will be completely empty after reading from BamReader/BamMultiReader when
+ QueryBases is empty.
*/
/*! \var BamAlignment::Qualities
\brief FASTQ qualities (ASCII characters, not numeric values)
+
+ \note Setting this field to "*" indicates to BamWriter that the quality scores are not to be stored,
+ but instead will be output as a sequence of '0xFF'. Otherwise, QueryBases must not be a "*" and
+ the length of this field should equal the length of QueryBases.
*/
/*! \var BamAlignment::TagData
\brief tag data (use the provided methods to query/modify)
\brief constructor
*/
BamAlignment::BamAlignment(void)
- : RefID(-1)
+ : Length(0)
+ , RefID(-1)
, Position(-1)
+ , Bin(0)
+ , MapQuality(0)
+ , AlignmentFlag(0)
, MateRefID(-1)
, MatePosition(-1)
, InsertSize(0)
const unsigned int tagDataLength = dataLength - tagDataOffset;
// check offsets to see what char data exists
- const bool hasSeqData = ( seqDataOffset < dataLength );
- const bool hasQualData = ( qualDataOffset < dataLength );
+ const bool hasSeqData = ( seqDataOffset < qualDataOffset );
+ const bool hasQualData = ( qualDataOffset < tagDataOffset );
const bool hasTagData = ( tagDataOffset < dataLength );
- // set up char buffers
- const char* allCharData = SupportData.AllCharData.data();
- const char* seqData = ( hasSeqData ? (((const char*)allCharData) + seqDataOffset) : (const char*)0 );
- const char* qualData = ( hasQualData ? (((const char*)allCharData) + qualDataOffset) : (const char*)0 );
- char* tagData = ( hasTagData ? (((char*)allCharData) + tagDataOffset) : (char*)0 );
-
// store alignment name (relies on null char in name as terminator)
- Name.assign((const char*)(allCharData));
+ Name.assign(SupportData.AllCharData.data());
// save query sequence
QueryBases.clear();
if ( hasSeqData ) {
+ const char* seqData = SupportData.AllCharData.data() + seqDataOffset;
QueryBases.reserve(SupportData.QuerySequenceLength);
for ( size_t i = 0; i < SupportData.QuerySequenceLength; ++i ) {
const char singleBase = Constants::BAM_DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
}
}
- // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
+ // save qualities
+
Qualities.clear();
if ( hasQualData ) {
- Qualities.reserve(SupportData.QuerySequenceLength);
- for ( size_t i = 0; i < SupportData.QuerySequenceLength; ++i ) {
- const char singleQuality = static_cast<const char>(qualData[i]+33);
- Qualities.append(1, singleQuality);
+ const char* qualData = SupportData.AllCharData.data() + qualDataOffset;
+
+ // if marked as unstored (sequence of 0xFF) - don't do conversion, just fill with 0xFFs
+ if ( qualData[0] == (char)0xFF )
+ Qualities.resize(SupportData.QuerySequenceLength, (char)0xFF);
+
+ // otherwise convert from numeric QV to 'FASTQ-style' ASCII character
+ else {
+ Qualities.reserve(SupportData.QuerySequenceLength);
+ for ( size_t i = 0; i < SupportData.QuerySequenceLength; ++i )
+ Qualities.append(1, qualData[i]+33);
}
}
// if QueryBases has data, build AlignedBases using CIGAR data
// otherwise, AlignedBases will remain empty (this case IS allowed)
- if ( !QueryBases.empty() ) {
+ if ( !QueryBases.empty() && QueryBases != "*" ) {
// resize AlignedBases
AlignedBases.reserve(SupportData.QuerySequenceLength);
// save tag data
TagData.clear();
if ( hasTagData ) {
+
+ char* tagData = (((char*)SupportData.AllCharData.data()) + tagDataOffset);
+
if ( IsBigEndian ) {
size_t i = 0;
while ( i < tagDataLength ) {
return false;
}
+/*! \fn bool BamAlignment::GetArrayTagType(const std::string& tag, char& type) const
+ \brief Retrieves the BAM tag type-code for the array elements associated with requested tag name.
+
+ \param[in] tag 2-character tag name
+ \param[out] type retrieved (1-character) type-code
+
+ \return \c true if found. False if not found, or if tag is not an array type.
+ \sa \samSpecURL for more details on reserved tag names, supported tag types, etc.
+*/
+bool BamAlignment::GetArrayTagType(const std::string& tag, char& type) const {
+
+ // skip if alignment is core-only
+ if ( SupportData.HasCoreOnly ) {
+ // TODO: set error string?
+ return false;
+ }
+
+ // skip if no tags present
+ if ( TagData.empty() ) {
+ // TODO: set error string?
+ return false;
+ }
+
+ // localize the tag data
+ char* pTagData = (char*)TagData.data();
+ const unsigned int tagDataLength = TagData.size();
+ unsigned int numBytesParsed = 0;
+
+ // if tag not found, return failure
+ if ( !FindTag(tag, pTagData, tagDataLength, numBytesParsed) ){
+ // TODO: set error string?
+ return false;
+ }
+
+ // check that tag type code is array
+ type = *(pTagData - 1);
+ if ( type != Constants::BAM_TAG_TYPE_ARRAY ) {
+ // TODO: set error string
+ return false;
+ }
+
+ // fetch element type
+ const char elementType = *pTagData;
+ switch ( elementType ) {
+
+ // allowable types
+ case (Constants::BAM_TAG_TYPE_INT8) :
+ case (Constants::BAM_TAG_TYPE_UINT8) :
+ case (Constants::BAM_TAG_TYPE_INT16) :
+ case (Constants::BAM_TAG_TYPE_UINT16) :
+ case (Constants::BAM_TAG_TYPE_INT32) :
+ case (Constants::BAM_TAG_TYPE_UINT32) :
+ case (Constants::BAM_TAG_TYPE_FLOAT) :
+ type = elementType;
+ break;
+
+ default:
+ //TODO: set error string
+ return false;
+ }
+
+ // if we get here, return success
+ return true;
+}
+
+
/*! \fn int BamAlignment::GetEndPosition(bool usePadded = false, bool closedInterval = false) const
\brief Calculates alignment end position, based on its starting position and CIGAR data.
return ErrorString;
}
+/*! \fn bool BamAlignment::GetSoftClips(std::vector<int>& clipSizes, std::vector<int>& readPositions, std::vector<int>& genomePositions, bool usePadded = false) const
+ \brief Identifies if an alignment has a soft clip. If so, identifies the
+ sizes of the soft clips, as well as their positions in the read and reference.
+
+ \param[out] clipSizes vector of the sizes of each soft clip in the alignment
+ \param[out] readPositions vector of the 0-based read locations of each soft clip in the alignment.
+ These positions are basically indexes within the read, not genomic positions.
+ \param[out] genomePositions vector of the 0-based genome locations of each soft clip in the alignment
+ \param[in] usePadded inserted bases affect reported position. Default is false, so that
+ reported position stays 'sync-ed' with reference coordinates.
+
+ \return \c true if any soft clips were found in the alignment
+*/
+bool BamAlignment::GetSoftClips(vector<int>& clipSizes,
+ vector<int>& readPositions,
+ vector<int>& genomePositions,
+ bool usePadded) const
+{
+ // initialize positions & flags
+ int refPosition = Position;
+ int readPosition = 0;
+ bool softClipFound = false;
+ bool firstCigarOp = true;
+
+ // iterate over cigar operations
+ vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
+ vector<CigarOp>::const_iterator cigarEnd = CigarData.end();
+ for ( ; cigarIter != cigarEnd; ++cigarIter) {
+ const CigarOp& op = (*cigarIter);
+
+ switch ( op.Type ) {
+
+ // increase both read & genome positions on CIGAR chars [DMXN=]
+ case Constants::BAM_CIGAR_DEL_CHAR :
+ case Constants::BAM_CIGAR_MATCH_CHAR :
+ case Constants::BAM_CIGAR_MISMATCH_CHAR :
+ case Constants::BAM_CIGAR_REFSKIP_CHAR :
+ case Constants::BAM_CIGAR_SEQMATCH_CHAR :
+ refPosition += op.Length;
+ readPosition += op.Length;
+ break;
+
+ // increase read position on insertion, genome position only if @usePadded is true
+ case Constants::BAM_CIGAR_INS_CHAR :
+ readPosition += op.Length;
+ if ( usePadded )
+ refPosition += op.Length;
+ break;
+
+ case Constants::BAM_CIGAR_SOFTCLIP_CHAR :
+
+ softClipFound = true;
+
+ //////////////////////////////////////////////////////////////////////////////
+ // if we are dealing with the *first* CIGAR operation
+ // for this alignment, we increment the read position so that
+ // the read and genome position of the clip are referring to the same base.
+ // For example, in the alignment below, the ref position would be 4, yet
+ // the read position would be 0. Thus, to "sync" the two,
+ // we need to increment the read position by the length of the
+ // soft clip.
+ // Read: ATCGTTTCGTCCCTGC
+ // Ref: GGGATTTCGTCCCTGC
+ // Cigar: SSSSMMMMMMMMMMMM
+ //
+ // NOTE: This only needs to be done if the soft clip is the _first_ CIGAR op.
+ //////////////////////////////////////////////////////////////////////////////
+ if ( firstCigarOp )
+ readPosition += op.Length;
+
+ // track the soft clip's size, read position, and genome position
+ clipSizes.push_back(op.Length);
+ readPositions.push_back(readPosition);
+ genomePositions.push_back(refPosition);
+
+ // any other CIGAR operations have no effect
+ default :
+ break;
+ }
+
+ // clear our "first pass" flag
+ firstCigarOp = false;
+ }
+
+ // return whether any soft clips found
+ return softClipFound;
+}
+
+/*! \fn std::vector<std::string> BamAlignment::GetTagNames(void) const
+ \brief Retrieves the BAM tag names.
+
+ When paired with GetTagType() and GetTag(), this method allows you
+ to iterate over an alignment's tag data without knowing the names (or types)
+ beforehand.
+
+ \return \c vector containing all tag names found (empty if none available)
+ \sa \samSpecURL for more details on reserved tag names, supported tag types, etc.
+*/
+std::vector<std::string> BamAlignment::GetTagNames(void) const {
+
+ std::vector<std::string> result;
+ if ( SupportData.HasCoreOnly || TagData.empty() )
+ return result;
+
+ char* pTagData = (char*)TagData.data();
+ const unsigned int tagDataLength = TagData.size();
+ unsigned int numBytesParsed = 0;
+ while ( numBytesParsed < tagDataLength ) {
+
+ // get current tag name & type
+ const char* pTagName = pTagData;
+ const char* pTagType = pTagData + 2;
+ pTagData += 3;
+ numBytesParsed +=3;
+
+ // store tag name
+ result.push_back( std::string(pTagName, 2) );
+
+ // find the next tag
+ if ( *pTagType == '\0' ) break;
+ if ( !SkipToNextTag(*pTagType, pTagData, numBytesParsed) ) break;
+ if ( *pTagData == '\0' ) break;
+ }
+
+ return result;
+}
+
/*! \fn bool BamAlignment::GetTagType(const std::string& tag, char& type) const
\brief Retrieves the BAM tag type-code associated with requested tag name.
projects / bamtools.git / blobdiff