Generate a VariantCollection with random variants that overlap
    at least one complete coding transcript.
763976
def random_variants(
        count,
        genome_name="GRCh38",
        deletions=True,
        insertions=True,
        random_seed=None):
    
    rng = random.Random(random_seed)
    ensembl = genome_for_reference_name(genome_name)

    if ensembl in _transcript_ids_cache:
        transcript_ids = _transcript_ids_cache[ensembl]
    else:
        transcript_ids = ensembl.transcript_ids()
        _transcript_ids_cache[ensembl] = transcript_ids

    variants = []

    # we should finish way before this loop is over but just in case
    # something is wrong with PyEnsembl we want to avoid an infinite loop
    for _ in range(count * 100):
        if len(variants) < count:
            transcript_id = rng.choice(transcript_ids)
            transcript = ensembl.transcript_by_id(transcript_id)

            if not transcript.complete:
                continue

            exon = rng.choice(transcript.exons)
            base1_genomic_position = rng.randint(exon.start, exon.end)
            transcript_offset = transcript.spliced_offset(base1_genomic_position)
            seq = transcript.sequence

            ref = str(seq[transcript_offset])
            if transcript.on_backward_strand:
                ref = reverse_complement(ref)

            alt_nucleotides = [x for x in STANDARD_NUCLEOTIDES if x != ref]

            if insertions:
                nucleotide_pairs = [
                    x + y
                    for x in STANDARD_NUCLEOTIDES
                    for y in STANDARD_NUCLEOTIDES
                ]
                alt_nucleotides.extend(nucleotide_pairs)
            if deletions:
                alt_nucleotides.append("")
            alt = rng.choice(alt_nucleotides)
            variant = Variant(
                transcript.contig,
                base1_genomic_position,
                ref=ref,
                alt=alt,
                ensembl=ensembl)
            variants.append(variant)
        else:
            return VariantCollection(variants)
    raise ValueError(
        ("Unable to generate %d random variants, "
         "there may be a problem with PyEnsembl") % count)
1068187
function generateFeatures() {
  var collection = {
    type: 'FeatureCollection',
    features: []
  };

  for (var i = 0; i < options.number; i++) {
    collection.features.push(generateFeature(randomType()));
  }
  return collection;
}
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
1085960
def export(ctx, outfile, variant_type):
    
    adapter = ctx.obj['adapter']
    version = ctx.obj['version']
    
    LOG.info("Export the variants from {0}".format(adapter))
    nr_cases = 0

    is_sv = variant_type == 'sv'
    existing_chromosomes = set(adapter.get_chromosomes(sv=is_sv))
    
    ordered_chromosomes = []
    for chrom in CHROMOSOME_ORDER:
        if chrom in existing_chromosomes:
            ordered_chromosomes.append(chrom)
            existing_chromosomes.remove(chrom)
    for chrom in existing_chromosomes:
        ordered_chromosomes.append(chrom)
    
    nr_cases = adapter.cases().count()
    LOG.info("Found {0} cases in database".format(nr_cases))

    head = HeaderParser()
    head.add_fileformat("VCFv4.3")
    head.add_meta_line("NrCases", nr_cases)
    head.add_info("Obs", '1', 'Integer', "The number of observations for the variant")
    head.add_info("Hom", '1', 'Integer', "The number of observed homozygotes")
    head.add_info("Hem", '1', 'Integer', "The number of observed hemizygotes")
    head.add_version_tracking("loqusdb", version, datetime.now().strftime("%Y-%m-%d %H:%M"))
    
    if variant_type == 'sv':
        head.add_info("END", '1', 'Integer', "End position of the variant")
        head.add_info("SVTYPE", '1', 'String', "Type of structural variant")
        head.add_info("SVLEN", '1', 'Integer', "Length of structural variant")
        
        
    for chrom in ordered_chromosomes:
        length = adapter.get_max_position(chrom)
        head.add_contig(contig_id=chrom, length=str(length))

    print_headers(head, outfile=outfile)
    
    for chrom in ordered_chromosomes:
        if variant_type == 'snv':
            LOG.info("Collecting all SNV variants")
            variants = adapter.get_variants(chromosome=chrom)
        else:
            LOG.info("Collecting all SV variants")
            variants = adapter.get_sv_variants(chromosome=chrom)
        LOG.info("{} variants found".format(variants.count()))
        for variant in variants:
            variant_line = format_variant(variant, variant_type=variant_type)
            # chrom = variant['chrom']
            # pos = variant['start']
            # ref = variant['ref']
            # alt = variant['alt']
            # observations = variant['observations']
            # homozygotes = variant['homozygote']
            # hemizygotes = variant['hemizygote']
            # info = "Obs={0}".format(observations)
            # if homozygotes:
            #     info += ";Hom={0}".format(homozygotes)
            # if hemizygotes:
            #     info += ";Hem={0}".format(hemizygotes)
            # variant_line = "{0}\t{1}\t.\t{2}\t{3}\t.\t.\t{4}\n".format(
            #     chrom, pos, ref, alt, info)
            print_variant(variant_line=variant_line, outfile=outfile)
