def get_name(seq):
    if 'SeqRecord' in seq.__class__.__name__:
        seq_class = SEQRECORD
    else:
        seq_class = seq.kind
        seq = seq.object
    if seq_class == SEQITEM:
        name = seq.name
    elif seq_class == SEQRECORD:
        name = seq.id
    return name def get_length(seq):
    return len(get_str_seq(seq)) def read_seqs(fhands, out_format=None, prefered_seq_classes=None):
    'It returns a stream of seqs in different codings: seqrecords, seqitems...'

    if not prefered_seq_classes:
        prefered_seq_classes = [SEQITEM, SEQRECORD]
    try:
        in_format = get_format(fhands[0])
    except FileIsEmptyError:
        return []
    # seqitems is incompatible with different input and output formats
    # or when in_format != a fasta or fastq
    if ((out_format not in (None, GUESS_FORMAT) and in_format != out_format
         and SEQITEM in prefered_seq_classes) or
        (in_format not in ('fasta',) + SANGER_FASTQ_FORMATS +
         ILLUMINA_FASTQ_FORMATS)):
        prefered_seq_classes.pop(prefered_seq_classes.index(SEQITEM))

    if not prefered_seq_classes:
        msg = 'No valid seq class left or prefered'
        raise ValueError(msg)

    for seq_class in prefered_seq_classes:
        if seq_class == SEQITEM:
            try:
                return _read_seqitems(fhands)
            except NotImplementedError:
                continue
        elif seq_class == SEQRECORD:
            try:
                seqs = _read_seqrecords(fhands)
                return assing_kind_to_seqs(SEQRECORD, seqs, None)
            except NotImplementedError:
                continue
        else:
            raise ValueError('Unknown class for seq: ' + seq_class)
    raise RuntimeError('We should not be here, fixme') class IntCounter(Counter):
    '''This is a subclass  of the Counter on python collections.

    It adds some statistical functionalities to the class'''

    def __init__(self, *args, **kwargs):
        'It inialices a '
        super(IntCounter, self).__init__(*args, **kwargs)
        self.labels = LABELS.copy()

    @property
    def min(self):
        'Get the minimun value'
        return min(self.keys())

    @property
    def max(self):
        'Get the maximun value'
        return max(self.keys())

    @property
    def count(self):
        'It returns the count of the values stored in the array'
        return sum(self.values())

    @property
    def median(self):
        'It calculates the median of the values appended'
        quotient, remainder = divmod(self.count, 2)
        if remainder == 0:
            val1 = self._get_value_for_index(quotient - 1)
            val2 = self._get_value_for_index(quotient)
            return (val1 + val2) / 2
        else:
            return self._get_value_for_index(quotient)

    @property
    def sum(self):
        'It gets the sum of the values'
        sum_ = 0
        for index, value in self.items():
            sum_ += (index * value)
        return int(sum_)

    @property
    def average(self):
        'It calculates the average'
        count = self.count
        sum_ = self.sum
        return sum_ / count

    @property
    def variance(self):
        'It gets the variance of the values'
        mean = self.average
        sum_ = 0
        for index, counts in self.items():
            sum_ += ((index - mean) ** 2) * counts
        return sum_ / self.count

    @property
    def quartiles(self):
        'It returns the quartiles'
        num_items = self.count
        if num_items < 4:
            msg = 'At least 4 values are required to calculate the quartiles'
            raise RuntimeError(msg)
        # quartile 1
        quotient, remainder = divmod(num_items + 1, 4)
        if not remainder:
            quartile1 = self._get_value_for_index(quotient - 1)
        else:
            val1 = self._get_value_for_index(quotient - 1)
            val2 = self._get_value_for_index(quotient)
            quartile1 = (val1 + val2) / 2
        # quartile 3
        quotient, remainder = divmod((num_items + 1) * 3, 4)
        if not remainder:
            quartile3 = self._get_value_for_index(quotient - 1)
        else:
            val1 = self._get_value_for_index(quotient - 1)
            val2 = self._get_value_for_index(quotient)
            quartile3 = (val1 + val2) / 2
        return quartile1, self.median, quartile3

    @property
    def irq(self):
        'It gets the interquartile range'
        # pylint: disable=W0612
        quart1 = self.quartiles[0]
        quart3 = self.quartiles[2]
        return quart3 - quart1

    @property
    def outlier_limits(self):
        'It returns the intercuartile'
        # pylint: disable=W0612
        quart1 = self.quartiles[0]
        quart3 = self.quartiles[2]

        iqr = self.irq
        limit_distance = round(iqr * 1.5)

        start = int(quart1 - limit_distance)
        end = int(quart3 + limit_distance)
        return (start, end)

    def _get_value_for_index(self, position):
        '''It takes a position and it returns the value for the given index'''
        cum_count = 0
        for index in sorted(self.keys()):
            count = self[index]
            cum_count += count
            if position <= cum_count - 1:
                return index
        else:
            if position >= cum_count:
                raise IndexError('You asked for an index beyond the scope')
            return index

    def _calculate_dist_range(self, min_, max_, outlier_threshold):
        'it calculates the range for the histogram'
        if ((min_ is not None or max_ is not None) and
            outlier_threshold is not None):
            msg = 'You can not pass max, min and outlier_threslhosld to '
            msg += 'calculate distribution range'
            raise ValueError(msg)

        if min_ is None:
            min_ = self.min
        if max_ is None:
            max_ = self.max
        if outlier_threshold:
            left_limit = self.count * outlier_threshold / 100
            rigth_limit = self.count - left_limit
            left_value = self._get_value_for_index(left_limit)
            rigth_value = self._get_value_for_index(rigth_limit)

            if min_ < left_value:
                min_ = left_value
            if max_ > rigth_value:
                max_ = rigth_value
        return min_, max_

    def calculate_bin_edges(self, min_, max_, n_bins=None):
        'It calculates the bin_edges'
        min_bins = get_setting('MIN_BINS')
        max_bins = get_setting('MAX_BINS')
        if n_bins is None:
            num_values = int(max_ - min_)
            if num_values == 0:
                n_bins = 1
            elif num_values < min_bins:
                n_bins = num_values
            else:
                n_bins = int(self.count / get_setting('MEAN_VALUES_IN_BIN'))
                if n_bins < min_bins:
                    n_bins = min_bins
                if n_bins > max_bins:
                    n_bins = max_bins
                if n_bins > num_values:
                    n_bins = num_values

        # now we can calculate the bin edges
        distrib_span = max_ - min_ if max_ != min_ else 1

        if distrib_span % n_bins:
            distrib_span = distrib_span + n_bins - (distrib_span % n_bins)
        bin_span = distrib_span // n_bins
        bin_edges = [min_ + bin_ * bin_span for bin_ in range(n_bins + 1)]
        return bin_edges

    def calculate_distribution(self, bins=None, min_=None, max_=None,
                               outlier_threshold=None):
        'It returns an histogram with the given range and bin'
        if self.count == 0:
            raise RuntimeError('No items in IntCounter')
        distrib = []
        min_, max_ = self._calculate_dist_range(min_, max_, outlier_threshold)
        if min_ is None or max_ is None:
            return None
        bin_edges = self.calculate_bin_edges(min_, max_, bins)
        for bin_index, left_edge in enumerate(bin_edges):
            try:
                rigth_edge = bin_edges[bin_index + 1]
            except IndexError:
                break
            sum_values = 0

            for index2 in sorted(self.keys()):
                value = self[index2]
                if index2 > rigth_edge:
                    break

                elif (left_edge <= index2 and index2 < rigth_edge or
                      left_edge <= index2 and index2 == max_):
                    sum_values += value

            distrib.append(sum_values)
        return {'counts': distrib, 'bin_limits': bin_edges}

    def update_labels(self, labels):
        'It prepares the labels for output files'
        self.labels.update(labels)

    def count_relative_to_value(self, value, comparison):
        'It counts the ints greater, equal, etc, relative to the given value.'
        return sum([counts for val, counts in self.items()
                    if comparison(val, value)])

    def __add__(self, other):
        'Add counts from two counters.'
        counter_python = super(IntCounter, self).__add__(other)
        return self.__class__(counter_python)

    def __str__(self):
        'It writes some basic stats of the values'
        if self.count != 0:
            labels = self.labels
            # now we write some basic stats
            format_num = lambda x: '{:,d}'.format(x) if isinstance(x, int) \
                                                                else '%.2f' % x
            text = '{}: {}\n'.format(labels['minimum'], format_num(self.min))
            text += '{}: {}\n'.format(labels['maximum'], format_num(self.max))
            text += '{}: {}\n'.format(labels['average'],
                                      format_num(self.average))

            if labels['variance'] is not None:
                text += '{}: {}\n'.format(labels['variance'],
                                          format_num(self.variance))
            if labels['sum'] is not None:
                text += '{}: {}\n'.format(labels['sum'],
                                          format_num(self.sum))
            if labels['items'] is not None:
                text += '{}: {}\n'.format(labels['items'], self.count)
            text += '\n'

            distrib = self.calculate_distribution()
            text += draw_histogram_ascii(distrib['bin_limits'],
                                         distrib['counts'])
            return text
        return '' class IntBoxplot(object):
    'It represents a Box and whisker plot'
    def __init__(self):
        'The init'
        self.counts = {}

    def append(self, category, value):
        'It appends a value to the distribution corresponding to a category'
        counts = self.counts
        try:
            cat_counts = counts[category]
        except KeyError:
            counts[category] = IntCounter()
            cat_counts = counts[category]
        cat_counts[value] += 1

    @property
    def aggregated_array(self):
        'It returns the IntSummarizedArray of all appended values.'
        aggregated_dist = None
        for cat_count in self.counts.viewvalues():
            if aggregated_dist is None:
                aggregated_dist = cat_count
            else:
                aggregated_dist += cat_count
        return aggregated_dist

    @property
    def ascii_plot(self):
        'It returns an string with an ASCII representation'
        distributions = self.counts
        categories = sorted(distributions.viewkeys())
        distrib_descriptions = {}
        for category in categories:
            distrib = distributions[category]
            min_ = distrib.min
            max_ = distrib.max
            try:
                quart1, median, quart3 = distrib.quartiles
                distrib_descriptions[category] = {'min': min_, 'max': max_,
                                            'quart1': quart1, 'median': median,
                                            'quart3': quart3}
            except RuntimeError:
                distrib_descriptions[category] = None

        min_value = min((d['min'] for d in distrib_descriptions.values()
                                                             if d is not None))
        max_value = max((d['max'] for d in distrib_descriptions.values()
                                                             if d is not None))

        len_str_float = lambda x: len('{:.1f}'.format(x))

        # how wide is the widest number of each type (min, max, median, etc)
        widths = {}
        for distrib in distrib_descriptions.values():
            if  distrib is None:
                continue
            for key in distrib.keys():
                value = len_str_float(distrib[key])
                if key not in widths or widths[key] < value:
                    widths[key] = value

        widths['labels'] = max([len(str(d)) for d in distrib_descriptions])

        category_format = '{:>' + str(widths['labels']) + 's}'
        header_format = category_format + ':{:>' + str(widths['min']) + '.1f}'
        header_format += ',{:>' + str(widths['quart1']) + '.1f}'
        header_format += ',{:>' + str(widths['median']) + '.1f}'
        header_format += ',{:>' + str(widths['quart3']) + '.1f}'
        header_format += ',{:>' + str(widths['max']) + '.1f}'

        header_len = None
        headers, distribs = [], []
        for category in categories:
            distrib = distrib_descriptions[category]
            if distrib is not None:
                header = header_format.format(str(category), distrib['min'],
                                          distrib['quart1'], distrib['median'],
                                          distrib['quart3'], distrib['max'])
                header += ' '
                if header_len is None or header_len < len(header):
                    header_len = len(header)
            else:
                header = None
            headers.append(header)
            distribs.append(distrib)

        plot_width = (get_setting('MAX_WIDTH_ASCII_PLOT') - widths['labels'] -
                      len(str(max_value)) - header_len)
        val_per_pixel = (max_value - min_value) / plot_width
        if val_per_pixel == 0:
            val_per_pixel = max_value / plot_width

        to_axis_scale = lambda x: int((x - min_value) / val_per_pixel)

        result = ''

        for distrib, header in zip(distribs, headers):
            if distrib is not None:
                line = header
                min_ = to_axis_scale(distrib['min'])
                line += ' ' * (min_ - 1)
                line += '<'  # min
                quart1 = to_axis_scale(distrib['quart1'])
                line += '-' * ((quart1 - min_) - 1)
                line += '['  # quartil 1
                median = to_axis_scale(distrib['median'])
                line += '=' * ((median - quart1) - 1)
                line += '|'  # median 1
                quart3 = to_axis_scale(distrib['quart3'])
                line += '=' * ((quart3 - median) - 1)
                line += ']'  # quart 2
                max_ = to_axis_scale(distrib['max'])
                line += '-' * ((max_ - quart3) - 1)
                line += '>'  # max
            else:
                line = category_format.format(str(category), 0, 0, 0, 0, 0)
            line += '\n'
            result += line

        axis1 = ' ' * header_len + '+' + '-' * (plot_width - 2) + '+' + '\n'
        axis2 = ' ' * header_len + str(min_value)
        axis2 += ' ' * (plot_width - len(str(min_value)) - len(str(max_value))
                       + 1)
        axis2 += str(max_value) + '\n'

        result += axis1
        result += axis2
        return result

    def __nonzero__(self):
        'It returns True if the object holds any counts'
        return bool(self.counts) def get_fig_and_canvas(num_rows=1, num_cols=1, figsize=None):
    if figsize is None:
        height = 5.0 * num_rows
        width = 7.5 * num_cols
        if height > 320.0:
            height = 320.0
        figsize = (width, height)
    try:
        fig = Figure(figsize=figsize)
        canvas = FigureCanvas(fig)
    except NameError:
        msg = 'Matplotlib module is required to draw graphical histograms'
        raise OptionalRequirementError(msg)
    return fig, canvas def draw_int_boxplot(boxplot, fhand=None, axes=None, title=None,
                     xlabel=None, ylabel=None):

    if axes is None:
        canvas, axes = get_canvas_and_axes()
        using_given_axes = False
    elif fhand is None:
        using_given_axes = True
    else:
        msg = 'If an axes is not given the fhand also required'
        raise NotImplementedError(msg)

    bar_width = 0.8

    x_vals = sorted(boxplot.counts.keys())

    # for the following line to work we would have to create all points in
    # memory, this is why we have reimplemented the boxplot
    # axes.boxplot([list(boxplot.counts[x_val].elements()) for x_val in x_vals]

    # we don't want more than 50 bars
    max_n_boxes = 40
    n_boxes = len(x_vals) - 1
    xpurge = n_boxes // max_n_boxes + 1

    xticks_lables = []
    xticks = []
    xpos = []
    for index, x_val in enumerate(x_vals):
        xpos.append(index + 0.5)
        if index % xpurge != 0:
            continue

        xticks_lables.append(str(x_val))
        # we add 0.5 to have the x ticks in  the
        # middle of the bar
        xticks.append(index + 0.5)

        intcounter = boxplot.counts[x_val]
        if intcounter.count < 5:
            # at least for values are required to calculate the qualties
            continue
        quart = intcounter.quartiles
        min_ = intcounter.min
        max_ = intcounter.max
        axes.bar(index + 0.5, quart[2] - quart[0], bottom=quart[0],
                 width=bar_width, facecolor='none', align='center')
        # median
        axes.plot([index + 0.5 - bar_width / 2, index + 0.5 + bar_width / 2],
                  [quart[1], quart[1]], color='red')
        # max
        axes.plot([index + 0.5 - bar_width / 4, index + 0.5 + bar_width / 4],
                  [max_, max_], color='black')
        axes.plot([index + 0.5, index + 0.5], [quart[0], min_], '--',
                  color='black')
        axes.plot([index + 0.5 - bar_width / 4, index + 0.5 + bar_width / 4],
                  [min_, min_], color='black')
        axes.plot([index + 0.5, index + 0.5], [quart[2], max_], '--',
                  color='black')

    max_x = len(x_vals) - 1
    axes.set_xlim(0, max_x + 1)

    axes.set_xticks(xticks)
    axes.set_xticklabels(xticks_lables)

    if title:
        axes.set_title(title)
    if xlabel:
        axes.set_xlabel(xlabel)
    if ylabel:
        axes.set_ylabel(ylabel)

    if not using_given_axes:
        canvas.print_figure(fhand)
        fhand.flush()
    return x_vals, xpos, xticks, xticks_lables VARSCAN = 'VarScan'
