TEST_DATA_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '..',
                                             '..', 'test', 'test_data')) VCF_BIN_DIR = os.path.join(BIN_DIR, 'vcf') class CloseToSnv(BaseAnnotator):
    '''Filter snps with other close snvs.

    Allowed snv_types:  [snp, indel, unknown] '''

    def __init__(self, distance=60, max_maf_depth=None, snv_type=None):
        self.distance = distance
        self.random_reader = None
        self.max_maf_depth = max_maf_depth
        self.snv_type = snv_type
        self.conf = {'distance': distance, 'max_maf_depth': max_maf_depth,
                     'snv_type': snv_type}

    def __call__(self, snv):
        if self.random_reader is None:
            self.random_reader = self._create_reader_from_snv(snv)

        self._clean_filter(snv)
        chrom = snv.chrom
        pos = snv.pos
        start = pos - self.distance if pos - self.distance > 0 else 0
        end = pos + self.distance
        snv_type = self.snv_type
        max_maf_depth = self.max_maf_depth
        passed_snvs = 0
        for snv_in_window in self.random_reader.fetch_snvs(chrom, start, end):
            if snv_in_window.pos == pos:
                continue
            if max_maf_depth is None and snv_type is None:
                passed_snvs += 1
            elif max_maf_depth is None and snv_type is not None:
                if snv_type == snv_in_window.var_type:
                    passed_snvs += 1
            elif max_maf_depth is not None and snv_type is None:
                calculated_maf = snv_in_window.maf_depth
                if calculated_maf and calculated_maf <= max_maf_depth:
                    passed_snvs += 1
            else:
                calculated_maf = snv_in_window.maf_depth
                if (calculated_maf and calculated_maf <= max_maf_depth and
                    snv_type == snv_in_window.kind):
                    passed_snvs += 1
        if passed_snvs:
            snv.add_filter(self.name)

        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        snv_type = '' if self.snv_type is None else self.snv_type[0]
        maf = self.max_maf_depth
        maf_str = '' if maf is None else '_{:.2f}'.format(maf)
        return 'cs{}{}{}'.format(snv_type, self.distance, maf_str)

    @property
    def description(self):
        snv_type = 'snv' if self.snv_type is None else self.snv_type
        maf = self.max_maf_depth
        maf_str = ', with maf:{:.2f}'.format(maf) if maf else ''
        desc = 'The snv is closer than {} nucleotides to another {}{} :: {}'
        return desc.format(self.distance, snv_type, maf_str, self.encode_conf)

    @property
    def is_filter(self):
        return True

    @property
    def item_wise(self):
        return True class HighVariableRegion(BaseAnnotator):
    'Filter depending on the variability of the region'

    def __init__(self, max_variability, window_in_bp, ref_fpath):
        self.max_variability = max_variability

        self._lengths = _get_lengths(open(ref_fpath))

        if not window_in_bp % 2:
            raise ValueError('Window in bp must be a odd number')
        self.window_in_bp = window_in_bp

        max_num_snps = (int(window_in_bp * max_variability) + 1) * 2
        if not max_num_snps % 2:
            max_num_snps += 1
        self._max_num_snps = max_num_snps

        self.conf = {'max_variability': max_variability,
                     'window_in_bp': window_in_bp}

    def __call__(self, snvs):
        # TODO: Randon acess iterator based on physical distance
        # RandomAcessRegionIterator(items, location_getter)
        # once we do this we can remove max_num_snps
        snvs = RandomAccessIterator(snvs, self._max_num_snps)
        half_win = (snvs._rnd_access_win - 1) // 2
        half_win_in_bp = (self.window_in_bp - 1) // 2
        for idx, snv in enumerate(snvs):
            self._clean_filter(snv)
            chrom = snv.chrom
            pos = snv.pos
            start = idx - half_win
            if start < 0:
                start = 0
            end = idx + half_win
            snvs_in_win = snvs[start:end]

            def snv_is_close(snv2):
                if snv2.chrom != chrom:
                    return False
                if abs(snv2.pos - pos) < half_win_in_bp:
                    return True
                else:
                    return False

            close_snvs = filter(snv_is_close, snvs_in_win)

            num_snvs = len(close_snvs)

            win_len = self.window_in_bp
            # The studied window could be smaller than expected if it is
            # located at the beginning of the chromosome
            dist_from_0 = pos
            if dist_from_0 < half_win_in_bp:
                win_len -= (half_win_in_bp - dist_from_0)

            # The studied window could be smaller than expected if it is
            # located at the end of the chromosome
            ref_len = self._lengths[chrom]
            end = ref_len - 1
            len_not_studied_at_end = pos + half_win_in_bp - end
            if len_not_studied_at_end > 0:
                win_len -= len_not_studied_at_end

            freq = num_snvs / win_len
            if freq >= self.max_variability:
                snv.add_filter(self.name)
            yield snv

    @property
    def name(self):
        return 'hv{}'.format(self.max_variability)

    @property
    def description(self):
        desc = 'The region has more than {} snvs per {} bases:: {}'
        return desc.format(int(self.max_variability * 100),
                           self.window_in_bp, self.encode_conf)

    @property
    def is_filter(self):
        return True

    @property
    def item_wise(self):
        return False class CloseToLimit(BaseAnnotator):
    'Filter snps close to chromosome limits'

    def __init__(self, distance, ref_fpath):
        self.distance = distance
        self._lengths = _get_lengths(open(ref_fpath))
        self.conf = {'distance': distance}

    def __call__(self, snv):
        self._clean_filter(snv)
        chrom = snv.chrom
        pos = snv.pos
        seq_len = self._lengths[chrom]
        if pos < self.distance or seq_len - pos < self.distance:
            snv.add_filter(self.name)
        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'cl{}'.format(self.distance)

    @property
    def description(self):
        desc = 'The snv is closer than {} nucleotides to the reference edge'
        desc += ':: {}'
        return desc.format(self.distance, self.encode_conf)

    @property
    def is_filter(self):
        return True

    @property
    def item_wise(self):
        return True class MafDepthLimit(BaseAnnotator):
    'Filter by maf'

    def __init__(self, max_maf, samples=None):
        self.max_maf = max_maf
        self.samples = samples
        self.conf = {'max_maf': max_maf, 'samples': samples}

    def __call__(self, snv):
        self._clean_filter(snv)
        maf = snv.maf_depth
        if max and maf > self.max_maf:
            snv.add_filter(self.name)
        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'maf{}'.format(self.max_maf)

    @property
    def description(self):
        samples = 'all' if self.samples is None else ','.join(self.samples)
        desc = 'The most frequent allele in samples: {}.'
        desc += 'frequency greater than {} :: {}'
        return desc.format(samples, self.max_maf, self.encode_conf)

    @property
    def is_filter(self):
        return True

    @property
    def item_wise(self):
        return True class CapEnzyme(BaseAnnotator):
    'Filters by detectablity by restriction enzymes'

    def __init__(self, all_enzymes, ref_fpath):
        self.all_enzymes = all_enzymes
        self.ref_index = seq_index(ref_fpath, 'fasta')
        self.conf = {'all_enzymes': all_enzymes}
        self._last_chrom = None

    def __call__(self, snv):
        self._clean_filter(snv)
        alleles = snv.alleles
        enzymes = set()
        # we have to make all the posible conbinations
        chrom = snv.chrom
        last_chrom = self._last_chrom
        if last_chrom is not None and chrom == self._last_chrom[0]:
            ref = last_chrom[1]
        else:
            ref = self.ref_index[snv.chrom]
            self._last_chrom = chrom, ref
        used_combinations = []
        for i_index in range(len(alleles)):
            for j_index in range(len(alleles)):
                allelei = alleles[i_index]
                allelej = alleles[j_index]
                if (i_index == j_index or
                    (allelei, allelej) in used_combinations or
                    (allelej, allelei) in used_combinations):
                    continue
                used_combinations.append((allelei, allelej))
                i_j_enzymes = _cap_enzymes_between_alleles(allelei, allelej,
                                                           ref, snv.pos,
                                                           snv.end,
                                                  all_enzymes=self.all_enzymes)
                enzymes = enzymes.union(i_j_enzymes)
        if not enzymes:
            snv.add_filter(self.name)
        else:
            enzymes = [str(e) for e in list(enzymes)]
            snv.add_info(info=self.info['id'], value=','.join(enzymes))

        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'ce{}'.format('t' if self.all_enzymes else 'f')

    @property
    def description(self):
        enzs = 'all' if self.all_enzymes else 'cheap_ones'
        desc = "SNV is not a CAP detectable by the enzymes: {}:: {}"
        return desc.format(enzs, self.encode_conf)

    @property
    def info(self):
        return {'id': 'CAP', 'num': 1, 'type': 'String',
                'desc': 'Enzymes that can be use to differentiate the alleles'}

    @property
    def is_filter(self):
        return True

    @property
    def item_wise(self):
        return True class AminoChangeAnnotator(BaseAnnotator):
    'The aminoacid changes with the alternative allele'

    def __init__(self, ref_fpath, orf_seq_fpath):
        self.orf_suffix = '_orf_seq'
        self.ref_index = seq_index(ref_fpath, 'fasta')
        self.orf_seq_index = seq_index(orf_seq_fpath, 'fasta')

        self.conf = {}

    def __call__(self, snv):
        self._clean_filter(snv)
        seq_name = snv.chrom
        seq_ref = self.ref_index[seq_name]
        aminos = None
        try:
            seq_estscan = self.orf_seq_index[seq_name + self.orf_suffix]
        except KeyError:
            return
        try:
            aminos = get_amino_change(seq_ref, seq_estscan, snv)
        except IsIndelError:
            snv.add_filter(self.name)
            snv.add_info(info=self.info_id, value='INDEL')
        except BetweenSegments:
            pass
        except OutsideAlignment:
            pass
        if aminos is None:
            return

        if set(aminos['ref_amino']) != set(aminos['alt_amino']):
            snv.add_filter(self.name)
            info_val = '{}->{}'.format(aminos['ref_amino'],
                                       ','.join(aminos['alt_amino']))
            snv.add_info(info=self.info_id, value=info_val)

        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'caa'

    @property
    def description(self):
        return "Alt alleles change the amino acid"

    @property
    def info(self):
        return {'id': 'AAC', 'num': 1, 'type': 'String',
                'desc': 'Amino acid change'}

    @property
    def is_filter(self):
        return False

    @property
    def item_wise(self):
        return True class IsVariableAnnotator(BaseAnnotator):
    _ids = count(0)

    def __init__(self, samples=None, consider_reference=False, filter_id=None,
                 min_samples_for_non_var=1):
        if filter_id is None:
            filter_id = self._ids.next()
        self.min_samples_for_non_var = min_samples_for_non_var
        self.filter_id = filter_id
        self.samples = samples
        self.consider_reference = consider_reference
        self.conf = {'samples': samples,
                     'consider_reference': consider_reference}

    def __call__(self, snv):
        samples = self.samples
        if samples is None:
            calls = snv.calls
        else:
            calls = [snv.get_call(sample) for sample in samples]

        alleles = set()
        if self.consider_reference:
            alleles.add(0)

        n_samples_called = 0
        for call in calls:
            int_alleles = call.int_alleles
            if int_alleles:
                n_samples_called += 1
            alleles.update(int_alleles)

        if not n_samples_called:
            is_variable = None
        elif len(alleles) == 1:
            if len(samples) >= self.min_samples_for_non_var:
                is_variable = False
            else:
                is_variable = None
        else:
            is_variable = True

        snv.add_info(info=self.info_id, value=str(is_variable))

        if self.return_modified_snv:
            return snv

    @property
    def info(self):
        if self.samples is None:
            samples = 'all'
        else:
            samples = ','.join(self.samples)
        desc = 'True if the genotypes are variable. False if they are not. '
        desc += 'None if there is not enough data in the samples: {samples}'
        desc += ' :: {conf}'
        desc = desc.format(samples=samples, conf=self.encode_conf)

        return {'id': 'IV{}'.format(self.filter_id), 'num': 1,
                'type': 'String', 'desc': desc}

    @property
    def is_filter(self):
        return False

    @property
    def item_wise(self):
        return True class IsVariableDepthAnnotator(BaseAnnotator):
    'Variable in readgroup using depths and not genotypes'
    _ids = count(0)

    def __init__(self, samples, filter_id=None, max_maf_depth=None,
                 min_reads=MIN_READS, in_union=True,
                 min_reads_per_allele=MIN_READS_PER_ALLELE,
                 in_all_groups=True, reference_free=True):
        if filter_id is None:
            filter_id = self._ids.next()
        self.filter_id = filter_id
        self.max_maf = max_maf_depth
        self.samples = samples
        self.min_reads = min_reads
        self.min_reads_per_allele = min_reads_per_allele
        self.in_union = in_union
        self.in_all_groups = in_all_groups
        self.reference_free = reference_free
        self.conf = {'max_maf': max_maf_depth, 'samples': samples,
                     'min_reads': min_reads, 'in_union': in_union,
                     'min_reads_per_allele': min_reads_per_allele,
                     'in_all_groups': in_all_groups,
                     'reference_free': reference_free}

    def __call__(self, snv):
        is_variable = variable_in_samples_depth(snv, self.samples,
                                                self.in_union,
                                                self.max_maf,
                                                self.in_all_groups,
                                                self.reference_free,
                                                self.min_reads,
                                                self.min_reads_per_allele)

        snv.add_info(info=self.info_id, value=str(is_variable))

        if self.return_modified_snv:
            return snv

    @property
    def info(self):
        samples = ','.join(self.samples)
        desc = 'True if the samples are variable. False if they are not. '
        desc += 'None if there is not enougth data in the samples: {samples}'
        desc += ' :: {conf}'
        desc = desc.format(samples=samples, conf=self.encode_conf)

        return {'id': 'IVD{}'.format(self.filter_id), 'num': 1,
                'type': 'String', 'desc': desc}

    @property
    def is_filter(self):
        return False

    @property
    def item_wise(self):
        return True class AminoSeverityChangeAnnotator(BaseAnnotator):
    'Check if the change is severe  or not'

    def __init__(self, ref_fpath, orf_seq_fpath):
        self.orf_suffix = '_orf_seq'
        self.ref_index = seq_index(ref_fpath, 'fasta')
        self.orf_seq_index = seq_index(orf_seq_fpath, 'fasta')
        blossum_path = join(DATA_DIR, 'blossum90.csv')
        self.blosum = _parse_blossum_matrix(blossum_path)

        self.conf = {}

    def _is_severe(self, ref_aa, alt_aa):
        return True if self.blosum[(ref_aa, alt_aa)] < 0 else False

    def __call__(self, snv):
        self._clean_filter(snv)
        seq_name = snv.chrom
        seq_ref = self.ref_index[seq_name]
        aminos = None
        try:
            seq_estscan = self.orf_seq_index[seq_name + self.orf_suffix]
        except KeyError:
            return
        try:
            aminos = get_amino_change(seq_ref, seq_estscan, snv)
        except IsIndelError:
            snv.add_filter(self.name)
        except BetweenSegments:
            pass
        except OutsideAlignment:
            pass

        if aminos is None:
            return

        ref_aa = aminos['ref_amino']
        alt_aas = aminos['alt_amino']
        if set(ref_aa) == set(alt_aas):
            return snv
        if any([self._is_severe(ref_aa, alt_aa) for alt_aa in alt_aas]):
            snv.add_filter(self.name)

        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'cas'

    @property
    def description(self):
        return "Alt alleles change the amino acid and the change is severe"

    @property
    def is_filter(self):
        return False

    @property
    def item_wise(self):
        return True class HeterozigoteInSamples(BaseAnnotator):

    def __init__(self, samples=None, filter_id=None, min_percent_het_gt=0,
                 gq_threshold=0, min_num_called=0):
        self._samples = samples
        self._min_percent_het_gt = min_percent_het_gt
        self._gq_threshold = gq_threshold
        self._min_num_called = min_num_called
        if filter_id is None:
            filter_id = self._ids.next()
        self.filter_id = filter_id
        self.conf = {'samples': samples,
                     'min_percent_het_get': min_percent_het_gt,
                     'gq_threslhold': gq_threshold,
                     'min_num_called': min_num_called}

    def __call__(self, snv):
        call_is_het = []
        for call in snv.calls:
            sample_name = call.sample
            if self._samples and sample_name not in self._samples:
                continue
            if call.gt_qual >= self._gq_threshold:
                call_is_het.append(call.is_het)
        num_calls = len(call_is_het)
        num_hets = len(filter(bool, call_is_het))
        if not num_calls or num_calls < self._min_num_called:
            result = None
        else:
            percent = int((num_hets / num_calls) * 100)
            if percent >= self._min_percent_het_gt:
                result = True
            else:
                result = False

        snv.add_info(info=self.info_id, value=str(result))

        if self.return_modified_snv:
            return snv

    @property
    def info(self):
        if self._samples is None:
            samples = 'all'
        else:
            samples = ','.join(self._samples)
        description = 'True if at least {min_perc}% of the called samples are '
        description += 'het. False if not. None if not enough data in the '
        description += 'samples {samples} :: {conf}'
        description = description.format(min_perc=self._min_percent_het_gt,
                                         samples=samples,
                                         conf=self.encode_conf)

        return {'id': 'HIS{}'.format(self.filter_id), 'num': 1,
                'type': 'String', 'desc': description}

    @property
    def is_filter(self):
        return False

    @property
    def item_wise(self):
        return True class LowComplexityRegionAnnotator(BaseAnnotator):
    def __init__(self, ref_fpath, dust_threshold=DEFATULT_DUST_THRESHOLD,
                 window=DEF_SNP_DUST_WINDOW):
        self.threshold = dust_threshold
        self.window = window
        self.ref_index = seq_index(ref_fpath, 'fasta')
        self._last_chrom = None
        self._scores = array('f')

    def __call__(self, snv):
        self._clean_filter(snv)
        # we have to make all the posible conbinations
        chrom = snv.chrom
        last_chrom = self._last_chrom
        if last_chrom is not None and chrom == self._last_chrom[0]:
            ref = last_chrom[1]
        else:
            ref = self.ref_index[snv.chrom]
            self._last_chrom = chrom, ref

        start, end = calculate_window(snv.pos, snv.end, self.window, len(ref))

        snv_win_seq = SeqWrapper(SEQRECORD, ref[start:end], None)
        score = calculate_dust_score(snv_win_seq)
        if score > self.threshold:
            snv.add_filter(self.name)

        self._scores.append(score)

        if self.return_modified_snv:
            return snv

    @property
    def name(self):
        return 'lcr'

    def draw_hist(self, fhand):
        fig = Figure()
        axes = fig.add_subplot(111)
        axes.hist(self._scores, fill=True, log=True, bins=20,
                  rwidth=1)
        axes.axvline(x=self.threshold)
        axes.set_xlabel('% complexity score')
        axes.set_ylabel('num. SNPs')
        _print_figure(axes, fig, fhand, plot_legend=False)

    @property
    def is_filter(self):
        return True

    @property
    def description(self):
        return "SNV is located in a low complexity region"

    @property
    def item_wise(self):
        return True class HihgCoverageRegionAnnotator(BaseAnnotator):
    def __init__(self, bam_fpaths, coverage_threshold, window, min_mapq=None):
        self._coverage_threshold = coverage_threshold
        self._window = window
        self._min_mapq = min_mapq
        self._bam_coverage = BamCoverages2(bam_fpaths, min_mapq=min_mapq,
                                           window=window)
        self._scores = {samp: IntCounter()
                        for samp in self._bam_coverage.samples}

    def __call__(self, snv):
        coverages = self._bam_coverage.calculate_coverage_in_pos(snv.chrom,
                                                                 snv.pos)
        for sample, coverage in coverages.items():
            self._scores[sample][coverage] += 1
        mean_coverage = sum(coverages.values()) / len(coverages)
        if mean_coverage > self._coverage_threshold:
            snv.add_filter(self.name)

        if self.return_modified_snv:
            return snv

    def draw_hist(self, fhand):
        scores = self._scores
        titles = ['aaaa'] * len(scores)
        plots_per_chart = 3
        num_cols = 2 if len(scores) > plots_per_chart else 1
        plot = HistogramPlotter(scores.values(), ylabel='Num. SNPs',
                                xlabel='Coverage', titles=titles, kind=LINE,
                                distrib_labels=scores.keys(),
                                num_cols=num_cols,
                                plots_per_chart=plots_per_chart)
        for axe in plot.axes:
            axe.axvline(self._coverage_threshold)
        plot.write_figure(fhand)

    @property
    def name(self):
        return 'hcr'

    @property
    def is_filter(self):
        return True

    @property
    def description(self):
        return "SNV is located in a high coverage region"

    @property
    def item_wise(self):
        return True class VCFReader(object):
    def __init__(self, fhand, compressed=None, filename=None,
                 min_calls_for_pop_stats=DEF_MIN_CALLS_FOR_POP_STATS):
        self.fhand = fhand
        self.pyvcf_reader = pyvcfReader(fsock=fhand, compressed=compressed,
                                        filename=filename)
        self.min_calls_for_pop_stats = min_calls_for_pop_stats
        self._snpcaller = None
        self._samples = self.pyvcf_reader.samples
        self._header_lines = self.pyvcf_reader._header_lines
        self._column_headers = self.pyvcf_reader._column_headers

    def parse_snvs(self):
        min_calls_for_pop_stats = self.min_calls_for_pop_stats
        last_snp = None
        try:
            counter =0
            for snp in self.pyvcf_reader:
                counter +=1
                snp = SNV(snp, reader=self,
                          min_calls_for_pop_stats=min_calls_for_pop_stats)
                last_snp = snp
                yield snp
        except Exception:
            from traceback import print_exception
            exc_type, exc_value, exc_traceback = sys.exc_info()

            print_exception(exc_type, exc_value, exc_traceback,
                            limit=20, file=sys.stderr)

            if last_snp is not None:
                chrom = str(last_snp.chrom)
                pos = last_snp.pos
                msg = 'Last parsed SNP was: {} {}\n'.format(chrom, pos + 1)
                sys.stderr.write(msg)
                raise

    def fetch_snvs(self, chrom, start, end=None):
        min_calls_for_pop_stats = self.min_calls_for_pop_stats
        try:
            snvs = self.pyvcf_reader.fetch(chrom, start + 1, end=end)
        except KeyError:
            snvs = []
        if snvs is None:
            snvs = []

        for snp in snvs:
            snp = SNV(snp, reader=self,
                      min_calls_for_pop_stats=min_calls_for_pop_stats)
            yield snp

    def sliding_windows(self, size, step=None, ref_fhand=None,
                        min_num_snps=DEF_MIN_NUM_SNPS_IN_WIN):
        random_snp_reader = VCFReader(open(self.fhand.name))
        sliding_window = _SNPSlidingWindow(snp_reader=random_snp_reader,
                                           win_size=size, win_step=step,
                                           min_num_snps=min_num_snps,
                                           ref_fhand=ref_fhand)
        for window in sliding_window.windows():
            yield window

    @property
    def snpcaller(self):
        if self._snpcaller is not None:
            return self._snpcaller

        metadata = self.pyvcf_reader.metadata
        if 'source' in metadata:
            if 'VarScan2' in metadata['source']:
                snpcaller = VARSCAN
            elif 'freebayes' in metadata['source'][0].lower():
                snpcaller = FREEBAYES
            else:
                snpcaller = GENERIC
        elif 'UnifiedGenotyper' in metadata:
            snpcaller = GATK
        else:
            snpcaller = GENERIC
        self._snpcaller = snpcaller
        return snpcaller

    @property
    def samples(self):
        return self._samples

    @property
    def filters(self):
        return self.pyvcf_reader.filters

    @property
    def infos(self):
        return self.pyvcf_reader.infos

    def _build_header(self, samples=None):
        if samples is None:
            samples = self.samples
        header = '\n'.join(self._header_lines)
        header += '\n#' + '\t'.join(self._column_headers)
        header += '\t' + '\t'.join(samples)
        return header

    @property
    def header(self):
        return self._build_header()

    def create_template_reader(self, samples=None):
        vcf = self._build_header(samples)
        fake_vcf = BytesIO(vcf)
        return VCFReader(fake_vcf) class TemporaryDir(object):
    '''This class creates temporary directories '''
    def __init__(self, suffix='', prefix='', directory=None):
        '''It initiates the class.'''
        self._name = tempfile.mkdtemp(suffix=suffix, prefix=prefix,
                                      dir=directory)

    @property
    def name(self):
        'Returns path to the dict'
        return self._name

    def close(self):
        '''It removes the temp dir'''
        if os.path.exists(self._name):
            shutil.rmtree(self._name)
