<commit_before>/*****************************************************************************
  multiBamCov.cpp

  (c) 2009 - <NAME>
  Hall Laboratory
  Department of Biochemistry and Molecular Genetics
  University of Virginia
  <EMAIL>

  Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#include "lineFileUtilities.h"
#include "multiBamCov.h"
#include "api/BamMultiReader.h"


/*
    Constructor
*/
MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly) 
:
_bam_files(bam_files),
_bed_file(bed_file),
_minQual(minQual),
_properOnly(properOnly)
{
	_bed = new BedFile(_bed_file);
    LoadBamFileMap();
}


/*
    Destructor
*/
MultiCovBam::~MultiCovBam(void) 
{}



void MultiCovBam::CollectCoverage()
{
    BamMultiReader reader;
    reader.SetIndexCacheMode(BamIndex::NoIndexCaching);
    
    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)
                {
                    reader.SetRegion(reader.GetReferenceID(bed.chrom),
                                           (int) bed.start,
                                           reader.GetReferenceID(bed.chrom),
                                           (int) bed.end);
                                           
                    // everything checks out, just iterate through specified region, counting alignments
                    vector<int> counts(_bam_files.size());
                    BamAlignment al;
                    while ( reader.GetNextAlignment(al) )
                    {
                        // map qual must exceed minimum
                        if (al.MapQuality >= _minQual) {
                            // ignore if not properly paired and we actually care.
                            if (_properOnly && !al.IsProperPair())
                                continue;

                            // 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);
        }
    }
}


void MultiCovBam::LoadBamFileMap(void) 
{
    for (size_t i = 0; i < _bam_files.size(); ++i)
    {
        bamFileMap[_bam_files[i]] = i;
    }
}

void MultiCovBam::ReportCounts(const vector<int> &counts) <commit_msg>New region checks for multiBamCov.  Thanks to <NAME>.<commit_after>/*****************************************************************************
  multiBamCov.cpp

  (c) 2009 - <NAME>
  Hall Laboratory
  Department of Biochemistry and Molecular Genetics
  University of Virginia
  <EMAIL>

  Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#include "lineFileUtilities.h"
#include "multiBamCov.h"
#include "api/BamMultiReader.h"


/*
    Constructor
*/
MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly) 
:
_bam_files(bam_files),
_bed_file(bed_file),
_minQual(minQual),
_properOnly(properOnly)
{
	_bed = new BedFile(_bed_file);
    LoadBamFileMap();
}


/*
    Destructor
*/
MultiCovBam::~MultiCovBam(void) 
{}



void MultiCovBam::CollectCoverage()
{
    BamMultiReader reader;
    reader.SetIndexCacheMode(BamIndex::NoIndexCaching);
    
    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)
                {
                    // initialize counts for each file to 0
                    vector<int> counts(_bam_files.size(), 0);
                    // get the BAM refId for this chrom.
                    int refId = reader.GetReferenceID(bed.chrom);
                    // set up a BamRegion to which to attempt to jump
                    BamRegion region(refId, (int)bed.start, refId, (int)bed.end);
                    
                    // everything checks out, just iterate through specified region, counting alignments
                    if ( (refId != -1) && (reader.SetRegion(region)) ) {
                        BamAlignment al;
                        while ( reader.GetNextAlignment(al) )
                        {
                            // map qual must exceed minimum
                            if (al.MapQuality >= _minQual) {
                                // ignore if not properly paired and we actually care.
                                if (_properOnly && !al.IsProperPair())
                                    continue;

                                // 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);
        }
    }
}


void MultiCovBam::LoadBamFileMap(void) 
{
    for (size_t i = 0; i < _bam_files.size(); ++i)
    {
        bamFileMap[_bam_files[i]] = i;
    }
}

void MultiCovBam::ReportCounts(const vector<int> &counts) 
{'fraction_non_alphanumeric': 0.13166454532597935, 'fraction_numerical': 0.002969876962240136, 'mean_word_length': 1.5475504322766571, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '340564', 'n_tokens_mistral': 1964, 'n_tokens_neox': 1827, 'n_words': 476}