441367
def make_separate_indels_and_one_alt_with_all_snps_no_combinations(self, ref_seq):
        ''''''
        final_start_position = min([x.POS for x in self.vcf_records])
        final_end_position = max([x.ref_end_pos() for x in self.vcf_records])
        snps = []
        new_vcf_records = []

        for record in self.vcf_records:
            if record.is_snp():
                snps.append(copy.copy(record))
            else:
                new_record = copy.copy(record)
                new_record.add_flanking_seqs(ref_seq, final_start_position, final_end_position)
                new_vcf_records.append(new_record)

        if len(snps):
            new_record = copy.copy(snps[0])
            for snp in snps[1:]:
                merged = new_record.merge(snp, ref_seq)
                if merged is not None:
                    new_record = merged
            new_record.add_flanking_seqs(ref_seq, final_start_position, final_end_position)
            new_vcf_records.append(new_record)

        alts = ','.join(sorted(list(set([x.ALT[0] for x in new_vcf_records]))))
        new_record = vcf_record.VcfRecord('\t'.join([self.vcf_records[0].CHROM, str(final_start_position + 1), '.', new_vcf_records[0].REF, alts, '.', 'PASS', '.']))

        return new_record
9577
def _get_files_variantcall(sample):
    
    out = []
    algorithm = sample["config"]["algorithm"]
    out = _maybe_add_summary(algorithm, sample, out)
    out = _maybe_add_alignment(algorithm, sample, out)
    out = _maybe_add_callable(sample, out)
    out = _maybe_add_disambiguate(algorithm, sample, out)
    out = _maybe_add_variant_file(algorithm, sample, out)
    out = _maybe_add_sv(algorithm, sample, out)
    out = _maybe_add_hla(algorithm, sample, out)
    out = _maybe_add_heterogeneity(algorithm, sample, out)

    out = _maybe_add_validate(algorithm, sample, out)
    return _add_meta(out, sample)
266672
def add_variant(self, chrom, pos, rs_id, ref, alt, qual, filt, info, form=None, genotypes=[]):
        
        variant_info = [chrom, pos, rs_id, ref, alt, qual, filt, info]
        if form:
            variant_info.append(form)
        for individual in genotypes:
            variant_info.append(individual)
        
        variant_line = '\t'.join(variant_info)
        variant = format_variant(
            line = variant_line, 
            header_parser = self.metadata, 
            check_info = self.check_info
        )
        
        if not (self.split_variants and len(variant['ALT'].split(',')) > 1):
            self.variants.append(variant)
            
        # If multiple alternative and split_variants we must split the variant                 
        else:
            for splitted_variant in split_variants(
                                                    variant_dict=variant, 
                                                    header_parser=self.metadata, 
                                                    allele_symbol=self.allele_symbol):
                self.variants.append(splitted_variant)