import unittest
from os.path import join
from tempfile import NamedTemporaryFile
from subprocess import check_output
from StringIO import StringIO
from crumbs.utils.test_utils import TEST_DATA_DIR
from crumbs.utils.bin_utils import VCF_BIN_DIR
from crumbs.vcf.annotation import (CloseToSnv, HighVariableRegion,
                                   CloseToLimit, MafDepthLimit, CapEnzyme,
                                   AminoChangeAnnotator, IsVariableAnnotator,
                                   IsVariableDepthAnnotator,
                                   AminoSeverityChangeAnnotator,
                                   HeterozigoteInSamples,
                                   LowComplexityRegionAnnotator,
                                   HihgCoverageRegionAnnotator)
from crumbs.vcf.snv import VCFReader
from crumbs.utils.file_utils import TemporaryDir
        assert snv.infos[info_id] == 'True'

        annotator = IsVariableDepthAnnotator(filter_id=2,
                                             samples=['sample06_gbs'])
        records = list(VCFReader(open(FREEBAYES_MULTI_VCF_PATH),
                                 min_calls_for_pop_stats=1).parse_snvs())
        snv = records[0].copy()
        annotator(snv)
        info_id = annotator.info_id
        assert snv.infos[info_id] == 'False'

    def test_is_variable_annotator(self):

        vcf = '''#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 1 2 3 4 5 6 7
20\t14\t.\tA\tA\t29\tPASS\t.\tGT\t./.\t./.\t1/1\t0/0\t1/.\t1/1\t1/1'''
        vcf = VCF_HEADER + vcf
        fhand = StringIO(vcf)
        snvs = list(VCFReader(fhand).parse_snvs())
        snv = snvs[0]

        annotator = IsVariableAnnotator(samples=['1', '2'])

        snvc = snv.copy()
        annotator(snvc)
        assert snvc.infos['IV0'] == 'None'

        snvc = snv.copy()
        annotator = IsVariableAnnotator(filter_id=0)
        annotator(snvc)
        assert snvc.infos['IV0'] == 'True'

        snvc = snv.copy()
        annotator = IsVariableAnnotator(filter_id=0, samples=['6', '7'])
        annotator(snvc)
        assert snvc.infos['IV0'] == 'False'

        # min num samples for not variable
        snvc = snv.copy()
        annotator = IsVariableAnnotator(filter_id=0, samples=['6', '7'],
                                        min_samples_for_non_var=3)
        annotator(snvc)
        assert snvc.infos['IV0'] == 'None'

        # with reference as one allele
        snvc = snv.copy()
        annotator = IsVariableAnnotator(filter_id=0, samples=['6', '7'],
                                        consider_reference=True)
        annotator(snvc)
        assert snvc.infos['IV0'] == 'True'

        snvc = snv.copy()
        annotator = IsVariableAnnotator(filter_id=0, samples=['1', '2'],
                                        consider_reference=True)
        annotator(snvc)
        assert snvc.infos['IV0'] == 'None'