GATK = 'gatk'
FREEBAYES = 'freebayes'
GENERIC = 'generic'
HOM_REF = 0
HET = 1
HOM_ALT = 2
HOM = 3
MISSING_ALLELE_CHAR = '.'
DEF_MIN_CALLS_FOR_POP_STATS = 10
DEF_MIN_NUM_SNPS_IN_WIN = 5
class _SNPQueue(object):
class _SNPSlidingWindow(object):
class VCFReader(object):
class VCFWriter(pyvcfWriter):
class SNV(object):
class Call(object):
    def __init__(self):
    def empty(self):
    def pop(self, location):
    def extend(self, snps):
    def __init__(self, snp_reader, win_size, win_step, min_num_snps,
                 ref_fhand=None):
    def _snp_in_window(self, snp, win):
    def _get_chrom_lengths(self):
    def windows(self):
    def __init__(self, fhand, compressed=None, filename=None,
                 min_calls_for_pop_stats=DEF_MIN_CALLS_FOR_POP_STATS):
    def parse_snvs(self):
    def fetch_snvs(self, chrom, start, end=None):
    def sliding_windows(self, size, step=None, ref_fhand=None,
                        min_num_snps=DEF_MIN_NUM_SNPS_IN_WIN):
    def snpcaller(self):
    def samples(self):
    def filters(self):
    def infos(self):
    def _build_header(self, samples=None):
    def header(self):
    def create_template_reader(self, samples=None):
    def __init__(self, stream, template_reader, lineterminator="\n"):
    def write_snv(self, snv):
    def write_snvs(self, snvs):
    def flush(self):
