Normalizes variants and reruns SnpEFF for resulting VCF
518852
def normalize(in_file, data, passonly=False, normalize_indels=True, split_biallelic=True,
              rerun_effects=True, remove_oldeffects=False, nonrefonly=False, work_dir=None):
    
    if remove_oldeffects:
        out_file = "%s-noeff-nomultiallelic%s" % utils.splitext_plus(in_file)
    else:
        out_file = "%s-nomultiallelic%s" % utils.splitext_plus(in_file)
    if work_dir:
        out_file = os.path.join(work_dir, os.path.basename(out_file))
    if not utils.file_exists(out_file):
        if vcfutils.vcf_has_variants(in_file):
            ready_ma_file = _normalize(in_file, data, passonly=passonly,
                                       normalize_indels=normalize_indels,
                                       split_biallelic=split_biallelic,
                                       remove_oldeffects=remove_oldeffects,
                                       nonrefonly=nonrefonly,
                                       work_dir=work_dir)
            if rerun_effects:
                ann_ma_file, _ = effects.add_to_vcf(ready_ma_file, data)
                if ann_ma_file:
                    ready_ma_file = ann_ma_file
            utils.symlink_plus(ready_ma_file, out_file)
        else:
            utils.symlink_plus(in_file, out_file)
    return vcfutils.bgzip_and_index(out_file, data["config"])
97878
def standardize():
    

    def f(G, bim):
        G_out = standardize_snps(G)
        return G_out, bim

    return f
691958
def translateIndelGenotypeFromVCF(vcf_genotypes, ref):
    ''''''

    # indels
    def getPrefix(s1, s2):
        '''get common prefix of strings s1 and s2.'''
        n = min(len(s1), len(s2))
        for x in range(n):
            if s1[x] != s2[x]:
                return s1[:x]
        return s1[:n]

    def getSuffix(s1, s2):
        '''get common sufix of strings s1 and s2.'''
        n = min(len(s1), len(s2))
        if s1[-1] != s2[-1]:
            return ""
        for x in range(-2, -n - 1, -1):
            if s1[x] != s2[x]:
                return s1[x + 1:]
        return s1[-n:]

    def getGenotype(variant, ref):

        if variant == ref:
            return "*", 0

        if len(ref) > len(variant):
            # is a deletion
            if ref.startswith(variant):
                return "-%s" % ref[len(variant):], len(variant) - 1
            elif ref.endswith(variant):
                return "-%s" % ref[:-len(variant)], -1
            else:
                prefix = getPrefix(ref, variant)
                suffix = getSuffix(ref, variant)
                shared = len(prefix) + len(suffix) - len(variant)
                # print "-", prefix, suffix, ref, variant, shared, len(prefix), len(suffix), len(ref)
                if shared < 0:
                    raise ValueError()
                return "-%s" % ref[len(prefix):-(len(suffix) - shared)], len(prefix) - 1

        elif len(ref) < len(variant):
            # is an insertion
            if variant.startswith(ref):
                return "+%s" % variant[len(ref):], len(ref) - 1
            elif variant.endswith(ref):
                return "+%s" % variant[:len(ref)], 0
            else:
                prefix = getPrefix(ref, variant)
                suffix = getSuffix(ref, variant)
                shared = len(prefix) + len(suffix) - len(ref)
                if shared < 0:
                    raise ValueError()

                return "+%s" % variant[len(prefix):-(len(suffix) - shared)], len(prefix)
        else:
            assert 0, "snp?"

    # in pileup, the position refers to the base
    # after the coordinate, hence subtract 1
    # pos -= 1

    genotypes, offsets = [], []
    is_error = True

    for variant in vcf_genotypes:
        try:
            g, offset = getGenotype(variant, ref)
        except ValueError:
            break

        genotypes.append(g)
        if g != "*":
            offsets.append(offset)

    else:
        is_error = False

    if is_error:
        raise ValueError()

    assert len(set(offsets)) == 1, "multiple offsets for indel"
    offset = offsets[0]

    genotypes = "/".join(genotypes)
    return genotypes, offset
