Top level entry to regional de-novo based variant calling with cortex_var.
186091
def run_cortex(align_bams, items, ref_file, assoc_files, region=None,
               out_file=None):
    
    raise NotImplementedError("Cortex currently out of date and needs reworking.")
    if len(align_bams) == 1:
        align_bam = align_bams[0]
        config = items[0]["config"]
    else:
        raise NotImplementedError("Need to add multisample calling for cortex_var")
    if out_file is None:
        out_file = "%s-cortex.vcf" % os.path.splitext(align_bam)[0]
    if region is not None:
        work_dir = safe_makedir(os.path.join(os.path.dirname(out_file),
                                             region.replace(".", "_")))
    else:
        work_dir = os.path.dirname(out_file)
    if not file_exists(out_file):
        bam.index(align_bam, config)
        variant_regions = config["algorithm"].get("variant_regions", None)
        if not variant_regions:
            raise ValueError("Only support regional variant calling with cortex_var: set variant_regions")
        target_regions = subset_variant_regions(variant_regions, region, out_file)
        if os.path.isfile(target_regions):
            with open(target_regions) as in_handle:
                regional_vcfs = [_run_cortex_on_region(x.strip().split("\t")[:3], align_bam,
                                                       ref_file, work_dir, out_file, config)
                                 for x in in_handle]

            combine_file = "{0}-raw{1}".format(*os.path.splitext(out_file))
            _combine_variants(regional_vcfs, combine_file, ref_file, config)
            _select_final_variants(combine_file, out_file, config)
        else:
            vcfutils.write_empty_vcf(out_file)
    return out_file
66493
def use_vqsr(algs, call_file=None):
    
    from bcbio.variation import vcfutils
    vqsr_callers = set(["gatk", "gatk-haplotype"])
    vqsr_sample_thresh = 50
    vqsr_supported = collections.defaultdict(int)
    coverage_intervals = set([])
    for alg in algs:
        callers = alg.get("variantcaller")
        if isinstance(callers, six.string_types):
            callers = [callers]
        if not callers:  # no variant calling, no VQSR
            continue
        if "vqsr" in (alg.get("tools_off") or []):  # VQSR turned off
            continue
        for c in callers:
            if c in vqsr_callers:
                if "vqsr" in (alg.get("tools_on") or []):  # VQSR turned on:
                    vqsr_supported[c] += 1
                    coverage_intervals.add("genome")
                # Do not try VQSR for gVCF inputs
                elif call_file and vcfutils.is_gvcf_file(call_file):
                    pass
                else:
                    coverage_intervals.add(alg.get("coverage_interval", "exome").lower())
                    vqsr_supported[c] += 1
    if len(vqsr_supported) > 0:
        num_samples = max(vqsr_supported.values())
        if "genome" in coverage_intervals or num_samples >= vqsr_sample_thresh:
            return True
    return False
981159
@Override
	public void initVariantProviderStrategy(GeneratorContext context, Map<String, VariantSet> variantsSetMap) {

		List<Map<String, String>> variantMapStrategies = new ArrayList<>();
		Map<String, String> ctxVariantMap = context.getVariantMap();
		VariantSet skinVariantSet = variantsSetMap.get(JawrConstant.SKIN_VARIANT_TYPE);
		String skinVariant = ctxVariantMap.get(JawrConstant.SKIN_VARIANT_TYPE);

		VariantSet localeVariantSet = variantsSetMap.get(JawrConstant.LOCALE_VARIANT_TYPE);
		String localeVariant = ctxVariantMap.get(JawrConstant.LOCALE_VARIANT_TYPE);

		variantMapStrategies.add(getVariantMap(skinVariant, localeVariant));
		if (localeVariantSet != null) {
			variantMapStrategies.add(getVariantMap(skinVariant, localeVariantSet.getDefaultVariant()));
		}
		variantMapStrategies.add(getVariantMap(skinVariantSet.getDefaultVariant(), localeVariant));
		if (localeVariantSet != null) {
			variantMapStrategies
					.add(getVariantMap(skinVariantSet.getDefaultVariant(), localeVariantSet.getDefaultVariant()));
		}

		variantMapStrategyIterator = variantMapStrategies.iterator();
	}
