Structural variant workflow.
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
900380
def generate_random_missense_variants(num_variants=10, max_search=100000, reference="GRCh37"):
    
    variants = []
    for i in range(max_search):
        bases = ["A", "C", "T", "G"]
        random_ref = choice(bases)
        bases.remove(random_ref)
        random_alt = choice(bases)
        random_contig = choice(["1", "2", "3", "4", "5"])
        random_variant = Variant(contig=random_contig, start=randint(1, 1000000),
                                 ref=random_ref, alt=random_alt, ensembl=reference)
        try:
            effects = random_variant.effects()
            for effect in effects:
                if isinstance(effect, Substitution):
                    variants.append(random_variant)
                    break
        except:
            continue
        if len(variants) == num_variants:
            break
    return VariantCollection(variants)
1285435
def _add_variant_sameas_relationships(self, patient_var_map, rs_map):
        
        model = Model(self.graph)
        for patient in patient_var_map:
            for variant_id, variant in patient_var_map[patient].items():
                variant_bnode = self.make_id("{0}".format(variant_id), "_")
                build = variant['build']
                chromosome = variant['chromosome']
                position = variant['position']
                reference_allele = variant['reference_allele']
                variant_allele = variant['variant_allele']
                if build and chromosome and position\
                        and reference_allele and variant_allele:
                    if re.fullmatch(r'[ATCG]', reference_allele)\
                            and re.fullmatch(r'[ATCG]', variant_allele):
                        # variation is snp
                        rs_id = self._get_rs_id(variant, rs_map, 'snp')
                        if rs_id:
                            dbsnp_curie = 'dbSNP:rs{0}'.format(rs_id)
                            model.addSameIndividual(variant_bnode, dbsnp_curie)

                    elif re.fullmatch(r'\-', reference_allele)\
                            or re.fullmatch(r'\-', variant_allele):
                        rs_id = self._get_rs_id(variant, rs_map, 'indel')
                        if rs_id is not None:
                            dbsnp_curie = 'dbSNP:rs{0}'.format(rs_id)
                            model.addSameIndividual(variant_bnode, dbsnp_curie)
                    else:
                        rs_id = self.\
                            _get_rs_id(variant, rs_map, 'indel')
                        if rs_id is not None:
                            dbsnp_curie = 'dbSNP:rs{0}'.format(rs_id)
                            model.addSameIndividual(variant_bnode, dbsnp_curie)
        return
478395
def _compatible_small_variants(data, items):
    
    from bcbio import heterogeneity
    VarFile = collections.namedtuple("VarFile", ["name", "sample", "normal"])
    out = []
    paired = vcfutils.get_paired(items)
    for v in heterogeneity.get_variants(data, include_germline=not paired):
        vrn_file = v["vrn_file"]
        base, ext = utils.splitext_plus(os.path.basename(vrn_file))
        if paired:
            out.append(VarFile(vrn_file, paired.tumor_name, paired.normal_name))
        else:
            out.append(VarFile(vrn_file, dd.get_sample_name(data), None))
    return out
456220
private static List<String> getVariantKeys(String variantKeyPrefix, Collection<String> variants) {

		List<String> variantKeys = new ArrayList<>();
		for (String variant : variants) {
			if (variant == null) {
				variant = "";
			}
			if (variantKeyPrefix == null) {
				variantKeys.add(variant);
			} else {
				variantKeys.add(variantKeyPrefix + variant);
			}
		}

		return variantKeys;
	}
