<commit_before>#include "ResultData.h"
#include "TryCatch.h"

using namespace std;



vcf::Variant& Results::vcf(
        vcf::Variant& var, // variant to update
        long double pHom,
        long double bestComboOddsRatio,
        //long double alleleSamplingProb,
        Samples& samples,
        string refbase,
        vector<Allele>& altAllelesIncludingNulls,
        map<string, int> repeats,
        vector<string>& sampleNames,
        int coverage,
        GenotypeCombo& genotypeCombo,
        bool bestOverallComboIsHet,
        map<string, vector<Allele*> >& alleleGroups,
        map<int, vector<Genotype> >& genotypesByPloidy,
        vector<string>& sequencingTechnologies,
        AlleleParser* parser) {

    Parameters& parameters = parser->parameters;

    GenotypeComboMap comboMap;
    genotypeCombo2Map(genotypeCombo, comboMap);

    // set up the reported reference allele
    long int referencePosition = (long int) parser->currentPosition; // 0-based

    // remove alt alleles
    vector<Allele> altAlleles;
    for (vector<Allele>::iterator aa = altAllelesIncludingNulls.begin(); aa != altAllelesIncludingNulls.end(); ++aa) {
        if (!aa->isNull()) {
            altAlleles.push_back(*aa);
        }
    }

    // adjust reference position, reference sequence, and alt alleles by
    // stripping invariant bases off the beginning and end of the alt alleles
    // first we find the minimum start and end matches
    vector<Allele> adjustedAltAlleles;
    int minStartMatch = 0;
    int minEndMatch = 0;
    for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
        vector<pair<int, string> > cigar = splitCigar(aa->cigar);
        int startmatch = 0;
        int endmatch = 0;
        if (cigar.front().second == "M") {
            startmatch = cigar.front().first;
        }
        if (cigar.back().second == "M") {
            endmatch = cigar.back().first;
        }
	// check excludes complex alleles of the form, e.g. 1X3I
        if (cigar.size() > 1 && cigar.front().second == "M" && (cigar.at(1).second == "D" || cigar.at(1).second == "I")) {
            startmatch -= 1; // require at least one base flanking deletions
        }
        if (aa == altAlleles.begin()) {
            minStartMatch = startmatch;
            minEndMatch = endmatch;
        } else {
            minStartMatch = min(minStartMatch, startmatch);
            minEndMatch = min(minEndMatch, endmatch);
        }
    }
    // if either is non-zero, we have to adjust cigars and alternate sequences to be printed
    // this is done solely for reporting, so the altAlleles structure is used
    // for stats generation out of the ML genotype combination
    map<string, string> adjustedCigar;
    if (minStartMatch || minEndMatch) {
        for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
            // subtract the minStartMatch and minEndMatch bases from the allele start and end
            adjustedAltAlleles.push_back(*aa);
            Allele& allele = adjustedAltAlleles.back();
            vector<pair<int, string> > cigar = splitCigar(allele.cigar);
            // TODO clean this up by writing a wrapper for Allele::subtract() (?)
            if (cigar.front().second == "M") {
                cigar.front().first -= minStartMatch;
                allele.alternateSequence = allele.alternateSequence.substr(minStartMatch);
            }
            if (cigar.back().second == "M") {
                cigar.back().first -= minEndMatch;
                allele.alternateSequence = allele.alternateSequence.substr(0, allele.alternateSequence.size() - minEndMatch);
            }
            allele.cigar = joinCigar(cigar);
            allele.position += minStartMatch;
            allele.referenceLength -= minStartMatch + minEndMatch;
            adjustedCigar[aa->base()] = allele.cigar;
        }
        referencePosition += minStartMatch;
    } else {
        adjustedAltAlleles = altAlleles;
        for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
            adjustedCigar[aa->base()] = aa->cigar;
        }
    }

    var.ref = parser->referenceSubstr(referencePosition, 1);

    // the reference sequence should be able to encompass all events at the site, +1bp on the left
    for (vector<Allele>::iterator aa = adjustedAltAlleles.begin(); aa != adjustedAltAlleles.end(); ++aa) {

        Allele& altAllele = *aa;
        switch (altAllele.type) {
            case ALLELE_SNP:
            case ALLELE_REFERENCE:
            case ALLELE_MNP:
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            case ALLELE_DELETION:
                // extend the reference sequence
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            case ALLELE_INSERTION:
                break;
            case ALLELE_COMPLEX:
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            default:
                cerr << "Unhandled allele type: " << altAllele.typeStr() << endl;
                break;
        }

    }

    for (vector<Allele>::iterator aa = adjustedAltAlleles.begin(); aa != adjustedAltAlleles.end(); ++aa) {
        Allele& altAllele = *aa;
        string altseq;
        switch (altAllele.type) {
            case ALLELE_REFERENCE:
                break;
            case ALLELE_SNP:
            case ALLELE_MNP:
                altseq = var.ref;
                altseq.replace(0, altAllele.alternateSequence.size(), altAllele.alternateSequence);
                var.alt.push_back(altseq);
                break;
            case ALLELE_DELETION:
            case ALLELE_INSERTION: // XXX is this correct???
            case ALLELE_COMPLEX:
                var.alt.push_back(altAllele.alternateSequence);
                break;
            default:
                cerr << "Unhandled allele type: " << altAllele.typeStr() << endl;
                break;
        }
    }

    assert(!var.ref.empty());
    for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
        assert(!a->empty());
        if (*a == var.ref) {
	    cerr << "variant at " << parser->currentSequenceName << ":" << referencePosition + 1 << endl;
	    cerr << "alt is the same as the reference" << endl;
	    cerr << *a << " == " << var.ref << endl;
	    assert(*a != var.ref);
	}
    }


    // get the required size of the reference sequence


    // set up VCF record-wide variables

    var.sequenceName = parser->currentSequenceName;
    // XXX this should be the position of the matching reference haplotype
    var.position = referencePosition + 1;
    var.id = ".";
    var.filter = ".";
    // XXX this should be the size of the maximum deletion + 1bp on the left end
    var.quality = float2phred(pHom);