1327429
def get_variants_in_region(self, chrom, start, end):
        """"""
        if self.chrom is not None and chrom == self.chrom:
            # We are going to search for 'NA' since the chromosome was set
            chrom = "NA"

        iterator = self._bgen.iter_variants_in_region(
            CHROM_STR_DECODE.get(chrom, chrom), start, end,
        )
        for info, dosage in iterator:
            yield Genotypes(
                Variant(
                    info.name, CHROM_STR_ENCODE.get(info.chrom, info.chrom),
                    info.pos, [info.a1, info.a2],
                ),
                dosage,
                reference=info.a1,
                coded=info.a2,
                multiallelic=True,
            )
868830
def _run_variantcall_batch_multicore(items, regions, final_file):
    
    batch_name = _get_batch_name(items)
    variantcaller = _get_batch_variantcaller(items)
    work_bams = [dd.get_work_bam(d) or dd.get_align_bam(d) for d in items]
    def split_fn(data):
        out = []
        for region in regions:
            region = _region_to_coords(region)
            chrom, start, end = region
            region_str = "_".join(str(x) for x in region)
            out_file = os.path.join(dd.get_work_dir(items[0]), variantcaller, chrom,
                                    "%s-%s.vcf.gz" % (batch_name, region_str))
            out.append((region, work_bams, out_file))
        return final_file, out
    parallel = {"type": "local", "num_jobs": dd.get_num_cores(items[0]), "cores_per_job": 1}
    run_parallel = dmulti.runner(parallel, items[0]["config"])
    to_run = copy.deepcopy(items[0])
    to_run["sam_ref"] = dd.get_ref_file(to_run)
    to_run["group_orig"] = items
    parallel_split_combine([[to_run]], split_fn, run_parallel,
                           "variantcall_sample", "concat_variant_files",
                           "vrn_file", ["region", "sam_ref", "config"])
    return final_file
456041
function createSequenceVars()
    {
        if ( !$this->hasSequence() )
            return;

        $fName      = $this->Name;
        $uniqid     = $this->UniqID;
        $this->NewNodes[] = eZTemplateNodeTool::createCodePieceNode( "// creating sequence variables for \{$fName} loop" );
        $this->NewNodes[] = eZTemplateNodeTool::createVariableNode( false,
                                                                    $this->Parameters['sequence_array'],
                                                                    $this->NodePlacement,
                                                                    array( 'treat-value-as-non-object' => true, 'text-result' => false ),
                                                                    "${fName}_sequence_array_$uniqid" );
        $this->NewNodes[] = eZTemplateNodeTool::createCodePieceNode( "\$${fName}_sequence_var_$uniqid = current( \$${fName}_sequence_array_$uniqid );\n" );
    }
722209
def general_conv(x,
                 num_filters=64,
                 filter_size=7,
                 stride=1,
                 stddev=0.02,
                 padding="VALID",
                 name="conv",
                 do_norm="instance",
                 do_relu=True,
                 relufactor=0):
  """"""
  with tf.variable_scope(name):
    x = layers().Conv2D(
        num_filters,
        filter_size,
        stride,
        padding,
        activation=None,
        kernel_initializer=tf.truncated_normal_initializer(stddev=stddev),
        bias_initializer=tf.constant_initializer(0.0))(x)
    if do_norm == "layer":
      x = layer_norm(x)
    elif do_norm == "instance":
      x = instance_norm(x)

    if do_relu:
      if relufactor == 0:
        x = tf.nn.relu(x, "relu")
      else:
        x = lrelu(x, leak=relufactor)

    return x
172966
def run(call_file, ref_file, vrn_files, data):
    
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if includes_missingalt(data):
        logger.info("Removing variants with missing alts from %s." % call_file)
        call_file = gatk_remove_missingalt(call_file, data)

    if "gatkcnn" in dd.get_tools_on(data):
        return _cnn_filter(call_file, vrn_files, data)
    elif config_utils.use_vqsr(algs, call_file):
        if vcfutils.is_gvcf_file(call_file):
            raise ValueError("Cannot force gVCF output with joint calling using tools_on: [gvcf] and use VQSR. "
                             "Try using cutoff-based soft filtering with tools_off: [vqsr]")
        snp_file, indel_file = vcfutils.split_snps_indels(call_file, ref_file, data["config"])
        snp_filter_file = _variant_filtration(snp_file, ref_file, vrn_files, data, "SNP",
                                              vfilter.gatk_snp_cutoff)
        indel_filter_file = _variant_filtration(indel_file, ref_file, vrn_files, data, "INDEL",
                                                vfilter.gatk_indel_cutoff)
        orig_files = [snp_filter_file, indel_filter_file]
        out_file = "%scombined.vcf.gz" % os.path.commonprefix(orig_files)
        combined_file = vcfutils.combine_variant_files(orig_files, out_file, ref_file, data["config"])
        return combined_file
    else:
        snp_filter = vfilter.gatk_snp_cutoff(call_file, data)
        indel_filter = vfilter.gatk_indel_cutoff(snp_filter, data)
        return indel_filter