def get_ids(snp):
def get_or_create_id(snp, prefix=''):
    def __init__(self, record, reader, ploidy=2,
                 min_calls_for_pop_stats=DEF_MIN_CALLS_FOR_POP_STATS):
    def ids(self):
    def get_or_create_id(self, prefix=''):
    def alleles(self):
    def obs_het(self):
    def exp_het(self):
    def calls(self):
    def get_call(self, sample):
    def _calculate_maf_and_mac(self):
    def allele_counts(self):
    def genotype_counts(self):
    def genotype_freqs(self):
    def biallelic_genotype_counts(self):
    def biallelic_genotype_freqs(self):
    def maf(self):
    def mac(self):
    def inbreed_coef(self):
    def _calculate_mafs_dp(self):
    def maf_depth(self):
    def allele_depths(self):
    def depth(self):
    def chrom(self):
    def pos(self):
    def end(self):
    def qual(self):
    def ref(self):
    def num_called(self):
    def call_rate(self):
    def is_snp(self):
    def __str__(self):
    def __unicode__(self):
    def __repr__(self):
    def filters(self):
    def filters(self, value):
    def add_filter(self, filter_name):
    def infos(self):
    def add_info(self, *args, **kwargs):
    def kind(self):
    def is_indel(self):
    def calldata_class(self):
    def copy_mapping_calls(self, call_mapper):
    def copy(self):
    def remove_gt_from_het_calls(self):
        def call_mapper(call):
    def remove_gt_from_low_qual_calls(self, min_qual):
        def call_mapper(call):
    def filter_calls_by_sample(self, samples, reverse=False, keep_info=False):
    def allele_freqs(self):
    def is_polymorphic(self, freq_threshold=1):
    def __init__(self, call, snv):
    def _get_allele_depths(self):
    def _get_depths_generic(self):
    def _get_depths_gatk(self):
    def _get_depths_varscan(self):
    def _get_depths_freebayes(self):
    def depth(self):
    def gt_qual(self):
    def ref_depth(self):
    def alt_sum_depths(self):
    def allele_depths(self):
    def has_alternative_counts(self):
    def maf_depth(self):
    def sample(self):
    def called(self):
    def gt_type(self):
    def is_het(self):
    def _to_int(allele):
    def int_alleles(self):
    def __str__(self):
    def __unicode__(self):
    def __repr__(self):
    def copy(self, return_pyvcf_call=False):
    def copy_setting_gt(self, gt, return_pyvcf_call=False):
    def copy_calldata(self, dict_with_data, return_pyvcf_call=False): MSG = 'A python package to run this executable is required,'