<commit_msg>remove assertion (and thus exit) on alt == ref<commit_after>#include "ResultData.h"
#include "TryCatch.h"

using namespace std;



vcf::Variant& Results::vcf(
        vcf::Variant& var, // variant to update
        long double pHom,
        long double bestComboOddsRatio,
        //long double alleleSamplingProb,
        Samples& samples,
        string refbase,
        vector<Allele>& altAllelesIncludingNulls,
        map<string, int> repeats,
        vector<string>& sampleNames,
        int coverage,
        GenotypeCombo& genotypeCombo,
        bool bestOverallComboIsHet,
        map<string, vector<Allele*> >& alleleGroups,
        map<int, vector<Genotype> >& genotypesByPloidy,
        vector<string>& sequencingTechnologies,
        AlleleParser* parser) {

    Parameters& parameters = parser->parameters;

    GenotypeComboMap comboMap;
    genotypeCombo2Map(genotypeCombo, comboMap);

    // set up the reported reference allele
    long int referencePosition = (long int) parser->currentPosition; // 0-based

    // remove alt alleles
    vector<Allele> altAlleles;
    for (vector<Allele>::iterator aa = altAllelesIncludingNulls.begin(); aa != altAllelesIncludingNulls.end(); ++aa) {
        if (!aa->isNull()) {
            altAlleles.push_back(*aa);
        }
    }

    // adjust reference position, reference sequence, and alt alleles by
    // stripping invariant bases off the beginning and end of the alt alleles
    // first we find the minimum start and end matches
    vector<Allele> adjustedAltAlleles;
    int minStartMatch = 0;
    int minEndMatch = 0;
    for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
        vector<pair<int, string> > cigar = splitCigar(aa->cigar);
        int startmatch = 0;
        int endmatch = 0;
        if (cigar.front().second == "M") {
            startmatch = cigar.front().first;
        }
        if (cigar.back().second == "M") {
            endmatch = cigar.back().first;
        }
	// check excludes complex alleles of the form, e.g. 1X3I
        if (cigar.size() > 1 && cigar.front().second == "M" && (cigar.at(1).second == "D" || cigar.at(1).second == "I")) {
            startmatch -= 1; // require at least one base flanking deletions
        }
        if (aa == altAlleles.begin()) {
            minStartMatch = startmatch;
            minEndMatch = endmatch;
        } else {
            minStartMatch = min(minStartMatch, startmatch);
            minEndMatch = min(minEndMatch, endmatch);
        }
    }
    // if either is non-zero, we have to adjust cigars and alternate sequences to be printed
    // this is done solely for reporting, so the altAlleles structure is used
    // for stats generation out of the ML genotype combination
    map<string, string> adjustedCigar;
    if (minStartMatch || minEndMatch) {
        for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
            // subtract the minStartMatch and minEndMatch bases from the allele start and end
            adjustedAltAlleles.push_back(*aa);
            Allele& allele = adjustedAltAlleles.back();
            vector<pair<int, string> > cigar = splitCigar(allele.cigar);
            // TODO clean this up by writing a wrapper for Allele::subtract() (?)
            if (cigar.front().second == "M") {
                cigar.front().first -= minStartMatch;
                allele.alternateSequence = allele.alternateSequence.substr(minStartMatch);
            }
            if (cigar.back().second == "M") {
                cigar.back().first -= minEndMatch;
                allele.alternateSequence = allele.alternateSequence.substr(0, allele.alternateSequence.size() - minEndMatch);
            }
            allele.cigar = joinCigar(cigar);
            allele.position += minStartMatch;
            allele.referenceLength -= minStartMatch + minEndMatch;
            adjustedCigar[aa->base()] = allele.cigar;
        }
        referencePosition += minStartMatch;
    } else {
        adjustedAltAlleles = altAlleles;
        for (vector<Allele>::iterator aa = altAlleles.begin(); aa != altAlleles.end(); ++aa) {
            adjustedCigar[aa->base()] = aa->cigar;
        }
    }

    var.ref = parser->referenceSubstr(referencePosition, 1);

    // the reference sequence should be able to encompass all events at the site, +1bp on the left
    for (vector<Allele>::iterator aa = adjustedAltAlleles.begin(); aa != adjustedAltAlleles.end(); ++aa) {

        Allele& altAllele = *aa;
        switch (altAllele.type) {
            case ALLELE_SNP:
            case ALLELE_REFERENCE:
            case ALLELE_MNP:
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            case ALLELE_DELETION:
                // extend the reference sequence
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            case ALLELE_INSERTION:
                break;
            case ALLELE_COMPLEX:
                if (var.ref.size() < altAllele.referenceLength) {
                    var.ref = parser->referenceSubstr(referencePosition, altAllele.referenceLength);
                }
                break;
            default:
                cerr << "Unhandled allele type: " << altAllele.typeStr() << endl;
                break;
        }

    }

    for (vector<Allele>::iterator aa = adjustedAltAlleles.begin(); aa != adjustedAltAlleles.end(); ++aa) {
        Allele& altAllele = *aa;
        string altseq;
        switch (altAllele.type) {
            case ALLELE_REFERENCE:
                break;
            case ALLELE_SNP:
            case ALLELE_MNP:
                altseq = var.ref;
                altseq.replace(0, altAllele.alternateSequence.size(), altAllele.alternateSequence);
                var.alt.push_back(altseq);
                break;
            case ALLELE_DELETION:
            case ALLELE_INSERTION: // XXX is this correct???
            case ALLELE_COMPLEX:
                var.alt.push_back(altAllele.alternateSequence);
                break;
            default:
                cerr << "Unhandled allele type: " << altAllele.typeStr() << endl;
                break;
        }
    }

    assert(!var.ref.empty());
    for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
        assert(!a->empty());
        if (*a == var.ref) {
	    cerr << "variant at " << parser->currentSequenceName << ":" << referencePosition + 1 << endl;
	    cerr << "alt is the same as the reference" << endl;
	    cerr << *a << " == " << var.ref << endl;
	}
    }


    // get the required size of the reference sequence


    // set up VCF record-wide variables

    var.sequenceName = parser->currentSequenceName;
    // XXX this should be the position of the matching reference haplotype
    var.position = referencePosition + 1;
    var.id = ".";
    var.filter = ".";
    // XXX this should be the size of the maximum deletion + 1bp on the left end
    var.quality = float2phred(pHom);

{'fraction_non_alphanumeric': 0.11365359432076642, 'fraction_numerical': 0.0027569095044455167, 'mean_word_length': 2.1030795551753636, 'pattern_counts': {'":': 2, '<': 121, '<?xml version=': 0, '>': 77, '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': '2180075', 'n_tokens_mistral': 4393, 'n_tokens_neox': 4132, 'n_words': 1039}