258465
def _run_purecn(paired, work_dir):
    
    segfns = {"cnvkit": _segment_normalized_cnvkit, "gatk-cnv": _segment_normalized_gatk}
    out_base, out, all_files = _get_purecn_files(paired, work_dir)
    failed_file = out_base + "-failed.log"
    cnr_file = tz.get_in(["depth", "bins", "normalized"], paired.tumor_data)
    if not utils.file_uptodate(out["rds"], cnr_file) and not utils.file_exists(failed_file):
        cnr_file, seg_file = segfns[cnvkit.bin_approach(paired.tumor_data)](cnr_file, work_dir, paired)
        from bcbio import heterogeneity
        vcf_file = heterogeneity.get_variants(paired.tumor_data, include_germline=False)[0]["vrn_file"]
        vcf_file = germline.filter_to_pass_and_reject(vcf_file, paired, out_dir=work_dir)
        with file_transaction(paired.tumor_data, out_base) as tx_out_base:
            # Use UCSC style naming for human builds to support BSgenome
            genome = ("hg19" if dd.get_genome_build(paired.tumor_data) in ["GRCh37", "hg19"]
                      else dd.get_genome_build(paired.tumor_data))
            cmd = ["PureCN.R", "--seed", "42", "--out", tx_out_base, "--rds", "%s.rds" % tx_out_base,
                   "--sampleid", dd.get_sample_name(paired.tumor_data),
                   "--genome", genome,
                   "--vcf", vcf_file, "--tumor", cnr_file,
                   "--segfile", seg_file, "--funsegmentation", "Hclust", "--maxnonclonal", "0.3"]
            if dd.get_num_cores(paired.tumor_data) > 1:
                cmd += ["--cores", str(dd.get_num_cores(paired.tumor_data))]
            try:
                cmd = "export R_LIBS_USER=%s && %s && %s" % (utils.R_sitelib(), utils.get_R_exports(),
                                                             " ".join([str(x) for x in cmd]))
                do.run(cmd, "PureCN copy number calling")
            except subprocess.CalledProcessError as msg:
                if _allowed_errors(str(msg)):
                    logger.info("PureCN failed to find solution for %s: skipping" %
                                dd.get_sample_name(paired.tumor_data))
                    with open(failed_file, "w") as out_handle:
                        out_handle.write(str(msg))
                else:
                    logger.exception()
                    raise
            for f in all_files:
                if os.path.exists(os.path.join(os.path.dirname(tx_out_base), f)):
                    shutil.move(os.path.join(os.path.dirname(tx_out_base), f),
                                os.path.join(os.path.dirname(out_base), f))
    out = _get_purecn_files(paired, work_dir, require_exist=True)[1]
    return out if (out.get("rds") and os.path.exists(out["rds"])) else None