BIO = 'biopython'
BIO_BGZF = 'biopython with Bgzf support'
    NCBIXML = create_fake_class(MSG + BIO)
    NCBIWWW = create_fake_class(MSG + BIO)
def create_fake_class(msg):
        def __init__(self, *args, **kwargs):
def create_fake_funct(msg):
    def FakeRequiredfunct(*args, **kwargs):
    class FakePythonRequiredClass(object): class ReadRefCoord(object):
    def __init__(self, alig_read, sam, hard_clip_as_soft=False):
        self._alig_read = alig_read
        self._sam = sam
        self.hard_clip_as_soft = hard_clip_as_soft
        self._blocks = None
        self._read_len = None
        self._last_position = None

    @property
    def read_len(self):
        if self._read_len is None:
            self.blocks
        return self._read_len

    @property
    def blocks(self):
        alig_read = self._alig_read
        if self._blocks is not None:
            return self._block
        ref_start = alig_read.reference_start

        cigar_tuples = alig_read.cigartuples
        is_fwd = not(alig_read.is_reverse)
        if not is_fwd:
            cigar_tuples = reversed(cigar_tuples)

        read_pos = 0  # if is reversed this is not true, but we fixed later
        read_len = 0
        ref_pos = ref_start
        blocks = []
        for operation, length in alig_read.cigartuples:
            blk_start = ref_pos, read_pos
            if operation == 5:
                if self.hard_clip_as_soft:
                    operation = 4   # if like a soft clip
                else:
                    continue
            if operation in (0, 7, 8):  # match
                read_len += length
                ref_pos += length
                if is_fwd:
                    read_pos += length
                else:
                    read_pos -= length
            elif operation == 1:  # insertion in read
                read_len += length
                if is_fwd:
                    read_pos += length
                else:
                    read_pos -= length
            elif operation == 2 or operation == 3:    # deletion in read
                ref_pos += length
            elif operation == 4:    # soft clip
                # This is not a block to be returned because it is not aligned
                read_len += length
                if is_fwd:
                    read_pos += length
                else:
                    read_pos -= length
            elif operation == 6:
                continue
            blk_end = ref_pos - 1, read_pos - 1
            if operation != 4:
                ref_start = blk_start[0]
                ref_end = blk_end[0]
                read_start = blk_start[1]
                read_end = blk_end[1]
                if operation in (2, 3):
                    read_start, read_end = None, None
                elif operation == 1:
                    ref_start, ref_end = None, None
                block = Block(ref_start, ref_end, read_start, read_end)
                blocks.append(block)

        if not is_fwd:
            rev_blocks = []
            for blk in blocks:
                if blk.read_stop is not None and blk.read_start is not None:
                    blk_read_start = blk.read_start + read_len - 1
                    blk_read_end = blk.read_stop + read_len + 1
                else:
                    blk_read_start = None
                    blk_read_end = None
                block = Block(blk.ref_start, blk.ref_stop, blk_read_start,
                              blk_read_end)
                rev_blocks.append(block)
            blocks = rev_blocks
        self._block = blocks
        self._read_len = read_len
        return blocks

    def get_read_pos(self, ref_pos):
        if (self._last_position is not None and
            self._last_position[0] == ref_pos):
            return self._last_position[1]
        alig_read = self._alig_read
        sam = self._sam
        chrom = sam.getrname(alig_read.reference_id)
        if chrom != ref_pos[0]:
            msg = 'The aligned read is not aligned to the given chrom: '
            msg += chrom
            raise ValueError(msg)
        ref_pos = ref_pos[1]

        blocks = self.blocks

        block_with_ref_pos = None
        for block in blocks:
            if block.ref_start <= ref_pos <= block.ref_stop:
                block_with_ref_pos = block
                break
        if block_with_ref_pos is None:
            read_pos = None
        else:
            read_start = block_with_ref_pos.read_start
            read_stop = block_with_ref_pos.read_stop

            if read_start is None or read_stop is None:
                read_pos = None
            else:
                ref_start = block_with_ref_pos.ref_start
                ref_stop = block_with_ref_pos.ref_stop

                if alig_read.is_reverse:
                    read_pos = read_start - (ref_pos - ref_start)
                else:
                    read_pos = read_stop - (ref_stop - ref_pos)

        self._last_position = ref_pos, read_pos
        return read_pos

    def get_read_pos_counting_from_end(self, ref_pos):
        if self._last_position is None:
            read_pos = self.get_read_pos(ref_pos)
        else:
            read_pos = self._last_position[1]
        read_len = self.read_len
        if read_pos is not None:
            return read_pos - read_len
