<commit_before>
#include <stdio.h>
#include <string>
#include <iostream>
#include <fstream>
#include <vector>
#include "api/BamReader.h"
#include "transrate-pileup.h"

using namespace BamTools;

using namespace std;

double nullprior = 0.5;

struct ContigRecord {
  int bases_mapped;
  int p_seq_true;
  int bridges;
  int length;
  string name;
  int fragments_mapped;
  int both_mapped;
  int properpair;
  int good;
  int bases_uncovered;
  double p_unique;
  double p_not_segmented;
};

class BetterBam {
    int realistic_distance;
    int seq_count;
    int i,j;
    uint32_t nm_tag;
    int ldist;
    int rdist;
    int maxL=0;
    std::string file;
    BamReader reader;
    BamAlignment alignment;
  public:
    std::vector<ContigRecord> array;
    BetterBam (std::string);

    // set realistic distance between pairs to have 0.03% false positive rate
    void set_fragment_size(int size, int sd) {
      realistic_distance = size + 3 * sd;
    }

    int load_bam() {
      if (!reader.Open(file)) {
        cerr << "Could not open BAM file" << endl;
        return 1;
      }
      // get sam header
      SamSequenceDictionary dictionary = reader.GetHeader().Sequences;
      seq_count = dictionary.Size();
      array.resize(seq_count);
      // get an iterator for looping over the sequences in the header
      std::vector<SamSequence>::iterator it = dictionary.Begin();
      // fill the vector with intial values
      for (i = 0; i < seq_count; i++) {
        array[i].bases_mapped = 0;
        array[i].p_seq_true = 0;
        array[i].bridges = 0;
        if (it[i].HasLength()) {
          array[i].length = atoi(it[i].Length.c_str());
          if (array[i].length > maxL) {
            maxL = array[i].length;
          }
        } else {
          array[i].length = 0;
        }
        array[i].name = it[i].Name;
        array[i].fragments_mapped = 0;
        array[i].both_mapped = 0;
        array[i].properpair = 0;
        array[i].good = 0;
        array[i].bases_uncovered = 0;
        array[i].p_unique = 0;
        array[i].p_not_segmented = 1;
      }
      // loop through bam file
      i = -2;
      TransratePileup pileup(maxL);
      int ref_length = -1;
      while (reader.GetNextAlignment(alignment)) {

        // read must be mapped
        if (!alignment.IsMapped()) {
          continue;
        }
        // check this read comes from the currently loaded contig
        // if not, load the new contig
        if (alignment.RefID != i) {
          if (i>=0) {
            array[i].bases_uncovered = pileup.getBasesUncovered();
            array[i].p_unique = pileup.getUniqueBases();
            array[i].p_not_segmented = pileup.p_not_segmented();
          }
          i = alignment.RefID;
          ref_length = array[i].length;
          pileup.clearCoverage(ref_length);
        }

        pileup.addAlignment(alignment);

        array[i].bases_mapped += alignment.Length;

        // store edit distance for sequence accuracy calculation
        if (alignment.HasTag("NM")) {
          if (alignment.GetTag("NM", nm_tag)) {
            array[i].p_seq_true += nm_tag;
          }
        }

        // count fragments where either or both mates mapped
        if (alignment.IsFirstMate() ||
          (alignment.IsSecondMate() && !alignment.IsMateMapped())) {
          array[i].fragments_mapped++;
        }

        // from now on ignore fragments unless both mates mapped
        if (!(alignment.IsFirstMate() && alignment.IsMateMapped())) {
          continue;
        }

        array[i].both_mapped++;

        // mates must align to same contig, otherwise we record a bridge
        if (alignment.RefID != alignment.MateRefID) {
          array[i].bridges++;
          continue;
        }

        // fragment length must be plausible
        ldist = max(alignment.Position-alignment.MatePosition,
                    alignment.MatePosition-alignment.Position);
        if (ldist > realistic_distance) {
          // mates are too far apart
          continue;
        }

        // read orientation must match the generated library
        // in this case we only test for FR/RF orientation,
        // that is - we expect mates to be on opposite strands
        bool is_reversed = alignment.IsReverseStrand();
        bool is_mate_reversed = alignment.IsMateReverseStrand();<commit_msg>still count proper pairs<commit_after>
#include <stdio.h>
#include <string>
#include <iostream>
#include <fstream>
#include <vector>
#include "api/BamReader.h"
#include "transrate-pileup.h"

using namespace BamTools;

using namespace std;

double nullprior = 0.5;

struct ContigRecord {
  int bases_mapped;
  int p_seq_true;
  int bridges;
  int length;
  string name;
  int fragments_mapped;
  int both_mapped;
  int properpair;
  int good;
  int bases_uncovered;
  double p_unique;
  double p_not_segmented;
};

class BetterBam {
    int realistic_distance;
    int seq_count;
    int i,j;
    uint32_t nm_tag;
    int ldist;
    int rdist;
    int maxL=0;
    std::string file;
    BamReader reader;
    BamAlignment alignment;
  public:
    std::vector<ContigRecord> array;
    BetterBam (std::string);

    // set realistic distance between pairs to have 0.03% false positive rate
    void set_fragment_size(int size, int sd) {
      realistic_distance = size + 3 * sd;
    }

    int load_bam() {
      if (!reader.Open(file)) {
        cerr << "Could not open BAM file" << endl;
        return 1;
      }
      // get sam header
      SamSequenceDictionary dictionary = reader.GetHeader().Sequences;
      seq_count = dictionary.Size();
      array.resize(seq_count);
      // get an iterator for looping over the sequences in the header
      std::vector<SamSequence>::iterator it = dictionary.Begin();
      // fill the vector with intial values
      for (i = 0; i < seq_count; i++) {
        array[i].bases_mapped = 0;
        array[i].p_seq_true = 0;
        array[i].bridges = 0;
        if (it[i].HasLength()) {
          array[i].length = atoi(it[i].Length.c_str());
          if (array[i].length > maxL) {
            maxL = array[i].length;
          }
        } else {
          array[i].length = 0;
        }
        array[i].name = it[i].Name;
        array[i].fragments_mapped = 0;
        array[i].both_mapped = 0;
        array[i].properpair = 0;
        array[i].good = 0;
        array[i].bases_uncovered = 0;
        array[i].p_unique = 0;
        array[i].p_not_segmented = 1;
      }
      // loop through bam file
      i = -2;
      TransratePileup pileup(maxL);
      int ref_length = -1;
      while (reader.GetNextAlignment(alignment)) {

        // read must be mapped
        if (!alignment.IsMapped()) {
          continue;
        }

        // check this read comes from the currently loaded contig
        // if not, load the new contig
        if (alignment.RefID != i) {
          if (i>=0) {
            array[i].bases_uncovered = pileup.getBasesUncovered();
            array[i].p_unique = pileup.getUniqueBases();
            array[i].p_not_segmented = pileup.p_not_segmented();
          }
          i = alignment.RefID;
          ref_length = array[i].length;
          pileup.clearCoverage(ref_length);
        }

        pileup.addAlignment(alignment);

        array[i].bases_mapped += alignment.Length;

        // store edit distance for sequence accuracy calculation
        if (alignment.HasTag("NM")) {
          if (alignment.GetTag("NM", nm_tag)) {
            array[i].p_seq_true += nm_tag;
          }
        }

        // count fragments where either or both mates mapped
        if (alignment.IsFirstMate() ||
          (alignment.IsSecondMate() && !alignment.IsMateMapped())) {
          array[i].fragments_mapped++;
        }

        // from now on ignore fragments unless both mates mapped
        if (!(alignment.IsFirstMate() && alignment.IsMateMapped())) {
          continue;
        }

        array[i].both_mapped++;

        // count proper pairs, although we have our own definition because
        // not all aligners use the same definition
        if (alignment.IsProperPair()) {
          ++array[i].properpair;
        }

        // mates must align to same contig, otherwise we record a bridge
        if (alignment.RefID != alignment.MateRefID) {
          array[i].bridges++;
          continue;
        }

        // fragment length must be plausible
        ldist = max(alignment.Position-alignment.MatePosition,
                    alignment.MatePosition-alignment.Position);
        if (ldist > realistic_distance) {
          // mates are too far apart
          continue;
        }

        // read orientation must match the generated library
        // in this case we only test for FR/RF orientation,
        // that is - we expect mates to be on opposite strands
        bool is_reversed = alignment.IsReverseStrand();
        bool is_mate_reversed = alignment.IsMateReverseStrand();
{'fraction_non_alphanumeric': 0.10935574229691877, 'fraction_numerical': 0.0056022408963585435, 'mean_word_length': 2.110104529616725, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1966281', 'n_tokens_mistral': 2811, 'n_tokens_neox': 2669, 'n_words': 700}