C++ BamAlignment类代码示例

ReadGroup::ReadGroup(BamAlignment &al, int max_isize, int isize_samples,
	string prefix, list<string> blacklist) :
	max_isize(max_isize),
	isize_samples(isize_samples),
	prefix(prefix),
	blacklisted(false)
{
	if (!al.GetReadGroup(name))
		name = "none";

	nreads = 0;

	/* Determine if this read group is in the blacklist */
	for (list<string>::iterator it = blacklist.begin();
	     it != blacklist.end(); ++it) {
		if (*it == name) {
			blacklisted = true;
			break;
		}
	}

	if (!blacklisted) {
		f1.open((prefix + "/" + name + "_1.fq.gz").c_str());
		f2.open((prefix + "/" + name + "_2.fq.gz").c_str());
	}

	witness(al);
}

void
clipAlignment(BamAlignment &al)
{
	int offset, length;
	CigarOp cop1 = al.CigarData[0];
	CigarOp cop2 = al.CigarData[al.CigarData.size() - 1];

	if (copcomp(cop2, cop1)) {
		offset = 0;
		length = min(al.Length, (signed)cop1.Length);
	} else {
		offset = al.Length - min(al.Length, (signed)cop2.Length);
		length = min(al.Length, (signed)cop2.Length);
	}

	try {
		al.Qualities = al.Qualities.substr(offset, length);
		al.QueryBases = al.QueryBases.substr(offset, length);
	} catch (exception &e) {
		cout << "ERROR: substr failed in clipAlignment()" << endl;
		cout << al.Name << " " << (al.IsReverseStrand() ? "(-)" : "(+)");
		cout << " offset: " << offset << " length: " << length
		     << " taglen: " << al.Length << endl;
		cout << "cop1: " << cop1.Length << cop1.Type << endl;
		cout << "cop2: " << cop2.Length << cop2.Type << endl;
		exit(1);
	}
}

pos_t VariantProcessor::processMatchOrMismatch(const BamAlignment& alignment, 
					       vector<VariantPtr>& read_variants, 
					       const uint32_t& op_length, const string& refseq, 
					       const pos_t& refpos, const pos_t& readpos) {
  // Process a matching or mismatching sequence in the CIGAR string,
  // adding any SNP variants present.
  int endpos = alignment.GetEndPosition();
  for (int i = 0; i < op_length; i++) {
    assert(alignment.Position + i < endpos);
    char query_base = alignment.QueryBases[readpos + i];
    assert(refpos + i < refseq.size());
    char ref_base = refseq[refpos + i];
    if (ref_base != query_base) {
      // SNP
      string ref(1, ref_base), alt(1, query_base);
      char qual_base = alignment.Qualities[refpos + i]; // TODO check

      VariantPtr snp(new Variant(VariantType::SNP, alignment.RefID, 
				 alignment.Position+i, 1, 0, ref, alt));
      block_variants.insert(snp);
      read_variants.push_back(snp);
      //cout << "mismatch at " << alignment.Position + i <<" refbase: " << ref_base << " querybase: " << query_base << endl;      
    }
  }
}

int main (int argc, char *argv[]) {

     if( (argc== 1) ||
    	(argc== 2 && string(argv[1]) == "-h") ||
    	(argc== 2 && string(argv[1]) == "-help") ||
    	(argc== 2 && string(argv[1]) == "--help") ){
	 cout<<"Usage:setAsUnpaired [in bam] [outbam]"<<endl<<"this program takes flags all paired sequences as singles"<<endl;
    	return 1;
    }

     string bamfiletopen = string(argv[1]);
     string bamFileOUT   = string(argv[2]);

     BamReader reader;
     BamWriter writer;

     if ( !reader.Open(bamfiletopen) ) {
    	cerr << "Could not open input BAM files." << endl;
    	return 1;
     }
    const SamHeader header = reader.GetHeader();
    const RefVector references = reader.GetReferenceData();
    if ( !writer.Open(bamFileOUT,header,references) ) {
    	cerr << "Could not open output BAM file "<<bamFileOUT << endl;
    	return 1;
    }

    BamAlignment al;
 
    while ( reader.GetNextAlignment(al) ) {
	if(al.IsMapped()){
	    cerr << "Cannot yet handle mapped reads " << endl;
	    return 1;
	}

	
	al.SetIsPaired (false);
	
	writer.SaveAlignment(al);    

    } //while al

    reader.Close();
    writer.Close();

    return 0;
}