from operator import itemgetter
from crumbs.seq.seq import get_name, get_length
from crumbs.seq.seqio import read_seqs
from crumbs.statistics import IntCounter, IntBoxplot
from crumbs.plot import get_fig_and_canvas, draw_int_boxplot
from crumbs.vcf.snv import (VARSCAN, GATK, FREEBAYES, HOM_REF, HET, HOM_ALT,
                            HOM, DEF_MIN_CALLS_FOR_POP_STATS, VCFReader,
                            pyvcfReader)
from crumbs.utils.optional_modules import Reader
from crumbs.bam.coord_transforms import ReadRefCoord
            data[rc][acs] = {}
        if gt not in data[rc][acs]:
            data[rc][acs][gt] = {'num_gt': 0, 'sum_gq': 0}
        data[rc][acs][gt]['num_gt'] += 1
        data[rc][acs][gt]['sum_gq'] += gq

    def _get_data_for_gt_type(self, ref_count, alt_count, gt_type):
        genotypes = self._data.get(ref_count, {}).get(alt_count, None)
        if genotypes is None:
            return None

        gts = self._genotypes[gt_type]
        gt_count = 0
        gq_quals = 0
        for gt in gts:
            if gt in genotypes:
                gt_count += genotypes[gt]['num_gt']
                gq_quals += genotypes[gt]['sum_gq']

        if not gt_count:
            return None

        return gt_count, gq_quals

    def get_gt_count(self, ref_count, alt_count, gt_type):
        count_and_qual = self._get_data_for_gt_type(ref_count, alt_count,
                                                    gt_type)
        if count_and_qual is None:
            return None
        else:
            return count_and_qual[0]

    def get_avg_gt_qual(self, ref_count, alt_count, gt_type):
        counts_and_quals = self._get_data_for_gt_type(ref_count, alt_count,
                                                      gt_type)
        if counts_and_quals is None:
            return None

        gt_count, gq_quals = counts_and_quals

        return gq_quals / gt_count

    def get_avg_gt_quals(self, gt_type, rc_max=None, ac_max=None):
        for rc, allele_counts in self._data.items():
            if rc_max and rc > rc_max:
                continue
            for ac in allele_counts:
                if ac_max and ac > ac_max:
                    continue
                avg_qual = self.get_avg_gt_qual(rc, ac, gt_type)
                if avg_qual is not None:
                    yield rc, ac, avg_qual

    def get_gt_counts(self, gt_type, rc_max=None, ac_max=None):
        for rc, allele_counts in self._data.items():
            if rc_max and rc > rc_max:
                continue
            for ac in allele_counts:
                if ac_max and ac > ac_max:
                    continue
                gt_count = self.get_gt_count(rc, ac, gt_type)
                if gt_count is not None:
                    yield rc, ac, gt_count

    def get_gt_depths_for_coverage(self, coverage):
        for ref_count in range(coverage + 1):
            alt_count = coverage - ref_count

            genotypes = self._data.get(ref_count, {}).get(alt_count, None)
            if genotypes is None:
                continue
            yield ref_count, genotypes