160592
def variants_to_vcf(variants, variant_to_metadata, out=sys.stdout):
    

    # TODO: The variant metadata dictionary (in the `VariantCollection`)
    # contains different data depending on the original input file format (VCF,
    # MAF, CSV).  It's easy to output variants originally in VCF format, but we
    # might want to consider how fields from MAF (for example) map to those in
    # VCF.

    # TODO: The VCF file we output doesn't contain any VCF metadata headers, as
    # the original headers were thrown away when the VCF file was parsed. We
    # may want to keep some of that information and/or infer some of the
    # headers based on the variants themselves. The former is difficult because
    # merge conflicts will inevitably occur; the latter is difficult because
    # the variants themselves don't contain all the information required for
    # these metadata headers (e.g., descriptions).
    #
    # As a side note, adding headers for certain fields will make them parse
    # correctly -- into an integer instead of a string (the default), for
    # example.

    # TODO: If we maintain headers (see above TODO), what should happen if
    # variant sources use different reference genomes?
    #
    # If we don't maintain headers, what should the default reference genome
    # be? This code chose one fairly arbitrarily.

    # TODO: If we end up needing more functions to "build" VCF record fields
    # (see functions below), we may want to abstract away the individual
    # functions and create a mapping from field to format function.

    def get_metadata_field(key, variant, default='.'):
        """Retrieve field from variant metadata dictionary."""
        val = variant_to_metadata[variant].get(key)
        if val is None:
            return default
        return val

    def build_filter_field(variant):
        """Build the filter field from the given variant.

        The `filter` field, as stored in the variants metadata dictionary,
        comes in 3 flavors:

            - empty list: 1+ filters were run and none failed
            - non-empty list: 1+ filters were run and 1+ failed
            - `None`: no filters were run

        This function maps each of these internal representations to their
        corresponding VCF representations.
        """
        filter_metadata = get_metadata_field('filter', variant)
        if type(filter_metadata) == list:
            return 'PASS' if filter_metadata == [] else ';'.join(filter_metadata)
        else:
            # TODO: Can the filter field ever be something other than the 3
            # possibilities described in this function's doc comment?
            bad_value_msg = (
                'Filter metadata field took on unexpected value `{}`. Update '
                'code to account for this value.').format(str(filter_metadata))
            assert filter_metadata == '.', bad_value_msg
            return filter_metadata

    def build_info_field(variant):
        """Build the info field from the given variant.

        Format is `<key>=<val>,...;<key>=<val>,...;<key>=<val>,...`.
        """

        def build_info_pair(key, val):
            """Build key/val pair for info field."""
            # Note: Different from `val == True`, which returns True when `val == 1`.
            if val is True:
                return key

            if type(val) == list:
                val = ','.join(map(str, val))
            else:
                val = str(val)
            return '{}={}'.format(key, val)

        info_dict = get_metadata_field('info', variant, default={})
        if not info_dict:
            return '.'

        return ';'.join(build_info_pair(k, v) for (k, v) in info_dict.items())

    def build_format_field(variant):
        """Build the sample format string from the given variant.

        Each sample column follows this format for the specified variant.
        """
        sample_info = get_metadata_field('sample_info', variant, default={})
        return ':'.join(list(sample_info.values())[0].keys()) if sample_info else '.'

    def build_sample_fields(variant):
        """Build the sample fields for the given variant."""

        def build_sample_field(sample):
            """Build a specific sample's field (i.e., one sample column)."""
            sample_vals = sample_info[sample].values()
            return ':'.join(build_sample_val(val) for val in sample_vals)

        def build_sample_val(sample_val):
            """Build a specific value for a sample (i.e., one value for one column)."""
            if type(sample_val) is list:
                return ','.join(map(str, sample_val))
            elif sample_val is not None:
                return str(sample_val)
            else:
                return '.'

        sample_info = get_metadata_field('sample_info', variant, default={})
        return list(build_sample_field(sample) for sample in sample_info)

    def build_vcf_record(variant, add_sample_info):
        """Return a list of all the variant's VCF record fields."""
        record = [
            str(variant.original_contig),
            str(variant.original_start),
            get_metadata_field('id', variant),
            variant.original_ref,
            variant.original_alt,
            str(get_metadata_field('qual', variant)),
            build_filter_field(variant),
            build_info_field(variant),
        ]
        if add_sample_info:
            record.append(build_format_field(variant))
            record.extend(build_sample_fields(variant))
        return record

    def merge_duplicate_variants():
        """Merge duplicate variants (according to ID) and return *list* of merged, sorted variants.

        Multiple `Variant`s can have the same VCF id (e.g., those variants which
        originally had > 1 alternate base), but we can't output a VCF file with multiple
        records having the same id. This function merges those duplicate variants.
        """

        def construct_id2variants():
            """Construct dictionary which maps variant IDs to `Variant`s."""
            id2variants = defaultdict(list)
            for variant in variants:
                # Note: For variants without IDs (e.g., those that came from
                # MAF files), we assign the missing value.
                variant_id = get_metadata_field('id', variant)
                id2variants[variant_id].append(variant)
            return id2variants

        def merge_variant_list(duplicates):
            """Merge duplicate variant list into one."""

            # TODO: Currently assumes that only alternate bases differ, but we may
            # want or need to merge variants that have the same ID but different
            # contigs, positions, or reference bases. If merging isn't possible (which
            # I think is the case for many situations), it may be easiest to assign
            # the "missing value" to the IDs and move on.

            assert len(duplicates) > 0
            assert all(v.original_contig == duplicates[0].original_contig and
                       v.original_start == duplicates[0].original_start and
                       v.original_ref == duplicates[0].original_ref for v in duplicates)
            import copy
            merged = copy.copy(duplicates[0])
            merged.original_alt = ','.join(duplicate.original_alt for duplicate in duplicates)
            return merged

        id2variants = construct_id2variants()
        variants_no_id = id2variants.pop('.', [])  # don't want to merge variants w/ no id
        merged_variants = list(map(merge_variant_list, id2variants.values())) + variants_no_id

        # groups variants by contig; relative ordering of contigs doesn't matter
        return sorted(
            merged_variants,
            key=lambda v: (str(v.original_contig), str(v.original_start)))

    def get_sample_names():
        """Return the sample names for all variants."""

        # TODO: For now, ensures that every variant has the same samples. If they didn't,
        # we'd have to use the missing value for VCF records with no data on a given
        # sample.  In and of itself, that isn't a problem until the format field contains
        # GT (genotype), which is required for all samples if present in the format field.
        #
        # A couple of ways to handle this:
        #
        #   1) Ignore this requirement of VCF files. I'd assume this would still be
        #   compatible with varcode, but may not be with other VCF parsers.
        #   2) Remove the GT field from those records where this rule would be violated.
        #   This is sad because we lose information.
        #
        # It's also important to note that if/when we combine all samples, we'll have to
        # add a `.` (missing value) for those samples in which a variant has no data.
        # Check the VCF spec to make sure this is valid; if not, we may have to write
        # `.:.:.:.` -- one `.` for each field in the format string.
        #
        # Moreover, we'll want to note whether we care about maintaining the relative
        # ordering of the sample names in the original VCF files. This probably isn't
        # important, but the code below does not do this (it effectively alphabetizes the
        # sample names because of the use of `set`).

        sample_names = set()
        for variant in variants:
            sample_info = get_metadata_field('sample_info', variant, default={})
            addl_sample_names = set(sample_info.keys())

            # Ensures all variants have the same samples.
            if sample_names and sample_names != addl_sample_names:
                return []

            sample_names.update(addl_sample_names)

        return list(sample_names)

    headers = ['#CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']
    sample_names = get_sample_names()
    if sample_names:
        headers += ['FORMAT'] + sample_names

    unique_variants_list = merge_duplicate_variants()

    # usually we have just one reference genome for the whole variant collection
    # but if variants from multiple sources have been merged then we might
    # not be able to write out a VCF since the individual variants may be using
    # different coordinate systems
    genome_names = list({v.ensembl.reference_name for v in unique_variants_list})
    if len(genome_names) > 1:
        raise ValueError(
            "Cannot create VCF for variants with multiple reference genomes: %s" % (
                genome_names,))
    genome_name = genome_names[0]
    print('##fileformat=VCFv4.2', file=out)
    print('##reference=%s' % genome_name, file=out)
    print('\t'.join(headers), file=out)
    for variant in unique_variants_list:
        record = build_vcf_record(variant, add_sample_info=bool(sample_names))
        print('\t'.join(record), file=out)