10599
def population_variant_regions(items, merged=False):
    
    def _get_variant_regions(data):
        out = dd.get_variant_regions(data) or dd.get_sample_callable(data)
        # Only need to merge for variant region inputs, not callable BED regions which don't overlap
        if merged and dd.get_variant_regions(data):
            merged_out = dd.get_variant_regions_merged(data)
            if merged_out:
                out = merged_out
            else:
                out = merge_overlaps(out, data)
        return out
    import pybedtools
    if len(items) == 1:
        return _get_variant_regions(items[0])
    else:
        paired = vcfutils.get_paired(items)
        if paired:
            return _get_variant_regions(paired.tumor_data)
        else:
            vrs = []
            for data in items:
                vr_bed = _get_variant_regions(data)
                if vr_bed:
                    vrs.append((pybedtools.BedTool(vr_bed).total_coverage(), vr_bed))
            vrs.sort(reverse=True)
            if vrs:
                return vrs[0][1]
984951
async def variant(self, elem=None, elem_type=None, params=None, wrapped=None):
        
        elem_type = elem_type if elem_type else elem.__class__
        if hasattr(elem_type, "serialize_archive"):
            elem = elem_type() if elem is None else elem
            return await elem.serialize_archive(
                self, elem=elem, elem_type=elem_type, params=params
            )

        if self.writing:
            return await self._dump_variant(
                self.iobj,
                elem=elem,
                elem_type=elem_type if elem_type else elem.__class__,
                params=params,
            )
        else:
            return await self._load_variant(
                self.iobj,
                elem_type=elem_type if elem_type else elem.__class__,
                params=params,
                elem=elem,
                wrapped=wrapped,
            )
92301
def variantAnnotationsGenerator(self, request):
        
        compoundId = datamodel.VariantAnnotationSetCompoundId.parse(
            request.variant_annotation_set_id)
        dataset = self.getDataRepository().getDataset(compoundId.dataset_id)
        variantSet = dataset.getVariantSet(compoundId.variant_set_id)
        variantAnnotationSet = variantSet.getVariantAnnotationSet(
            request.variant_annotation_set_id)
        iterator = paging.VariantAnnotationsIntervalIterator(
            request, variantAnnotationSet)
        return iterator