int main (int argc, char *argv[]) {

     if( (argc== 1) ||
    	(argc== 2 && string(argv[1]) == "-h") ||
    	(argc== 2 && string(argv[1]) == "-help") ||
    	(argc== 2 && string(argv[1]) == "--help") ){
	 cout<<"Usage:editDist [in bam]"<<endl<<"this program returns the NM field of all aligned reads"<<endl;
	 return 1;
     }

     string bamfiletopen = string(argv[1]);
     // cout<<bamfiletopen<<endl;
     BamReader reader;
     // cout<<"ok"<<endl;
     if ( !reader.Open(bamfiletopen) ) {
	 cerr << "Could not open input BAM files." << endl;
	 return 1;
     }

     BamAlignment al;
     // cout<<"ok"<<endl;
     while ( reader.GetNextAlignment(al) ) {
	 // cout<<al.Name<<endl;
	 if(!al.IsMapped())
	     continue;

	 if(al.HasTag("NM") ){
	     int editDist;
	     if(al.GetTag("NM",editDist) ){
		 cout<<editDist<<endl;
	     }else{
		 cerr<<"Cannot retrieve NM field for "<<al.Name<<endl;
		 return 1;
	     }
	 }else{
	     cerr<<"Warning: read "<<al.Name<<" is aligned but has no NM field"<<endl;
	 }

		    

     } //while al

     reader.Close();

     return 0;
}

int VariantProcessor::run() {

  int nmapped = 0;
  last_aln_pos = 0;
  bool stop;
  BamAlignment al; // TODO: mem copying issues?
  
  if (!reader.IsOpen()) {
    std::cerr << "error: BAM file '" << filename << "' is not open." << std::endl;
  }

  while (reader.GetNextAlignment(al)) {
    if (!al.IsMapped()) continue;

    // TODO add chromosome checking code here

    assert(al.Position >= last_aln_pos); // ensure is sorted

    if (al.Position > block_start) {
      // only check if we can stop if we've moved in the block
      stop = isBlockEnd(al);
    }

    if (stop) {
      // end of block; process all variants in this block and output
      // haplotype count statistics
      processBlockAlignments();
      
      // reset block
      blockReset((pos_t) al.Position);
      stop = false;
    } else {
      // process read
      last_aln_pos = processAlignment(al);
    }
    
    // for debug TODO
    for (set<VariantPtr>::const_iterator it = block_variants.begin(); it != block_variants.end(); ++it) {
      (*it)->print();
    }

    nmapped++;
  }

  return nmapped;
}