1096761
def _postprocess_variants(record_file, data, ref_file, out_file):
    
    if not utils.file_uptodate(out_file, record_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = ["dv_postprocess_variants.py", "--ref", ref_file,
                   "--infile", record_file, "--outfile", tx_out_file]
            do.run(cmd, "DeepVariant postprocess_variants %s" % dd.get_sample_name(data))
    return out_file
403685
def finalize_sv(orig_vcf, data, items):
    
    paired = vcfutils.get_paired(items)
    # For paired/somatic, attach combined calls to tumor sample
    if paired:
        sample_vcf = orig_vcf if paired.tumor_name == dd.get_sample_name(data) else None
    else:
        sample_vcf = "%s-%s.vcf.gz" % (utils.splitext_plus(orig_vcf)[0], dd.get_sample_name(data))
        sample_vcf = vcfutils.select_sample(orig_vcf, dd.get_sample_name(data), sample_vcf, data["config"])
    if sample_vcf:
        effects_vcf, _ = effects.add_to_vcf(sample_vcf, data, "snpeff")
    else:
        effects_vcf = None
    return effects_vcf or sample_vcf
1032462
def cli(variant_file, vep, split, outfile, verbose, silent, check_info,
        allele_symbol, logfile, loglevel):
    
    from vcf_parser import logger, init_log

    if not loglevel:
        if verbose:
            loglevel = 'INFO'

    init_log(logger, logfile, loglevel)
    nr_of_variants = 0
    start = datetime.now()
    
    # with open(variant_file, 'r', encoding="utf-8") as f:
    #     for line in f:
    #         if not line.startswith('#'):
    #             nr_of_variants += 1
    
    if variant_file == '-':
        logger.info("Start parsing variants from stdin")
        my_parser = VCFParser(
            fsock=sys.stdin, 
            split_variants=split,
            check_info=check_info, 
            allele_symbol=allele_symbol
        )
    else:
        logger.info("Start parsing variants from file {0}".format(variant_file))
        my_parser = VCFParser(
            infile = variant_file,
            split_variants=split, 
            check_info=check_info, 
            allele_symbol=allele_symbol
        )

    if outfile:
        f = open(outfile, 'w', encoding='utf-8')
        logger.info("Printing vcf to file {0}".format(outfile))
    
    if not silent:
        logger.info("Printing vcf to stdout")
    else:
        logger.info("Skip printing since silent is active")
    
    for line in my_parser.metadata.print_header():
        if outfile:
            f.write(line+'\n')
        else:
            if not silent:
                print(line)
    try:
        for variant in my_parser:
            variant_line = '\t'.join([variant[head] for head in my_parser.header])
            if outfile:
                f.write(variant_line + '\n')
            else:
                if not silent:
                    print(variant_line)
            nr_of_variants += 1
    except SyntaxError as e:
        print(e)

    logger.info('Number of variants: {0}'.format(nr_of_variants))
    logger.info('Time to parse file: {0}'.format(str(datetime.now() - start)))
117945
def annotate_snvs(adapter, vcf_obj):
    
    variants = {}
    
    for nr_variants, variant in enumerate(vcf_obj, 1):
        # Add the variant to current batch
        variants[get_variant_id(variant)] = variant
        # If batch len == 1000 we annotate the batch
        if (nr_variants % 1000) == 0:
            
            for var_obj in adapter.search_variants(list(variants.keys())):
                var_id = var_obj['_id']
                if var_id in variants:
                    annotate_variant(variants[var_id], var_obj)

            for variant_id in variants:
                yield variants[variant_id]
            
            variants = {}
    
    for var_obj in adapter.search_variants(list(variants.keys())):
        var_id = var_obj['_id']
        if var_id in variants:
            annotate_variant(variants[var_id], var_obj)

    for variant_id in variants:
        yield variants[variant_id]
632702
def _run_vardict_paired(align_bams, items, ref_file, assoc_files,
                          region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-paired-variants.vcf.gz" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            vrs = bedutils.population_variant_regions(items)
            target = shared.subset_variant_regions(vrs, region,
                                                   out_file, items=items, do_merge=True)
            paired = vcfutils.get_paired_bams(align_bams, items)
            if not _is_bed_file(target):
                vcfutils.write_empty_vcf(tx_out_file, config,
                                         samples=[x for x in [paired.tumor_name, paired.normal_name] if x])
            else:
                if not paired.normal_bam:
                    ann_file = _run_vardict_caller(align_bams, items, ref_file,
                                                   assoc_files, region, out_file)
                    return ann_file
                vardict = get_vardict_command(items[0])
                vcfstreamsort = config_utils.get_program("vcfstreamsort", config)
                compress_cmd = "| bgzip -c" if out_file.endswith("gz") else ""
                freq = float(utils.get_in(config, ("algorithm", "min_allele_fraction"), 10)) / 100.0
                # merge bed file regions as amplicon VarDict is only supported in single sample mode
                opts, var2vcf_opts = _vardict_options_from_config(items, config, out_file, target)
                fix_ambig_ref = vcfutils.fix_ambiguous_cl()
                fix_ambig_alt = vcfutils.fix_ambiguous_cl(5)
                remove_dup = vcfutils.remove_dup_cl()
                if any("vardict_somatic_filter" in tz.get_in(("config", "algorithm", "tools_off"), data, [])
                       for data in items):
                    somatic_filter = ""
                    freq_filter = ""
                else:
                    var2vcf_opts += " -M "  # this makes VarDict soft filter non-differential variants
                    somatic_filter = ("| sed 's/\\\\.*Somatic\\\\/Somatic/' "
                                      "| sed 's/REJECT,Description=\".*\">/REJECT,Description=\"Not Somatic via VarDict\">/' "
                                      """| %s -c 'from bcbio.variation import freebayes; freebayes.call_somatic("%s", "%s")' """
                                      % (sys.executable, paired.tumor_name, paired.normal_name))
                    freq_filter = ("| bcftools filter -m '+' -s 'REJECT' -e 'STATUS !~ \".*Somatic\"' 2> /dev/null "
                                   "| %s -x 'bcbio.variation.vardict.add_db_germline_flag(x)' "
                                   "| %s "
                                   "| %s -x 'bcbio.variation.vardict.depth_freq_filter(x, %s, \"%s\")'" %
                                   (os.path.join(os.path.dirname(sys.executable), "py"),
                                    _lowfreq_linear_filter(0, True),
                                    os.path.join(os.path.dirname(sys.executable), "py"),
                                    0, bam.aligner_from_header(paired.tumor_bam)))
                jvm_opts = _get_jvm_opts(items[0], tx_out_file)
                py_cl = os.path.join(utils.get_bcbio_bin(), "py")
                setup = ("%s && unset JAVA_HOME &&" % utils.get_R_exports())
                contig_cl = vcfutils.add_contig_to_header_cl(ref_file, tx_out_file)
                cmd = ("{setup}{jvm_opts}{vardict} -G {ref_file} "
                       "-N {paired.tumor_name} -b \"{paired.tumor_bam}|{paired.normal_bam}\" {opts} "
                       "| awk 'NF>=48' | testsomatic.R "
                       "| var2vcf_paired.pl -P 0.9 -m 4.25 {var2vcf_opts} "
                       "-N \"{paired.tumor_name}|{paired.normal_name}\" "
                       "| {contig_cl} {freq_filter} "
                       "| bcftools filter -i 'QUAL >= 0' "
                       "{somatic_filter} | {fix_ambig_ref} | {fix_ambig_alt} | {remove_dup} | {vcfstreamsort} "
                       "{compress_cmd} > {tx_out_file}")
                do.run(cmd.format(**locals()), "Genotyping with VarDict: Inference", {})
    return out_file
830153
def vcf2pileup(vcf, sample):
    ''''''

    chromosome = vcf.contig
    pos = vcf.pos
    reference = vcf.ref
    allelles = [reference] + vcf.alt

    data = vcf[sample]

    # get genotype
    genotypes = data["GT"]
    if len(genotypes) > 1:
        raise ValueError("only single genotype per position, %s" % (str(vcf)))

    genotypes = genotypes[0]

    # not a variant
    if genotypes[0] == ".":
        return None

    genotypes = [allelles[int(x)] for x in genotypes if x != "/"]

    # snp_quality is "genotype quality"
    snp_quality = consensus_quality = data.get("GQ", [0])[0]
    mapping_quality = vcf.info.get("MQ", [0])[0]
    coverage = data.get("DP", 0)

    if len(reference) > 1 or max([len(x) for x in vcf.alt]) > 1:
        # indel
        genotype, offset = translateIndelGenotypeFromVCF(genotypes, reference)

        return PileupIndel(chromosome,
                           pos + offset,
                           "*",
                           genotype,
                           consensus_quality,
                           snp_quality,
                           mapping_quality,
                           coverage,
                           genotype,
                           "<" * len(genotype),
                           0,
                           0,
                           0)

    else:
        genotype = encodeGenotype("".join(genotypes))
        read_bases = ""
        base_qualities = ""

        return PileupSubstitution(chromosome, pos, reference,
                                  genotype, consensus_quality,
                                  snp_quality, mapping_quality,
                                  coverage, read_bases,
                                  base_qualities)
54111
def _dict_of_vars_to_vcf_file(variants, outfile):
    ''''''
    header_lines = [
        '##fileformat=VCFv4.2',
        '##source=cluster_vcf_records, version ' + cluster_vcf_records_version,
        '##fileDate=' + str(datetime.date.today()),
        '\t'.join(['#CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO'])
    ]

    with open(outfile, 'w') as f:
        print(*header_lines, sep='\n', file=f)

        for ref_name in sorted(variants):
            for pos in sorted(variants[ref_name]):
                for ref_string in sorted(variants[ref_name][pos]):
                    alts = sorted(list(variants[ref_name][pos][ref_string]))
                    print(ref_name, pos + 1, '.', ref_string, ','.join(alts), '.', 'PASS', 'SVTYPE=MERGED', sep='\t', file=f)
270084
def variants(ctx, variant_id, chromosome, end_chromosome, start, end, variant_type,
             sv_type):
    """"""
    if sv_type:
        variant_type = 'sv'

    adapter = ctx.obj['adapter']

    if (start or end):
        if not (chromosome and start and end):
            LOG.warning("Regions must be specified with chromosome, start and end")
            return

    if variant_id:
        variant = adapter.get_variant({'_id':variant_id})
        if variant:
            click.echo(variant)
        else:
            LOG.info("Variant {0} does not exist in database".format(variant_id))
        return

    if variant_type == 'snv':
        result = adapter.get_variants(
            chromosome=chromosome,
            start=start,
            end=end
        )
    else:
        LOG.info("Search for svs")
        result = adapter.get_sv_variants(
            chromosome=chromosome,
            end_chromosome=end_chromosome,
            sv_type=sv_type,
            pos=start,
            end=end
        )

    i = 0
    for variant in result:
        i += 1
        pp(variant)

    LOG.info("Number of variants found in database: %s", i)
1164831
def _upgrade_snpeff_data(galaxy_dir, args, remotes):
    
    snpeff_version = effects.snpeff_version(args)
    if not snpeff_version:
        return
    for dbkey, ref_file in genome.get_builds(galaxy_dir):
        resource_file = os.path.join(os.path.dirname(ref_file), "%s-resources.yaml" % dbkey)
        if os.path.exists(resource_file):
            with open(resource_file) as in_handle:
                resources = yaml.safe_load(in_handle)
            snpeff_db, snpeff_base_dir = effects.get_db({"genome_resources": resources,
                                                         "reference": {"fasta": {"base": ref_file}}})
            if snpeff_db:
                snpeff_db_dir = os.path.join(snpeff_base_dir, snpeff_db)
                if os.path.exists(snpeff_db_dir) and _is_old_database(snpeff_db_dir, args):
                    shutil.rmtree(snpeff_db_dir)
                if not os.path.exists(snpeff_db_dir):
                    print("Installing snpEff database %s in %s" % (snpeff_db, snpeff_base_dir))
                    dl_url = remotes["snpeff_dl_url"].format(
                        snpeff_ver=snpeff_version.replace(".", "_"),
                        genome=snpeff_db)
                    dl_file = os.path.basename(dl_url)
                    with utils.chdir(snpeff_base_dir):
                        subprocess.check_call(["wget", "--no-check-certificate", "-c", "-O", dl_file, dl_url])
                        subprocess.check_call(["unzip", dl_file])
                        os.remove(dl_file)
                    dl_dir = os.path.join(snpeff_base_dir, "data", snpeff_db)
                    shutil.move(dl_dir, snpeff_db_dir)
                    os.rmdir(os.path.join(snpeff_base_dir, "data"))
                if args.cwl:
                    create.directory_tarball(snpeff_db_dir)
1183613
def _replace_reference(self, var):
        """"""

        if var.type not in "cgmnr":
            raise HGVSUnsupportedOperationError("Can only update references for type c, g, m, n, r")

        if var.posedit.edit.type == "ins":
            # insertions have no reference sequence (zero-width), so return as-is
            return var
        if var.posedit.edit.type == "con":
            # conversions have no reference sequence (zero-width), so return as-is
            return var

        pos = var.posedit.pos
        if ((isinstance(pos.start, hgvs.location.BaseOffsetPosition) and pos.start.offset != 0)
                or (isinstance(pos.end, hgvs.location.BaseOffsetPosition) and pos.end.offset != 0)):
            _logger.info("Can't update reference sequence for intronic variant {}".format(var))
            return var

        # For c. variants, we need coords on underlying sequences
        if var.type == "c":
            tm = self._fetch_AlignmentMapper(
                tx_ac=var.ac, alt_ac=var.ac, alt_aln_method="transcript")
            pos = tm.c_to_n(var.posedit.pos)
        else:
            pos = var.posedit.pos
        seq = self.hdp.get_seq(var.ac, pos.start.base - 1, pos.end.base)

        edit = var.posedit.edit
        if edit.ref != seq:
            _logger.debug("Replaced reference sequence in {var} with {seq}".format(
                var=var, seq=seq))
            edit.ref = seq

        return var
190783
def _convert_to_sequence_variants(self, variant, acc):
        
        start = variant['start']
        insertion = ''.join(variant['ins'])
        deletion = ''.join(variant['del'])

        # defaults
        is_dup = False    # assume not dup
        fsext_len = None    # fs or ext length
        is_sub = False
        is_ext = False

        if start == 1:    # initial methionine is modified
            aa_start = aa_end = AAPosition(base=start, aa=deletion)
            ref = ''
            alt = ''
            self._is_init_met = True
            self._is_ambiguous = True    # side-effect

        if insertion and insertion.find("*") == 0:    # stop codon at variant position
            aa_start = aa_end = AAPosition(base=start, aa=deletion[0])
            ref = ''
            alt = '*'
            is_sub = True

        elif start == len(self._ref_seq):    # extension
            if self._alt_seq[-1] == '*':
                fsext_len = len(insertion) - len(deletion)    # don't include the former stop codon
            else:
                fsext_len = '?'
            subst_at_stop_codon = insertion[0]

            aa_start = aa_end = AAPosition(base=start, aa='*')
            ref = ''
            alt = subst_at_stop_codon
            is_ext = True

        elif self._is_frameshift:    # frameshift
            aa_start = aa_end = AAPosition(base=start, aa=deletion[0])
            ref = ''

            try:
                fsext_len = str(
                    insertion.index("*") + 1)    # start w/ 1st change; ends w/ * (inclusive)
            except ValueError:
                fsext_len = "?"

            alt = insertion[0]

        else:    # no frameshift - sub/delins/dup
            if insertion == deletion:    # silent
                aa_start = aa_end = AAPosition(base=start, aa=deletion)
                ref = alt = ''

            elif len(insertion) == len(deletion) == 1:    # substitution
                aa_start = aa_end = AAPosition(base=start, aa=deletion)
                ref = ''
                alt = insertion
                is_sub = True

            elif len(deletion) > 0:    # delins OR deletion OR stop codon at variant position
                ref = deletion
                end = start + len(deletion) - 1
                if len(insertion) > 0:    # delins
                    aa_start = AAPosition(base=start, aa=deletion[0])
                    if end > start:
                        aa_end = AAPosition(base=end, aa=deletion[-1])
                    else:
                        aa_end = aa_start
                    alt = insertion

                else:    # deletion OR stop codon at variant position
                    if len(deletion) + start == len(
                            self._ref_seq):    # stop codon at variant position
                        aa_start = AAPosition(base=start, aa=deletion[0])
                        aa_end = AAPosition(base=start, aa=deletion[0])
                        ref = ''
                        alt = '*'
                        is_sub = True
                    else:    # deletion
                        aa_start = AAPosition(base=start, aa=deletion[0])
                        if end > start:
                            aa_end = AAPosition(base=end, aa=deletion[-1])
                        else:
                            aa_end = aa_start
                        alt = None

            elif len(deletion) == 0:    # insertion OR duplication OR extension

                is_dup, dup_start = self._check_if_ins_is_dup(start, insertion)

                if is_dup:    # duplication
                    dup_end = dup_start + len(insertion) - 1
                    aa_start = AAPosition(base=dup_start, aa=insertion[0])
                    aa_end = AAPosition(base=dup_end, aa=insertion[-1])
                    ref = alt = None

                else:    # insertion
                    start -= 1
                    end = start + 1

                    aa_start = AAPosition(base=start, aa=self._ref_seq[start - 1])
                    aa_end = AAPosition(base=end, aa=self._ref_seq[end - 1])
                    ref = None
                    alt = insertion

            else:    # should never get here
                raise ValueError("unexpected variant: {}".format(variant))

        var_p = self._create_variant(
            aa_start,
            aa_end,
            ref,
            alt,
            fsext_len=fsext_len,
            is_dup=is_dup,
            acc=acc,
            is_ambiguous=self._is_ambiguous,
            is_sub=is_sub,
            is_ext=is_ext,
            is_init_met=self._is_init_met
        )

        return var_p
679855
def _seg_to_vcf(vals):
    
    call_to_cn = {"+": 3, "-": 1}
    call_to_type = {"+": "DUP", "-": "DEL"}
    if vals["CALL"] not in ["0"]:
        info = ["FOLD_CHANGE_LOG=%s" % vals["MEAN_LOG2_COPY_RATIO"],
                "PROBES=%s" % vals["NUM_POINTS_COPY_RATIO"],
                "SVTYPE=%s" % call_to_type[vals["CALL"]],
                "SVLEN=%s" % (int(vals["END"]) - int(vals["START"])),
                "END=%s" % vals["END"],
                "CN=%s" % call_to_cn[vals["CALL"]]]
        return [vals["CONTIG"], vals["START"], ".", "N", "<%s>" % call_to_type[vals["CALL"]], ".",
                ".", ";".join(info), "GT", "0/1"]
1245523
def to_vcf(in_file, caller, header_fn, vcf_fn, data, sep="\t"):
    
    out_file = "%s.vcf" % utils.splitext_plus(in_file)[0]
    if not utils.file_exists(out_file + ".gz") and not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            with open(in_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    out_handle.write(_vcf_header.format(caller=caller))
                    out_handle.write("\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL",
                                                "FILTER", "INFO", "FORMAT", dd.get_sample_name(data)]) + "\n")
                    header, in_handle = header_fn(in_handle)
                    for line in in_handle:
                        out = vcf_fn(dict(zip(header, line.strip().split(sep))))
                        if out:
                            out_handle.write("\t".join(out) + "\n")
    out_file = vcfutils.bgzip_and_index(out_file, data["config"])
    effects_vcf, _ = effects.add_to_vcf(out_file, data, "snpeff")
    return effects_vcf or out_file
1108010
def vcf_to_df_worker(arg):
    
    canvasvcf, exonbed, i = arg
    logging.debug("Working on job {}: {}".format(i, canvasvcf))
    samplekey = op.basename(canvasvcf).split(".")[0].rsplit('_', 1)[0]
    d = {'SampleKey': samplekey}

    exons = BedTool(exonbed)
    cn = parse_segments(canvasvcf)
    overlaps = exons.intersect(cn, wao=True)
    gcn_store = {}
    for ov in overlaps:
        # Example of ov.fields:
        # [u'chr1', u'11868', u'12227', u'ENSG00000223972.5',
        # u'ENST00000456328.2', u'transcribed_unprocessed_pseudogene',
        # u'DDX11L1', u'.', u'-1', u'-1', u'.', u'0']
        gene_name = "|".join((ov.fields[6], ov.fields[3], ov.fields[5]))
        if gene_name not in gcn_store:
            gcn_store[gene_name] = defaultdict(int)

        cn = ov.fields[-2]
        if cn == ".":
            continue
        cn = int(cn)
        if cn > 10:
            cn = 10
        amt = int(ov.fields[-1])
        gcn_store[gene_name][cn] += amt

    for k, v in sorted(gcn_store.items()):
        v_mean, v_median = counter_mean_and_median(v)
        d[k + ".avgcn"] = v_mean
        d[k + ".medcn"] = v_median
    cleanup()
    return d