175408
def _get_variant_file(x, key, suffix="", sample=None, ignore_do_upload=False):
    
    out = []
    fname = utils.get_in(x, key)
    upload_key = list(key)
    upload_key[-1] = "do_upload"
    do_upload = tz.get_in(tuple(upload_key), x, True)
    if fname and (ignore_do_upload or do_upload):
        if fname.endswith(".vcf.gz"):
            out.append({"path": fname,
                        "type": "vcf.gz",
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
            if utils.file_exists(fname + ".tbi"):
                out.append({"path": fname + ".tbi",
                            "type": "vcf.gz.tbi",
                            "index": True,
                            "ext": "%s%s" % (x["variantcaller"], suffix),
                            "variantcaller": x["variantcaller"]})
        elif fname.endswith((".vcf", ".bed", ".bedpe", ".bedgraph", ".cnr", ".cns", ".cnn", ".txt", ".tsv")):
            ftype = utils.splitext_plus(fname)[-1][1:]
            if ftype == "txt":
                extended_ftype = fname.split("-")[-1]
                if "/" not in extended_ftype:
                    ftype = extended_ftype
            out.append({"path": fname,
                        "type": ftype,
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
    if sample:
        out_sample = []
        for x in out:
            x["sample"] = sample
            out_sample.append(x)
        return out_sample
    else:
        return out
1362225
def variants(ctx, snpeff):
    """"""
    head = ctx.parent.head
    vcf_handle = ctx.parent.handle
    outfile = ctx.parent.outfile
    silent = ctx.parent.silent
    
    print_headers(head, outfile=outfile, silent=silent)
    
    for line in vcf_handle:
        print_variant(variant_line=line, outfile=outfile, silent=silent)
        if snpeff:
            variant_dict =  get_variant_dict(
                variant_line = line,
                header_line = head.header
            )
            #Create a info dict:
            info_dict = get_info_dict(
                info_line = variant_dict['INFO']
            )
            snpeff_string = info_dict.get('ANN')

            if snpeff_string:
                #Get the snpeff annotations
                snpeff_info = get_snpeff_info(
                    snpeff_string = snpeff_string,
                    snpeff_header = head.snpeff_columns
                )
1095406
def trim_variant_sequences(variant_sequences, min_variant_sequence_coverage):
    
    n_total = len(variant_sequences)
    trimmed_variant_sequences = [
        variant_sequence.trim_by_coverage(min_variant_sequence_coverage)
        for variant_sequence in variant_sequences
    ]
    collapsed_variant_sequences = collapse_substrings(trimmed_variant_sequences)
    n_after_trimming = len(collapsed_variant_sequences)
    logger.info(
        "Kept %d/%d variant sequences after read coverage trimming to >=%dx",
        n_after_trimming,
        n_total,
        min_variant_sequence_coverage)
    return collapsed_variant_sequences
774878
def finalize_vcf(in_file, variantcaller, items):
    
    out_file = "%s-annotated%s" % utils.splitext_plus(in_file)
    if not utils.file_uptodate(out_file, in_file):
        header_cl = _add_vcf_header_sample_cl(in_file, items, out_file)
        contig_cl = _add_contig_cl(in_file, items, out_file)
        cls = [x for x in (contig_cl, header_cl) if x]
        if cls:
            post_cl = " | ".join(cls) + " | "
        else:
            post_cl = None
        dbsnp_file = tz.get_in(("genome_resources", "variation", "dbsnp"), items[0])
        if dbsnp_file:
            out_file = _add_dbsnp(in_file, dbsnp_file, items[0], out_file, post_cl)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, items[0]["config"])
    else:
        return in_file
508087
def _variants_fills(fields, fills, info_types):
    """"""
    if fills is None:
        # no fills specified by user
        fills = dict()
    for f, vcf_type in zip(fields, info_types):
        if f == 'FILTER':
            fills[f] = False
        elif f not in fills:
            if f in config.STANDARD_VARIANT_FIELDS:
                fills[f] = config.DEFAULT_VARIANT_FILL[f]
            else:
                fills[f] = config.DEFAULT_FILL_MAP[vcf_type]
    # convert to tuple for zipping with fields
    fills = tuple(fills[f] for f in fields)
    return fills
62815
def variantcall_sample(data, region=None, align_bams=None, out_file=None):
    
    if out_file is None or not os.path.exists(out_file) or not os.path.lexists(out_file):
        utils.safe_makedir(os.path.dirname(out_file))
        ref_file = dd.get_ref_file(data)
        config = data["config"]
        caller_fns = get_variantcallers()
        caller_fn = caller_fns[config["algorithm"].get("variantcaller")]
        if len(align_bams) == 1:
            items = [data]
        else:
            items = multi.get_orig_items(data)
            assert len(items) == len(align_bams)
        assoc_files = tz.get_in(("genome_resources", "variation"), data, {})
        if not assoc_files: assoc_files = {}
        for bam_file in align_bams:
            bam.index(bam_file, data["config"], check_timestamp=False)
        out_file = caller_fn(align_bams, items, ref_file, assoc_files, region, out_file)
    if region:
        data["region"] = region
    data["vrn_file"] = out_file
    return [data]
52131
def _variant_filtration(in_file, ref_file, vrn_files, data, filter_type,
                        hard_filter_fn):
    
    # Algorithms multiplied by number of input files to check for large enough sample sizes
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if (not config_utils.use_vqsr(algs, in_file) or
          _already_cutoff_filtered(in_file, filter_type)):
        logger.info("Skipping VQSR, using cutoff-based filers: we don't have whole genome input data")
        return hard_filter_fn(in_file, data)
    elif not _have_training_data(vrn_files):
        logger.info("Skipping VQSR, using cutoff-based filers: genome build does not have sufficient training data")
        return hard_filter_fn(in_file, data)
    else:
        sensitivities = {"INDEL": "98.0", "SNP": "99.97"}
        recal_file, tranches_file = _run_vqsr(in_file, ref_file, vrn_files,
                                              sensitivities[filter_type], filter_type, data)
        if recal_file is None:  # VQSR failed
            logger.info("VQSR failed due to lack of training data. Using cutoff-based soft filtering.")
            return hard_filter_fn(in_file, data)
        else:
            return _apply_vqsr(in_file, ref_file, recal_file, tranches_file,
                               sensitivities[filter_type], filter_type, data)