// returns region state - whether alignment ends before, overlaps, or starts after currently specified region
// this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
BamReaderPrivate::RegionState BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {

    // if alignment is on any reference sequence before left bound
    if ( bAlignment.RefID < Region.LeftRefID ) return BEFORE_REGION;

    // if alignment starts on left bound reference
    else if ( bAlignment.RefID == Region.LeftRefID ) {

	// if alignment starts at or after left boundary
	if ( bAlignment.Position >= Region.LeftPosition) {

	    // if right boundary is specified AND
	    // left/right boundaries are on same reference AND
	    // alignment starts past right boundary
	    if ( Region.isRightBoundSpecified() &&
		 Region.LeftRefID == Region.RightRefID &&
		 bAlignment.Position > Region.RightPosition )
		return AFTER_REGION;

	    // otherwise, alignment is within region
	    return WITHIN_REGION;
	}

	// alignment starts before left boundary
	else {
	    // check if alignment overlaps left boundary
	    if ( bAlignment.GetEndPosition() >= Region.LeftPosition ) return WITHIN_REGION;
	    else return BEFORE_REGION;
	}
    }

    // alignment starts on a reference after the left bound
    else {

	// if region has a right boundary
	if ( Region.isRightBoundSpecified() ) {

	    // alignment is on reference between boundaries
	    if ( bAlignment.RefID < Region.RightRefID ) return WITHIN_REGION;

	    // alignment is on reference after right boundary
	    else if ( bAlignment.RefID > Region.RightRefID ) return AFTER_REGION;

	    // alignment is on right bound reference
	    else {
		// check if alignment starts before or at right boundary
		if ( bAlignment.Position <= Region.RightPosition ) return WITHIN_REGION;
		else return AFTER_REGION;
	    }
	}

	// otherwise, alignment is after left bound reference, but there is no right boundary
	else return WITHIN_REGION;
    }
}

    void getBamBlocks(const BamAlignment &bam, const RefVector &refs,
                      vector<BED> &blocks, bool breakOnDeletionOps) {

        CHRPOS currPosition = bam.Position;
        CHRPOS blockStart   = bam.Position;
        string chrom        = refs.at(bam.RefID).RefName;
        string name         = bam.Name;
        string strand       = "+";
        string score        = ToString(bam.MapQuality);
        char  prevOp        = '\0';
        if (bam.IsReverseStrand()) strand = "-";
        bool blocksFound = false;

        vector<CigarOp>::const_iterator cigItr = bam.CigarData.begin();
        vector<CigarOp>::const_iterator cigEnd = bam.CigarData.end();
        for ( ; cigItr != cigEnd; ++cigItr ) {
            if (cigItr->Type == 'M') {
                currPosition += cigItr->Length;
                // we only want to create a new block if the current M op
                // was preceded by an N op or a D op (and we are breaking on D ops)
                if ((prevOp == 'D' && breakOnDeletionOps == true) || (prevOp == 'N')) {
                    blocks.push_back( BED(chrom, blockStart, currPosition, name, score, strand) );
                    blockStart = currPosition;
                }
            }
            else if (cigItr->Type == 'D') {
                if (breakOnDeletionOps == false)
                    currPosition += cigItr->Length;
                else {
                    blocksFound = true;
                    currPosition += cigItr->Length;
                    blockStart    = currPosition;
                }
            }
            else if (cigItr->Type == 'N') {
                blocks.push_back( BED(chrom, blockStart, currPosition, name, score, strand) );
                blocksFound = true;
                currPosition += cigItr->Length;
                blockStart    = currPosition;
            }
            else if (cigItr->Type == 'S' || cigItr->Type == 'H' || cigItr->Type == 'P' || cigItr->Type == 'I') {
                // do nothing
            }
            else {
                cerr << "Input error: invalid CIGAR type (" << cigItr->Type
                    << ") for: " << bam.Name << endl;
                exit(1);
            }
            prevOp = cigItr->Type;
        }
        // if there were no splits, we just create a block representing the contiguous alignment.
        if (blocksFound == false) {
            blocks.push_back( BED(chrom, bam.Position, currPosition, name, score, strand) );
        }
    }

int DataStatisticsTool::Execute()
{
    // iterate over reads in BAM file(s)
    BamAlignment alignObj;
    while(bamReader.GetNextAlignment(alignObj))
    {
        if (alignObj.IsDuplicate()) continue;
        if (alignObj.IsFailedQC()) continue;
        if (!alignObj.IsMapped()) continue;
        if (!alignObj.IsPrimaryAlignment()) continue;
        if (alignObj.IsPaired() && !alignObj.IsProperPair()) continue;
        if (alignObj.IsPaired() && !alignObj.IsMateMapped()) continue;
        if (!alignObj.HasTag("MD")) continue;

//        // debug
//        GenericBamAlignmentTools::printBamAlignmentCigar(alignObj);
//        GenericBamAlignmentTools::printBamAlignmentMD(alignObj);

        // shift InDel
        GenericBamAlignmentTools::leftShiftInDel(alignObj);

//        // debug
//        GenericBamAlignmentTools::printBamAlignmentCigar(alignObj);
//        GenericBamAlignmentTools::printBamAlignmentMD(alignObj);

        // get the alignment sequences
        string alignRead;
        string alignGenome;
        GenericBamAlignmentTools::getAlignmentSequences(alignObj, alignRead, alignGenome);

        // update the statistics
        statistics.update(alignRead, alignGenome);
    }


    // print to screen
    cout << statistics << endl;
//    statistics.printMatchMismatch();

    // close BAM reader
    bamReader.Close();

    // close Fasta
    genomeFasta.Close();

    return 1;
}