1055887
public void addPostProcessVariant(Map<String, VariantSet> variants) {

		if (!searchingPostProcessorVariants) {
			throw new BundlingProcessException(
					"You are not allowed to define post process variants if we are not searching for post processor variants.");
		}
		postProcessVariants = VariantUtils.concatVariants(postProcessVariants, variants);
	}
593618
function variants(sequence, cb) {
        if (sequence.length === 2) {
            cb(sequence);
            cb(sequence[1] + sequence[0]);
        } else if (sequence.length > 2) {
            for(var i = 0; i < sequence.length; i++) {
                (function(sequence, ch) {
                    variants(sequence, function(v) {
                        cb(ch + v);
                    });
                })(sequence.substring(0, i) + sequence.substring(i+1), sequence[i]);
            }
        } else {
            cb(sequence);
        }
    }
1208848
private void mergeVariants(Variant current, List<Pair<Variant, List<AlternateCoordinate>>> varToAlts) {
        StudyEntry currentStudy = getStudy(current);

        // Build ALT index
        List<AlternateCoordinate> altList = buildAltsList(current, varToAlts);
        // Create a list of alternates as Integers to speed up the creation of the VariantAlternateRearranger
        List<Integer> altListHash = alternatesToHash(altList);
//        Map<AlternateCoordinate, Integer> altIdx = index(altList);

        // Check if the number of secondary alternates has increased. If so, rearrange the current study (files + samples)
        VariantAlternateRearranger currentStudyRearranger = null;
        if (altList.size() - 1 != currentStudy.getSecondaryAlternates().size()) {
            currentStudyRearranger = new VariantAlternateRearranger(alternatesToHash(buildAltList(current)), altListHash);
        }

        // Update SecALt list
        currentStudy.setSecondaryAlternates(altList.subList(1, altList.size()));

        // Find new formats
        final Map<String, Integer> newFormatPositions;
        final List<String> newFormat;
        if (expectedFormatsPosition.isEmpty()) {
            // Find all formats
            Set<String> allFormats = varToAlts.stream()
                    .map(Pair::getKey)
                    .map(this::getStudy)
                    .map(StudyEntry::getFormat)
                    .flatMap(Collection::stream)
                    .collect(Collectors.toSet());
            newFormatPositions = new HashMap<>(currentStudy.getFormatPositions());
//            newFormatPositions.putIfAbsent(getFilterKey(), newFormatPositions.size());
            for (String format : allFormats) {
                newFormatPositions.putIfAbsent(format, newFormatPositions.size());
            }
            newFormat = Arrays.asList(new String[newFormatPositions.size()]);
            newFormatPositions.forEach((format, formatIdx) -> {
                newFormat.set(formatIdx, format);
            });
        } else {
            // TODO: Try to avoid this copy
            newFormatPositions = new HashMap<>(expectedFormatsPosition);
            newFormat = new ArrayList<>(expectedFormatsPosition.keySet());
        }
        Map<String, Integer> extraFormats = newFormatPositions.entrySet()
                .stream()
                .filter(e -> !e.getKey().equals(getGtKey()) && !e.getKey().equals(getFilterKey()))
                .collect(Collectors.toMap(Map.Entry::getKey, Map.Entry::getValue));

        // Find new samples
        LinkedHashMap<String, Integer> newSamplesPosition;
        if (expectedSamplesPosition.isEmpty()) {
            newSamplesPosition = new LinkedHashMap<>(currentStudy.getSamplesPosition());
            varToAlts.stream()
                    .map(Pair::getKey)
                    .map(this::getStudy)
                    .map(StudyEntry::getOrderedSamplesName)
                    .flatMap(List::stream)
                    .forEach(sample -> newSamplesPosition.putIfAbsent(sample, newSamplesPosition.size()));
        } else {
            // TODO: Try to avoid this copy
            newSamplesPosition = new LinkedHashMap<>(expectedSamplesPosition);
        }

        Integer newGtIdx = newFormatPositions.get(getGtKey());
        Integer newFilterIdx = newFormatPositions.get(getFilterKey());

        // Create new Samples data
        List<List<String>> newSamplesData = newSamplesData(newSamplesPosition.size(), newFormatPositions);
        boolean[] alreadyMergedSamples = new boolean[newSamplesPosition.size()];
        // Copy current samples data into new samples data
        List<List<String>> currentSamplesData = currentStudy.getSamplesData();
        // If FT is required, but missing in the input variant, get FILTER from the File Attributes.
        String currentFilterValue = null;
        if (!currentStudy.getFormatPositions().containsKey(getFilterKey()) && newFormatPositions.containsKey(getFilterKey())) {
            currentFilterValue = currentStudy.getFiles().isEmpty() ? getDefaultValue(getFilterKey())
                    : currentStudy.getFiles().get(0).getAttributes().getOrDefault(getAnnotationFilterKey(), getDefaultValue(getFilterKey()));
        }
        int currentSampleIdx = 0;
        for (String sample : currentStudy.getOrderedSamplesName()) {
            List<String> currentSampleData = currentSamplesData.get(currentSampleIdx);
            Integer newSampleIdx = newSamplesPosition.get(sample);
            alreadyMergedSamples[newSampleIdx] = true;
            List<String> newSampleData = newSamplesData.get(newSampleIdx);
            int formatIdx = 0;
            int ploidy = 2; // Default ploidy. Required for the rearranger
            for (String format : currentStudy.getFormat()) {
                Integer newFormatIdx = newFormatPositions.get(format);
                if (newFormatIdx != null) {
                    if (currentSampleData.size() > formatIdx) {
                        String data = currentSampleData.get(formatIdx);
                        if (currentStudyRearranger != null) {
                            if (format.equals(VCFConstants.GENOTYPE_KEY)) {
                                Genotype genotype = new Genotype(data);
                                ploidy = genotype.getPloidy();
                                data = currentStudyRearranger.rearrangeGenotype(genotype).toString();
                            } else {
                                data = currentStudyRearranger.rearrange(format, data, ploidy);
                            }
                        }
                        newSampleData.set(newFormatIdx, data);
                    } else {
                        newSampleData.set(newFormatIdx, getDefaultValue(format));
                    }
                }
                formatIdx++;
            }
            // If FT is required, but missing in the input variant, add FILTER from the File Attributes.
            if (newFilterIdx != null && currentFilterValue != null) {
                newSampleData.set(newFilterIdx, currentFilterValue);
            }
            currentSampleIdx++;
        }
        if (currentStudyRearranger != null) {
            for (FileEntry file : currentStudy.getFiles()) {
                file.getAttributes().replaceAll(currentStudyRearranger::rearrange);
            }
        }

        for (Pair<Variant, List<AlternateCoordinate>> e : varToAlts) {
            Variant other = e.getKey();
            List<AlternateCoordinate> otherAlternates = e.getValue();


//            Map<Integer, AlternateCoordinate> otherAltIdx = index(alternates).entrySet().stream()
//                    .collect(Collectors.toMap(Map.Entry::getValue, Map.Entry::getKey));
            final StudyEntry otherStudy = getStudy(other);
            Map<String, Integer> otherStudyFormatPositions = otherStudy.getFormatPositions();
            checkForDuplicates(current, other, currentStudy, otherStudy, otherAlternates);

            VariantAlternateRearranger rearranger;
            Map<String, String> rearrangedGenotypesCache = new HashMap<>(5);
            // It may happen that the new list of alternates does not contains some of the other alternates.
            // In that case, use the rearranger
            if (altList.size() == 1 && altList.equals(otherAlternates)) {
                rearranger = null;
            } else {
                rearranger = new VariantAlternateRearranger(
                        alternatesToHash(otherAlternates),
                        altListHash, rearrangerConf);
            }
            // Add GT data for each sample to current Variant

            List<String> otherOrderedSamplesName = otherStudy.getOrderedSamplesName();
            for (int sampleIdx = 0; sampleIdx < otherOrderedSamplesName.size(); sampleIdx++) {
                String sampleName = otherOrderedSamplesName.get(sampleIdx);
                List<String> otherSampleData = otherStudy.getSamplesData().get(sampleIdx);
                Integer newSampleIdx = newSamplesPosition.get(sampleName);
                List<String> newSampleData = newSamplesData.get(newSampleIdx);

                boolean alreadyMergedSample = alreadyMergedSamples[newSampleIdx];
                alreadyMergedSamples[newSampleIdx] = true;

                // GT data
                int ploidy = -1;
                boolean isGtUpdated = false;
                List<Integer> updatedGtPositions = Collections.emptyList();
                if (newGtIdx != null) {
                    String gt = otherSampleData.get(otherStudyFormatPositions.get(getGtKey()));
                    if (StringUtils.isBlank(gt)) {
                        throw new IllegalStateException(String.format(
                                "No GT [%s] found for sample %s in \nVariant: %s\nOtherSe:%s\nOtherSp:%s",
                                getGtKey(), sampleName, other.getImpl(), otherStudy.getSamplesData(),
                                otherStudy.getSamplesPosition()));
                    }
                    // Use a cache with rearranged genotypes to reuse rearranged genotypes
                    String updatedGt = rearrangedGenotypesCache.get(gt);
                    if (updatedGt == null) {
                        Genotype genotype;
                        if (rearranger != null) {
                            genotype = rearranger.rearrangeGenotype(new Genotype(gt));
                        } else {
                            genotype = new Genotype(gt);
                        }
                        if (!genotype.isPhased()) {
                            genotype.normalizeAllelesIdx();
                        }
                        if (collapseDeletions) {
                            int[] allelesIdx = genotype.getAllelesIdx();
                            for (int i = 0; i < allelesIdx.length; i++) {
                                if (allelesIdx[i] < 0) {
                                    allelesIdx[i] = 0; // change to '0' for 'missing' reference (missing because change to '0' GT)
                                }
                            }
                        }

//                    updatedGt = updateGT(gt, altIdx, otherAltIdx);
                        updatedGt = genotype.toString();
                        rearrangedGenotypesCache.put(gt, updatedGt);
                        ploidy = genotype.getPloidy();
                    } else {
                        ploidy = Genotype.getPloidy(gt);
                    }

                    if (alreadyMergedSample) {
                        String currGT = newSampleData.get(newGtIdx);
                        List<String> gtlst;
                        if (currGT.contains(",")) {
                            gtlst = new ArrayList<>(Arrays.asList(currGT.split(",")));
                        } else {
                            gtlst = new ArrayList<>(1);
                            gtlst.add(currGT);
                        }
                        if (!gtlst.contains(updatedGt)) {
                            gtlst.add(updatedGt);
                            updatedGtPositions = collapseGT(gtlst);
                            updatedGt = StringUtils.join(updatedGtPositions.stream()
                                    .map(gtlst::get).collect(Collectors.toList()), ',');
                            isGtUpdated = true;
                        }
                    }
                    newSampleData.set(newGtIdx, updatedGt);
                }

                // Filter
                if (newFilterIdx != null) {
                    Integer otherFilterIdx = otherStudyFormatPositions.get(getFilterKey());
                    String filter;
                    if (otherFilterIdx != null) {
                        filter = otherSampleData.get(otherFilterIdx);
                    } else {
                        filter = otherStudy.getFiles().get(0).getAttributes()
                                .getOrDefault(StudyEntry.FILTER, getDefaultValue(getFilterKey()));
                    }

                    if (alreadyMergedSample && isGtUpdated) {
                        String currFilter = newSampleData.get(newFilterIdx);
                        if (currFilter != null && !currFilter.equals(getDefaultValue(getFilterKey()))) {
                            List<String> filterLst = new ArrayList<>(Arrays.asList(currFilter.split(",")));
                            filterLst.add(filter);
                            if (/*Objects.isNull(sampleToGt)*/updatedGtPositions.isEmpty()) {
                                filter = StringUtils.join(filterLst, ',');
                            } else {
                                filter = StringUtils.join(updatedGtPositions.stream()
                                        .map(filterLst::get).collect(Collectors.toList()), ',');
                            }
                        }
                    }
                    newSampleData.set(newFilterIdx, filter);
                }

                // Additional data
                for (Map.Entry<String, Integer> entry : extraFormats.entrySet()) {
                    Integer idx = otherStudyFormatPositions.get(entry.getKey());
                    String data = idx == null || idx >= otherSampleData.size() ? getDefaultValue(entry.getKey()) : otherSampleData.get(idx);
                    if (StringUtils.isNotEmpty(data)) {
                        if (rearranger != null) {
                            data = rearranger.rearrange(entry.getKey(), data, ploidy);
                        }
                        String old = newSampleData.set(entry.getValue(), data);
//                        if (old != null) {
//                            throw new IllegalStateException("TODO - merge additional formats!!!");
//                        }
                    }
                }

            }
            mergeFile(current, other, rearranger, currentStudy, otherStudy);
        }
        currentStudy.setSamplesData(newSamplesData);
        currentStudy.setSortedSamplesPosition(newSamplesPosition);
        currentStudy.setFormat(newFormat);
    }
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)
533208
def reads_supporting_variants(variants, samfile, **kwargs):
    
    for variant, allele_reads in reads_overlapping_variants(
            variants=variants,
            samfile=samfile,
            **kwargs):
        yield variant, filter_non_alt_reads_for_variant(variant, allele_reads)