class TestComplexityAnnotator(unittest.TestCase):
    def test_complexity_filter(self):
        records = list(VCFReader(open(VCF_PATH)).parse_snvs())
        low_complexity = LowComplexityRegionAnnotator(ref_fpath=REF_PATH)

        snp = records[0].copy()
        low_complexity(snp)
        assert low_complexity.name not in snp.filters

        snp1 = records[1]
        low_complexity(snp1)
        assert low_complexity.name not in snp1.filters
        fhand = NamedTemporaryFile(suffix='.png')
        low_complexity.draw_hist(fhand)


class HighCoverageRegionAnnotator(unittest.TestCase):
    def test_high_coverage_region(self):
        bam_fpath = join(TEST_DATA_DIR, 'seqs.bam')
        annotator = HihgCoverageRegionAnnotator([bam_fpath],
                                                coverage_threshold=5,
                                                window=21)
        snv = FakeClass()
        snv.chrom = 'reference1'
        snv.pos = 200
        snv.filters = []
        annotator(snv)
        assert snv.filters == ['hcr']

        annotator = HihgCoverageRegionAnnotator([bam_fpath],
                                                coverage_threshold=10,
                                                window=21)
        snv = FakeClass()
        snv.chrom = 'reference1'
        snv.pos = 200
        snv.filters = []
        annotator(snv)
        assert not snv.filters

        annotator = HihgCoverageRegionAnnotator([bam_fpath],
                                                coverage_threshold=10,
                                                window=21)
        for pos in range(189, 230):
            snv = FakeClass()
            snv.chrom = 'reference1'
            snv.pos = pos
            snv.filters = []
            annotator(snv)
        for pos in range(330, 400):
            snv = FakeClass()
            snv.chrom = 'reference2'
            snv.pos = pos
            snv.filters = []
            annotator(snv)

        with NamedTemporaryFile(suffix='.png') as fhand:
            annotator.draw_hist(fhand)
            # raw_input(fhand.name)


class TestInfoMappers(unittest.TestCase):

    def test_hetegorigot_percent(self):

        het_in_samples = HeterozigoteInSamples(filter_id=1)