// get next alignment (with character data fully parsed)
bool BamReaderPrivate::GetNextAlignment(BamAlignment& alignment) {

    // if valid alignment found
    if ( GetNextAlignmentCore(alignment) ) {

        // store alignment's "source" filename
        alignment.Filename = m_filename;

        // return success/failure of parsing char data
        if ( alignment.BuildCharData() )
            return true;
        else {
            const string alError = alignment.GetErrorString();
            const string message = string("could not populate alignment data: \n\t") + alError;
            SetErrorString("BamReader::GetNextAlignment", message);
            return false;
        }
    }

    // no valid alignment found
    return false;
}

//{{{ void process_pair(const BamAlignment &curr,
void
SV_Pair::
process_pair(const BamAlignment &curr,
             const RefVector refs,
             map<string, BamAlignment> &mapped_pairs,
             UCSCBins<SV_BreakPoint*> &r_bin,
             int weight,
             int ev_id,
             SV_PairReader *reader)
{
    if (mapped_pairs.find(curr.Name) == mapped_pairs.end())
        mapped_pairs[curr.Name] = curr;
    else {
        SV_Pair *new_pair = new SV_Pair(mapped_pairs[curr.Name],
                                        curr,
                                        refs,
                                        weight,
                                        ev_id,
                                        reader);
        //cerr << count_clipped(curr.CigarData) << "\t" <<
                //count_clipped(mapped_pairs[curr.Name].CigarData) << endl;
                
        if ( new_pair->is_sane() &&  
             new_pair->is_aberrant() &&
             (count_clipped(curr.CigarData) > 0) &&
             (count_clipped(mapped_pairs[curr.Name].CigarData) > 0) ) {
            SV_BreakPoint *new_bp = new_pair->get_bp();

#ifdef TRACE

            cerr << "READ\t" << 
                    refs.at(mapped_pairs[curr.Name].RefID).RefName << "," <<
                    mapped_pairs[curr.Name].Position << "," <<
                    (mapped_pairs[curr.Name].GetEndPosition(false, false) - 1)
                        << "\t" <<
                    refs.at(curr.RefID).RefName << "," <<
                    curr.Position << "," <<
                    (curr.GetEndPosition(false, false) - 1)
                        <<
                    endl;

            cerr << "\tPE\t" << *new_bp << endl;
#endif
            new_bp->cluster(r_bin);
        } else {
            delete(new_pair);
        }

        mapped_pairs.erase(curr.Name);
    }
}

string createReferenceSequence(const BamAlignment& alignment) {
  // Recreate a reference sequence for a particular alignment. This is
  // the reference sequence that is identical to the reference at this
  // spot. This means skipping insertions or soft clipped regions in
  // reads, adding deletions back in, and keeping read matches.
  const vector<CigarOp> cigar = alignment.CigarData;
  const string querybases = alignment.QueryBases;
  string md_tag;
  alignment.GetTag("MD", md_tag);
  
  vector<MDToken> tokens;
  string refseq, alignedseq; // final ref bases; aligned portion of ref bases
  int md_len = TokenizeMD(md_tag, tokens);

  // Create reference-aligned sequence of read; doesn't contain soft
  // clips or insertions. Then, deletions and reference alleles are
  // added onto this.
  int pos=0;
  for (vector<CigarOp>::const_iterator op = cigar.begin(); op != cigar.end(); ++op) {
    if (!(op->Type == 'S' || op->Type == 'I')) {
      alignedseq.append(querybases.substr(pos, op->Length));
      pos += op->Length;
    } else {
      pos += op->Length; // increment read position past skipped bases
    }
  }

  // the size of the aligned sequence MUST equal what is returned from
  // TokenizeMD: the number of aligned bases. Not the real reference
  // sequence is this length + deletions, which we add in below.
  assert(alignedseq.size() == md_len);

  pos = 0;
  for (vector<MDToken>::const_iterator it = tokens.begin(); it != tokens.end(); ++it) {
    if (it->type == MDType::isMatch) {
      refseq.append(alignedseq.substr(pos, it->length));
      pos += it->length;
    } else if (it->type == MDType::isSNP) {
      assert(it->length == it->seq.size());
      refseq.append(it->seq);
      pos += it->length;
    } else if (it->type == MDType::isDel) {
      // does not increment position in alignedseq
      assert(it->length == it->seq.size());
      refseq.append(it->seq);
    } else {
      assert(false);
    }
  }
  return refseq;
}

