New upstream version 1.3.1

#CFLAGS=-O3 -D_FILE_OFFSET_BITS=64 -g

#CXXFLAGS=-O3 -D_FILE_OFFSET_BITS=64 -g

all: vcflib/Makefile log

	cd src && $(MAKE)

## user manual and guide

[![Build Status](https://travis-ci.org/ekg/freebayes.svg)](https://travis-ci.org/ekg/freebayes)

[![Gitter](https://badges.gitter.im/Join Chat.svg)](https://gitter.im/ekg/freebayes?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)

--------

variant detection--- that identical sequences may have multiple possible 

alignments:

<img src="http://hypervolu.me/freebayes/figures/haplotype_calling.png" 

width=500/>

<img src="https://github.com/ekg/freebayes/raw/v1.3.0/paper/haplotype_calling.png" width=500/>

*FreeBayes* uses short-read alignments 

([BAM](http://samtools.sourceforge.net/SAMv1.pdf) files with 

    freebayes --help

## Getting the best results

freebayes provides many options to configure its operation.

These may be important in certain use cases, but they are not meant for standard use.

It is strongly recommended that you run freebayes with parameters that are as close to default as possible unless you can validate that the chosen parameters improve performance.

To achieve the desired tradeoff between sensitivity and specificity, the best approach is to filter the output using the `QUAL` field, which encodes the posterior estimate of the probability of variation.

Filtering the input with the provided options demonstrably hurts performance where truth sets can be used to evaluate results.

By removing information from the input, you can confuse the Bayesian model by making it appear that certain alleles frequently indicative of context specific error (such as indels in homopolymers) don't exist.

Many users apply these filters to force freebayes to not make haplotype calls.

This is almost always a mistake, as the haplotype calling process greatly improves the method's signal to noise ratio and normalizes differential alignment in complex regions.

If you wish to obtain a VCF that does not contain haplotype calls or complex alleles, first call with default parameters and then decompose the output with tools in vcflib, vt, vcf-tools, bcftools, GATK, and Picard.

Here we use a tool in vcflib that normalizes the haplotype calls into pointwise SNPs and indels:

    freebayes ... | vcfallelicprimitives -kg >calls.vcf

Note that this is not done by default as it makes it difficult to determine which variant calls freebayes completed.

The raw output faithfully describes exactly the calls that were made.

## Examples

    freebayes -f ref.fa -C 5 aln.bam >var.vcf

Skip over regions of high depth by discarding alignments overlapping positions where total read depth is greater than 200:

    freebayes -f ref.fa -g 200 aln.bam >var.vcf

Use a different ploidy:

    freebayes -f ref.fa -p 4 aln.bam >var.vcf

        // a dangerous game

        int start = pos - currentSequenceStart;

        double minEntropy = parameters.minRepeatEntropy;

        // check first that' wer'e actually ina repeat... TODO

        //cerr << "entropy of " << entropy(currentSequence.substr(start, repeatRightBoundary - pos)) << " is too low, " << endl;

        while (minEntropy > 0 && // ignore if turned off

               repeatRightBoundary - currentSequenceStart < currentSequence.size() && //guard

               // don't run off the end of the current sequence

               repeatRightBoundary - currentSequenceStart < currentSequence.size() &&

               // there is no point in going past the alignment end

               // because we won't make a haplotype call unless we have a covering observation from a read

               repeatRightBoundary < alignment.ENDPOSITION &&

               entropy(currentSequence.substr(start, repeatRightBoundary - pos)) < minEntropy) {

            //cerr << "entropy of " << entropy(currentSequence.substr(start, repeatRightBoundary - pos)) << " is too low, ";

            //cerr << "increasing rought boundary to ";

            ++repeatRightBoundary;

            //cerr << repeatRightBoundary << endl;

        }

        // now we

        //cachedRepeatCounts[pos] = repeatCounts(pos - currentSequenceStart, currentSequence, 12);

        // edge case, the indel is an insertion and matches the reference to the right

        // this means there is a repeat structure in the read, but not the ref

        if (currentSequence.substr(pos - currentSequenceStart, length) == readSequence) {

                if (parameters.baseQualityCap != 0) {

                    capBaseQuality(currentAlignment, parameters.baseQualityCap);

                }

                // do we exceed coverage anywhere?

                // do we touch anything where we had exceeded coverage?

                // if so skip this read, and mark and remove processed alignments and registered alleles overlapping the coverage capped position

                bool considerAlignment = true;

                if (parameters.skipCoverage > 0) {

                    for (unsigned long int i =  currentAlignment.POSITION; i < currentAlignment.ENDPOSITION; ++i) {

                        unsigned long int x = ++coverage[i];

                        if (x > parameters.skipCoverage && !gettingPartials) {

                            considerAlignment = false;

                            // we're exceeding coverage at this position for the first time, so clean up

                            if (!coverageSkippedPositions.count(i)) {

                                // clean up reads overlapping this position

                                removeCoverageSkippedAlleles(registeredAlleles, i);

                                removeCoverageSkippedAlleles(newAlleles, i);

                                // remove the alignments overlapping this position

                                removeRegisteredAlignmentsOverlappingPosition(i);

                                // record that the position is capped

                                coverageSkippedPositions.insert(i);

                            }

                        }

                    }

                }

                // decomposes alignment into a set of alleles

                // here we get the deque of alignments ending at this alignment's end position

                deque<RegisteredAlignment>& rq = registeredAlignments[currentAlignment.ENDPOSITION];

                //cerr << "parameters capcoverage " << parameters.capCoverage << " " << rq.size() << endl;

                if (considerAlignment) {

                    // and insert the registered alignment into that deque

                rq.push_front(RegisteredAlignment(currentAlignment, parameters));

                    rq.push_front(RegisteredAlignment(currentAlignment));

                    RegisteredAlignment& ra = rq.front();

                    registerAlignment(currentAlignment, ra, sampleName, sequencingTech);

                    // backtracking if we have too many mismatches

                        }

                    }

                }

            }

	    } while ((hasMoreAlignments = GETNEXT(bamMultiReader, currentAlignment))

                 && currentAlignment.POSITION <= position

                 && currentAlignment.REFID == currentRefID);

}

void AlleleParser::removeRegisteredAlignmentsOverlappingPosition(long unsigned int pos) {

    map<long unsigned int, deque<RegisteredAlignment> >::iterator f = registeredAlignments.begin();

    map<long unsigned int, set<deque<RegisteredAlignment>::iterator> > alignmentsToErase;

    set<Allele*> allelesToErase;

    while (f != registeredAlignments.end()) {

        for (deque<RegisteredAlignment>::iterator d = f->second.begin(); d != f->second.end(); ++d) {

            if (d->start <= pos && d->end > pos) {

                alignmentsToErase[f->first].insert(d);

                for (vector<Allele>::iterator a = d->alleles.begin(); a != d->alleles.end(); ++a) {

                    allelesToErase.insert(&*a);

                }

            }

        }

        ++f;

    }

    // clean up registered alleles--- maybe this should be done externally?

    for (vector<Allele*>::iterator a = registeredAlleles.begin(); a != registeredAlleles.end(); ++a) {

        if (allelesToErase.count(*a)) {

            *a = NULL;

        }

    }

    registeredAlleles.erase(remove(registeredAlleles.begin(), registeredAlleles.end(), (Allele*)NULL), registeredAlleles.end());

    if (alignmentsToErase.size()) {

        for (map<long unsigned int, set<deque<RegisteredAlignment>::iterator> >::iterator e = alignmentsToErase.begin();

             e != alignmentsToErase.end(); ++e) {

            deque<RegisteredAlignment> updated;

            map<long unsigned int, deque<RegisteredAlignment> >::iterator f = registeredAlignments.find(e->first);

            assert(f != registeredAlignments.end());

            for (deque<RegisteredAlignment>::iterator d = f->second.begin(); d != f->second.end(); ++d) {

                if (!e->second.count(d)) {

                    updated.push_back(*d);

                }

            }

            f->second = updated;

        }

    }

}

void AlleleParser::addToRegisteredAlleles(vector<Allele*>& alleles) {

    registeredAlleles.insert(registeredAlleles.end(),

                             alleles.begin(),

    alleles.erase(remove(alleles.begin(), alleles.end(), (Allele*)NULL), alleles.end());

}

void AlleleParser::removeCoverageSkippedAlleles(vector<Allele*>& alleles, long int position) {

    for (vector<Allele*>::iterator a = alleles.begin(); a != alleles.end(); ++a) {

        Allele* allele = *a;

        if (*a != NULL && allele->alignmentStart <= position && allele->alignmentEnd > position) {

            allele->processed = true;

            *a = NULL;

        }

    }

    alleles.erase(remove(alleles.begin(), alleles.end(), (Allele*)NULL), alleles.end());

}

// steps our position/beddata/reference pointers through all positions in all

// targets, returns false when we are finished

//

    DEBUG("to next target");

    clearRegisteredAlignments();

    coverageSkippedPositions.clear();

    cachedRepeatCounts.clear();

    coverage.clear();

    // reset haplotype length; there is no last call in this sequence; it isn't relevant

    lastHaplotypeLength = 0;

                || (registeredAlignments.empty() && currentRefID != currentAlignment.REFID)) {

                DEBUG("at end of sequence");

                clearRegisteredAlignments();

                coverageSkippedPositions.clear();

                cachedRepeatCounts.clear();

                coverage.clear();

                loadNextPositionWithAlignmentOrInputVariant(currentAlignment);

                justSwitchedTargets = true;

            }

    // done typically at each new read, but this handles the case where there is no data for a while

    //updateInputVariants(currentPosition, 1);

    //DEBUG2("updating registered alleles");

    //updateRegisteredAlleles(); // this removes unused left-flanking sequence

    //DEBUG2("updating prior variant alleles");

    //updatePriorAlleles();

    // remove past registered alleles

    DEBUG2("marking previous alleles as processed and removing from registered alleles");

    removePreviousAlleles(registeredAlleles, currentPosition);

        cachedRepeatCounts.erase(rc++);

    }

    DEBUG2("erasing old coverage cap");

    while (coverageSkippedPositions.size() && *coverageSkippedPositions.begin() < currentPosition) {

        coverageSkippedPositions.erase(coverageSkippedPositions.begin());

    }

    DEBUG2("erasing old coverage counts");

    map<long unsigned int, long unsigned int>::iterator cov = coverage.begin();

    while (cov != coverage.end() && cov->first < currentPosition) {

        coverage.erase(cov++);

    }

    return true;

}