321553
def _variant_sv(checkpoints):
    
    if not checkpoints.get("sv"):
        return [], []
    sv = [s("detect_sv", "batch-single",
            [["sv_batch_rec"]],
            [cwlout("sv_rec", "record",
                    fields=[cwlout(["sv", "variantcaller"], ["string", "null"]),
                            cwlout(["sv", "vrn_file"], ["File", "null"], [".tbi"]),
                            cwlout(["sv", "supplemental"], {"type": "array", "items": ["File"]}),
                            cwlout(["svvalidate", "summary"], ["File", "null"]),
                            cwlout("inherit", exclude=[["align_bam"], ["work_bam_plus"],
                                                       ["reference", "snpeff"]])])],
            "bcbio-vc", ["bedtools", "cnvkit", "delly", "duphold", "extract-sv-reads", "gsort",
                         "lumpy-sv;env=python2", "manta;env=python2", "break-point-inspector", "mosdepth", "samtools",
                         "smoove;env=python2", "pysam>=0.13.0",
                         "seq2c", "simple_sv_annotation;env=python2", "survivor", "svtools;env=python2",
                         "svtyper;env=python2",
                         "r=3.5.1", "r-base", "xorg-libxt", "vawk;env=python2"],
            disk={"files": 2.0})]
    sv_batch_inputs = [["analysis"], ["genome_build"],
                       ["work_bam_plus", "disc"], ["work_bam_plus", "sr"],
                       ["config", "algorithm", "background", "cnv_reference"],
                       ["config", "algorithm", "tools_on"],
                       ["config", "algorithm", "tools_off"],
                       ["config", "algorithm", "svprioritize"],
                       ["config", "algorithm", "svvalidate"], ["regions", "sample_callable"],
                       ["genome_resources", "variation", "gc_profile"],
                       ["genome_resources", "variation", "germline_het_pon"],
                       ["genome_resources", "aliases", "snpeff"], ["reference", "snpeff", "genome_build"],
                       ["sv_coverage_rec"]]
    if checkpoints.get("vc"):
        sv_batch_inputs.append(["variants", "samples"])
    steps = [s("calculate_sv_bins", "multi-combined",
               [["align_bam"], ["reference", "fasta", "base"],
                ["metadata", "batch"], ["metadata", "phenotype"],
                ["config", "algorithm", "background", "cnv_reference"],
                ["config", "algorithm", "callable_regions"],
                ["config", "algorithm", "coverage_interval"],
                ["config", "algorithm", "exclude_regions"],
                ["config", "algorithm", "sv_regions"],
                ["config", "algorithm", "variant_regions"],
                ["config", "algorithm", "variant_regions_merged"],
                ["config", "algorithm", "seq2c_bed_ready"],
                ["config", "algorithm", "svcaller"],
                ["depth", "variant_regions", "regions"],
                ["genome_resources", "variation", "lcr"], ["genome_resources", "variation", "polyx"],
                ["genome_resources", "variation", "encode_blacklist"],
                ["genome_resources", "rnaseq", "gene_bed"]],
               [cwlout("sv_bin_rec", "record",
                       fields=[cwlout(["regions", "bins", "target"], ["File", "null"]),
                               cwlout(["regions", "bins", "antitarget"], ["File", "null"]),
                               cwlout(["regions", "bins", "gcannotated"], ["File", "null"]),
                               cwlout(["regions", "bins", "group"], ["string", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["bedtools", "cnvkit"],
               disk={"files": 1.5}, cores=1),
             s("calculate_sv_coverage", "multi-parallel",
               [["sv_bin_rec"]],
               [cwlout("sv_rawcoverage_rec", "record",
                       fields=[cwlout(["depth", "bins", "target"], ["File", "null"]),
                               cwlout(["depth", "bins", "antitarget"], ["File", "null"]),
                               cwlout(["depth", "bins", "seq2c"], ["File", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["mosdepth", "cnvkit", "seq2c"],
               disk={"files": 1.5}),
             s("normalize_sv_coverage", "multi-combined",
               [["sv_rawcoverage_rec"]],
               [cwlout("sv_coverage_rec", "record",
                       fields=[cwlout(["depth", "bins", "normalized"], ["File", "null"]),
                               cwlout(["depth", "bins", "background"], ["File", "null"]),
                               cwlout("inherit")])],
               "bcbio-vc", ["cnvkit"],
               disk={"files": 1.5}),
             s("batch_for_sv", "multi-batch", sv_batch_inputs,
               [cwlout("sv_batch_rec", "record")],
               "bcbio-vc",
               unlist=[["config", "algorithm", "svcaller"]]),
             w("svcall", "multi-parallel", sv, []),
             s("summarize_sv", "multi-combined",
               [["sv_rec"]],
               [cwlout(["sv", "calls"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["sv", "supplemental"], {"type": "array", "items": ["File"]}),
                cwlout(["sv", "prioritize", "tsv"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["sv", "prioritize", "raw"], {"type": "array", "items": ["File", "null"]}),
                cwlout(["svvalidate", "grading_summary"], ["File", "null"]),
                cwlout(["svvalidate", "grading_plots"], {"type": "array", "items": ["File", "null"]})],
               "bcbio-vc", ["bcbio-prioritize"], disk={"files": 1.0}, cores=1)]
    final_outputs = [["sv", "calls"], ["svvalidate", "grading_summary"], ["sv", "prioritize", "tsv"],
                     ["sv", "prioritize", "raw"], ["sv", "supplemental"]]
    return steps, final_outputs
820509
def _call_variants(example_dir, region_bed, data, out_file):
    
    tf_out_file = "%s-tfrecord.gz" % utils.splitext_plus(out_file)[0]
    if not utils.file_exists(tf_out_file):
        with file_transaction(data, tf_out_file) as tx_out_file:
            model = "wes" if strelka2.coverage_interval_from_bed(region_bed) == "targeted" else "wgs"
            cmd = ["dv_call_variants.py", "--cores", dd.get_num_cores(data),
                   "--outfile", tx_out_file, "--examples", example_dir,
                   "--sample", dd.get_sample_name(data), "--model", model]
            do.run(cmd, "DeepVariant call_variants %s" % dd.get_sample_name(data))
    return tf_out_file
301357
def _dup_samples_by_variantcaller(samples, require_bam=True):
    
    samples = [utils.to_single_data(x) for x in samples]
    samples = germline.split_somatic(samples)
    to_process = []
    extras = []
    for data in samples:
        added = False
        for i, add in enumerate(handle_multiple_callers(data, "variantcaller", require_bam=require_bam)):
            added = True
            add = dd.set_variantcaller_order(add, i)
            to_process.append([add])
        if not added:
            data = _handle_precalled(data)
            data = dd.set_variantcaller_order(data, 0)
            extras.append([data])
    return to_process, extras
1072342
def make_V_and_J_mask_mapping(self, genV, genJ):
        
        #construct mapping between allele/gene names and index for custom V_usage_masks
        V_allele_names = [V[0] for V in genV]
        V_mask_mapping = {}
        for v in set([x.split('*')[0] for x in V_allele_names]):
            V_mask_mapping[v] = []
        for v in set(['V'.join((x.split('*')[0]).split('V')[1:]) for x in V_allele_names]):
            V_mask_mapping[v] = []
        for i, v in enumerate(V_allele_names):
            V_mask_mapping[v] = [i]
            V_mask_mapping['V'.join((v.split('*')[0]).split('V')[1:])].append(i)
            V_mask_mapping[v.split('*')[0]].append(i)      
        
        #construct mapping between allele/gene names and index for custom J_usage_masks
        J_allele_names = [J[0] for J in genJ]
        J_mask_mapping = {}
        for j in set([x.split('*')[0] for x in J_allele_names]):
            J_mask_mapping[j] = []
        for j in set(['J'.join((x.split('*')[0]).split('J')[1:]) for x in J_allele_names]):
            J_mask_mapping[j] = []
        for i, j in enumerate(J_allele_names):
            J_mask_mapping[j] = [i]
            J_mask_mapping['J'.join((j.split('*')[0]).split('J')[1:])].append(i)
            J_mask_mapping[j.split('*')[0]].append(i)
            
        self.V_allele_names = V_allele_names
        self.V_mask_mapping = V_mask_mapping
        self.J_allele_names = J_allele_names
        self.J_mask_mapping = J_mask_mapping
295940
def _handle_special_inputs(inputs, variables):
    
    from bcbio import structural
    optional = [["config", "algorithm", "coverage"],
                ["config", "algorithm", "variant_regions"],
                ["config", "algorithm", "sv_regions"],
                ["config", "algorithm", "validate"],
                ["config", "algorithm", "validate_regions"]]
    all_vs = set([get_base_id(v["id"]) for v in variables])
    out = []
    for input in inputs:
        if input == ["reference", "aligner", "indexes"]:
            for v in variables:
                vid = get_base_id(v["id"]).split("__")
                if vid[0] == "reference" and vid[1] in alignment.TOOLS:
                    out.append(vid)
        elif input == ["reference", "snpeff", "genome_build"]:
            found_indexes = False
            for v in variables:
                vid = get_base_id(v["id"]).split("__")
                if vid[0] == "reference" and vid[1] == "snpeff":
                    out.append(vid)
                    found_indexes = True
            assert found_indexes, "Found no snpEff indexes in %s" % [v["id"] for v in variables]
        elif input == ["config", "algorithm", "background", "cnv_reference"]:
            for v in variables:
                vid = get_base_id(v["id"]).split("__")
                if (vid[:4] == ["config", "algorithm", "background", "cnv_reference"] and
                      structural.supports_cnv_reference(vid[4])):
                    out.append(vid)
        elif input in optional:
            if _get_string_vid(input) in all_vs:
                out.append(input)
        else:
            out.append(input)
    return out
973044
def _create_genome_regions(data):
    
    work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "coverage", dd.get_sample_name(data)))
    variant_regions = os.path.join(work_dir, "target-genome.bed")
    with file_transaction(data, variant_regions) as tx_variant_regions:
        with open(tx_variant_regions, "w") as out_handle:
            for c in shared.get_noalt_contigs(data):
                out_handle.write("%s\t%s\t%s\n" % (c.name, 0, c.size))
    return variant_regions
1031806
def gru(name, input, state, kernel_r, kernel_u, kernel_c, bias_r, bias_u, bias_c, new_state, number_of_gates = 2):
    ''' 
    '''
    nn = Build(name)
    inputs = nn.concat(input, state)

    u = nn.sigmoid(nn.mad(inputs, kernel_u, bias_u))
    r = nn.sigmoid(nn.mad(inputs, kernel_r, bias_r))
    r_state = nn.mul(r, state)

    c = nn.tanh(nn.mad(kernel=kernel_c, bias=bias_c,
        x=nn.concat(input, r_state)))

    # new_h = u' * state + (1 - u') * c'
    #       = u' * state + c' - u' * c'

    # u' * state + c'
    nn.add(nn.mul(u, state), c)
    # - u' * c'
    nn.sub(nn._, nn.mul(u, c),
        out=new_state)

    return nn.layers;
103912
def _get_variant_effect(cls, variants, ref_sequence):
        ''''''
        assert len(variants) != 0

        var_types = [x.var_type for x in variants]
        if len(set(var_types)) != 1:
            return 'MULTIPLE', '.', '.'

        var_type = var_types[0]

        assert set([x.ref_name for x in variants]) == set([ref_sequence.id])
        codon_starts = [AssemblyVariants._get_codon_start(0, x.ref_start) for x in variants]
        assert len(set(codon_starts)) == 1
        codon_start = codon_starts[0]
        aa_start = codon_start // 3
        ref_codon = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start:codon_start+3])
        ref_aa = ref_codon.translate()

        if var_type == pymummer.variant.SNP:
            new_codon = list(ref_codon.seq)
            for v in variants:
                new_codon[v.ref_start - codon_start] = v.qry_base
            new_codon = pyfastaq.sequences.Fasta('new', ''.join(new_codon))
            qry_aa = new_codon.translate()

            if ref_aa.seq == qry_aa.seq:
                return ('SYN', '.', aa_start)
            elif qry_aa.seq == '*':
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            else:
                return ('NONSYN', ref_aa.seq + str(aa_start + 1) + qry_aa.seq, aa_start)
        elif var_type in [pymummer.variant.INS, pymummer.variant.DEL]:
            if len(variants) > 1:
                return 'INDELS', '.', aa_start

            var = variants[0]

            if var_type == pymummer.variant.INS:
                new_seq = pyfastaq.sequences.Fasta('seq', var.qry_base)
            else:
                new_seq = pyfastaq.sequences.Fasta('seq', var.ref_base)

            if len(new_seq) % 3 != 0:
                return ('FSHIFT', ref_aa.seq + str(aa_start + 1) + 'fs', aa_start)

            new_seq_aa = new_seq.translate()
            if '*' in new_seq_aa.seq:
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            elif var_type == pymummer.variant.INS:
                ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                aa_after_ins = ref_codon_after_ins.translate()
                return ('INS', ref_aa.seq + str(aa_start + 1) + '_' + aa_after_ins.seq + str(aa_start + 2) + 'ins' + new_seq_aa.seq , aa_start)
            else:
                if len(new_seq) == 3:
                    return ('DEL', ref_aa.seq + str(aa_start + 1) + 'del', aa_start)
                else:
                    assert len(new_seq) % 3 == 0
                    ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                    aa_after_ins = ref_codon_after_ins.translate()
                    return ('DEL', ref_aa.seq + str(aa_start + 1)+ '_' + aa_after_ins.seq + str(aa_start + 2) + 'del', aa_start)

        else:
            return ('UNKNOWN', '.', aa_start)