void CountDepth(Histogram& hist, BamMultiReader& reader, BamAlignment& al, int32_t refID, int64_t refLen)
{
    bool moreReads = (al.RefID == refID);

    int32_t maxReadLen = 1000;
    vector<int64_t> readEnds(maxReadLen);

    int64_t depth = 0;
    for(int64_t pos=0; pos<refLen; ++pos){
        while(moreReads and al.Position == pos){
            ++depth;
            assert(al.GetEndPosition() - pos < maxReadLen);
            ++readEnds[al.GetEndPosition() % maxReadLen];
            moreReads = GetNextAlignment(al, reader, refID);
        }
        depth -= readEnds[pos % maxReadLen];
        assert(depth >= 0);
        readEnds[pos % maxReadLen] = 0;
        if(depth >= hist.size())
            hist.resize(2 * depth);
        ++hist[depth];
    }
}

// get next alignment (with character data fully parsed)
bool BamReaderPrivate::GetNextAlignment(BamAlignment& alignment) {

    // if valid alignment found
    if ( GetNextAlignmentCore(alignment) ) {

        // store alignment's "source" filename
        alignment.Filename = m_filename;

        // return success/failure of parsing char data
        return alignment.BuildCharData();
    }

    // no valid alignment found
    return false;
}

bool 
processReadPair(const BamAlignment& al1, 
        const BamAlignment& al2, 
        const RefVector& refs, 
        const int32_t totalTail, 
        const int32_t critTail, 
        const bool diff_ref)
{
    if ((al1.IsFirstMate() && al2.IsFirstMate())
        || (al1.IsSecondMate() && al2.IsSecondMate())) {
        cerr << "Incompatible mate orders: name1 = " << al1.Name 
             << " is1stmate " << al1.IsFirstMate() << " is2ndmate " << al1.IsSecondMate()
             << " name2 = " << al2.Name 
             << " is1stmate " << al2.IsFirstMate() << " is2ndmate " << al2.IsSecondMate()
             << endl;
        exit(1);
    }

    int32_t total_tail = -1;
    if (! (total_tail = checkLinkPair(al1, al2, refs, totalTail, critTail, diff_ref))) {
        return false;  // reject all but link pairs
        // continue;
    }
    if (critTail && ! checkLinkPairCandidate(al1, refs, critTail)
        && ! checkLinkPairCandidate(al2, refs, critTail)) {
        return false;  // neither read was a link pair candidate
    }
    if (debug_processReadPair) cout << "---------------------------------" << endl;
    int32_t lpc_tail1 = checkLinkPairCandidate(al1, refs, critTail);
    int32_t lpc_tail2 = checkLinkPairCandidate(al2, refs, critTail);
    if (debug_processReadPair) {
        printAlignmentInfo(al1, refs);
        if (lpc_tail1) {
            cout << "LINK PAIR CANDIDATE ";
            cout << ((lpc_tail1 > 0) ? "--->" : "<---") << " " << lpc_tail1 << endl;
        }
        printAlignmentInfo(al2, refs);
        if (lpc_tail2) {
            cout << "LINK PAIR CANDIDATE ";
            cout << ((lpc_tail2 > 0) ? "--->" : "<---") << " " << lpc_tail2 << endl;
        }
        cout << "TOTAL TAIL " << (abs(readTail(al1, refs)) + abs(readTail(al2, refs))) << endl;
    }

    return true;
}