void AlleleParser::getAlleles(Samples& samples, int allowedAlleleTypes,

                              int haplotypeLength, bool getAllAllelesInHaplotype,

                              bool ignoreProcessedFlag) {

    Samples gvcf_held; // make some samples that by bass filtering for gvcf lines

    DEBUG2("getting alleles");

    for (Samples::iterator s = samples.begin(); s != samples.end(); ++s)

        s->second.clear();

    samples.clear();

    // Commenting this out and replacinf with .clear() to relly empty it, it is more aloc, but no major change

    //for (Samples::iterator s = samples.begin(); s != samples.end(); ++s)

    //    s->second.clear();

    // TODO ^^^ this should be optimized for better scanning performance

    // if we have targets and are outside of the current target, don't return anything

                  (allele.position == currentPosition)))

                ) ) {

            allele.update(haplotypeLength);

            if(parameters.gVCFout){

                gvcf_held[allele.sampleID][allele.currentBase].push_back(*a); // store things incase

            }

            if (allele.quality >= parameters.BQL0 && allele.currentBase != "N"

                && (allele.isReference() || !allele.alternateSequence.empty())) { // filters haplotype construction chaff

                //cerr << "keeping allele " << allele << endl;

            }

        }

    }

    if(samples.size() == 0 && parameters.gVCFout){

        samples = gvcf_held;  // if there are no non reference vals try to recover any allined values for gvcf if doing gvcf output!!

    }

    vector<string> samplesToErase;

    // now remove empty alleles from our return so as to not confuse processing

    int snpCount;

    int indelCount;

    int alleleTypes;

    Parameters parameters;

    RegisteredAlignment(BAMALIGN& alignment, Parameters parameters)

        //: alignment(alignment)

    RegisteredAlignment(BAMALIGN& alignment)

        : start(alignment.POSITION)

        , end(alignment.ENDPOSITION)

        , refid(alignment.REFID)

        , snpCount(0)

        , indelCount(0)

        , alleleTypes(0)

        , parameters(parameters)

    {

      FILLREADGROUP(readgroup, alignment);

    }

    vector<Allele*> registeredAlleles;

    map<long unsigned int, deque<RegisteredAlignment> > registeredAlignments;

    set<long unsigned int> coverageSkippedPositions;

    map<long unsigned int, long unsigned int> coverage;

    map<int, map<long int, vector<Allele> > > inputVariantAlleles; // all variants present in the input VCF, as 'genotype' alleles

    pair<int, long int> nextInputVariantPosition(void);

    void getInputVariantsInRegion(string& seq, long start = 0, long end = 0);

                                     int haplotypeLength = 1,

                                     bool getAllAllelesInHaplotype = false);

    void removePreviousAlleles(vector<Allele*>& alleles, long int position);

    void removeCoverageSkippedAlleles(vector<Allele*>& alleles, long int position);

    void removeRegisteredAlignmentsOverlappingPosition(long unsigned int pos);

    void removeFilteredAlleles(vector<Allele*>& alleles);

    void removeDuplicateAlleles(Samples& samples, map<string, vector<Allele*> >& alleleGroups,

                                int allowedAlleleTypes, int haplotypeLength, Allele& refallele);

################################################################################

# Compiler

CXXFLAGS=-O3 -D_FILE_OFFSET_BITS=64 -g -ggdb

CXX=g++ ${CXXFLAGS}

C=gcc

# Compiler flags

CFLAGS=-O3 -D_FILE_OFFSET_BITS=64 -g

#CFLAGS=-O3 -static -D VERBOSE_DEBUG  # enables verbose debugging via --debug2

SEQLIB_ROOT=../SeqLib

all: autoversion ../bin/freebayes ../bin/bamleftalign

static:

	$(MAKE) CXXFLAGS="$(CXXFLAGS) -static" all

	$(MAKE) CXXFLAGS="$(CXXFLAGS) -static -static-libstdc++ -static-libgcc -Wl,--allow-multiple-definition" all

debug:

	$(MAKE) CXXFLAGS="$(CXXFLAGS) -D VERBOSE_DEBUG -g -rdynamic" all