89495
def get_variants_in_region(self, chrom, start, end):
        """"""
        bim = self.bim.loc[
            (self.bim["chrom"] == CHROM_STR_TO_INT[chrom]) &
            (start <= self.bim["pos"]) &
            (self.bim["pos"] <= end)
        ]
        for i, g in enumerate(self.bed.iter_geno_marker(bim.index)):
            info = bim.iloc[i, :]
            name, geno = g
            yield Genotypes(
                Variant(info.name, CHROM_INT_TO_STR[info.chrom],
                        info.pos, [info.a1, info.a2]),
                self._normalize_missing(geno),
                reference=info.a2,
                coded=info.a1,
                multiallelic=info.multiallelic
            )
1231664
def iter_variants(self):
        """"""
        for variant in self._bgen.iter_variant_info():
            yield Variant(
                variant.name,
                CHROM_STR_ENCODE.get(variant.chrom, variant.chrom),
                variant.pos, [variant.a1, variant.a2],
            )
579009
def variation(fastas, sams, s2b, annots, qual, min_cov, min_per, no_sub, bp_table, print_con):
    
    genomes = {}

    # save base-pair data structures for each genome-sample pair
    for sam in sams:
        genomes, s2b = save_refs(sam, fastas, genomes, s2b)

    # get base-pair frequencies for each genome-sample pair
    for sam in sams:
        genomes = get_bp_stats(sam, genomes, s2b, qual)

    # calculate base-pair frequencies
    genomes = calc_frequencies(genomes, bp_table, min_cov, min_per)

    # print consensus genome
    if print_con is True:
        print('# saving consensus sequences', file = sys.stderr)
        genomes = print_consensus(genomes)

    if no_sub is True:
        return genomes

    # calculate substitution rates
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample in samples:
#                print(genomes.keys())
#                print(genomes[genome].keys())
#                print(genomes[genome][contig].keys())
#                print(annots.keys())
                genomes[genome][contig][sample]['sub_rates'] = \
                    sub_rates(genomes[genome][contig][sample]['bp_stats'], annots[contig])

    return genomes
