snippet from multi_bam_cov

multi_bam_cov.cpp

void MultiCovBam::CollectCoverage()
{
    BamMultiReader reader;

    if ( !reader.Open(_bam_files) )
    {
        cerr << "Could not open input BAM files." << endl; return;
    }
    else
    {
        // attempt to find index files
        reader.LocateIndexes();

        // if index data available for all BAM files, we can use SetRegion
        if ( reader.HasIndexes() ) {
            BED bed, nullBed;
            int lineNum = 0;
            BedLineStatus bedStatus;

            _bed->Open();
            // loop through each BED entry, jump to it, 
            // and collect coverage from each BAM
            while ((bedStatus = _bed->GetNextBed(bed, lineNum)) != BED_INVALID)
            {

                if (bedStatus == BED_VALID)
                {
                    // attempt to set region on reader
                    if ( !reader.SetRegion(reader.GetReferenceID(bed.chrom),
                                           (int) bed.start,
                                           reader.GetReferenceID(bed.chrom),
                                           (int) bed.end) )
                    {
                        cerr << "ERROR: set region failed. Check that REGION describes a valid range" << endl;
                        reader.Close();
                        exit(1);
                    }

                    // everything checks out, just iterate through specified region, counting alignments
                    vector<int> counts(_bam_files.size());
                    BamAlignment al;
                    while ( reader.GetNextAlignmentCore(al) )
                    {
                        // lookup the offset of the file name and tabulate coverage
                        // for the appropriate file
                        counts[bamFileMap[al.Filename]]++;
                    }
                    // report the cov at this interval for each file and reset
                    _bed->reportBedTab(bed);
                    ReportCounts(counts);
                    bed = nullBed;
                }
            }
            _bed->Close();
        }
        else {
            cerr << "Could not find indexes." << endl;
            reader.Close();
            exit(1);
        }
    }
}