class VcfStats(object):
    def __init__(self, vcf_fpath, gq_threshold=None, dp_threshold=100,
                 min_calls_for_pop_stats=DEF_MIN_CALLS_FOR_POP_STATS,
                 remarkable_coverages=None, window_size=WINDOWS_SIZE):
        if remarkable_coverages is None:
            remarkable_depths = REMARKABLE_DEPTHS
        self.remarkable_depths = remarkable_depths

        self._reader = VCFReader(open(vcf_fpath),
                               min_calls_for_pop_stats=min_calls_for_pop_stats)

        self._random_reader = pyvcfReader(filename=vcf_fpath)

        self.window_size = window_size
        self._gq_threshold = 0 if gq_threshold is None else gq_threshold

        self.dp_threshold = dp_threshold
        self._gt_qual_depth_counter = {HOM: IntBoxplot(), HET: IntBoxplot()}
        self._ac2d = _AlleleCounts2D()

        self.sample_dp_coincidence = {1: IntCounter()}
        for cov in remarkable_depths:
            self.sample_dp_coincidence[cov] = IntCounter()

        self.called_snvs = 0
        self.called_gts = IntCounter()

        # sample_counter
        self._sample_counters = {}

        for counter_name in SAMPLE_COUNTERS:
            if counter_name not in self._sample_counters:
                self._sample_counters[counter_name] = {}
            for sample in self._reader.samples:
                if counter_name in (GT_DEPTHS, GT_QUALS):
                    counters = {HOM: IntCounter(), HET: IntCounter()}
                else:
                    counters = IntCounter()
                self._sample_counters[counter_name][sample] = counters

        self._snv_counters = {MAFS: IntCounter(),
                              MACS: IntCounter(),
                              MAFS_DP: IntCounter(),
                              SNV_QUALS: IntCounter(),
                              HET_IN_SNP: IntCounter(),
                              SNV_DENSITY: IntCounter(),

                              INBREED_F_IN_SNP: IntCounter(),