901150
def version_variants(version):
    """"""

    result = [version]

    # Add trailing ".0" as needed to ensure that we have at least
    # major.minor.patch
    parts = version.split(".")
    while len(parts) < 3:
        parts.append("0")
        result.append(".".join(parts))

    return result
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)
849701
def _variant_checkpoints(samples):
    
    checkpoints = {}
    checkpoints["vc"] = any([dd.get_variantcaller(d) or d.get("vrn_file") for d in samples])
    checkpoints["sv"] = any([dd.get_svcaller(d) for d in samples])
    checkpoints["jointvc"] = any([(dd.get_jointcaller(d) or "gvcf" in dd.get_tools_on(d))
                                  for d in samples])
    checkpoints["hla"] = any([dd.get_hlacaller(d) for d in samples])
    checkpoints["align"] = any([(dd.get_aligner(d) or dd.get_bam_clean(d)) for d in samples])
    checkpoints["align_split"] = not all([(dd.get_align_split_size(d) is False or
                                           not dd.get_aligner(d))
                                          for d in samples])
    checkpoints["archive"] = any([dd.get_archive(d) for d in samples])
    checkpoints["umi"] = any([dd.get_umi_consensus(d) for d in samples])
    checkpoints["ensemble"] = any([dd.get_ensemble(d) for d in samples])
    checkpoints["cancer"] = any(dd.get_phenotype(d) in ["tumor"] for d in samples)
    return checkpoints