Multi-omics reveal phosphatidic acid phosphatases modify Niemann-Pick type C disease severity.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches

[1] § STAR★Methods › Method details › Human genome sequence analysis ↔ deepvariant/postprocess_variants.py, lines 611–646 · score 0.71 · Phred scaled, quality score, variant genotype, sequencing, Predictor, VCF
[2] § STAR★Methods › Quantification and statistical analysis ↔ deepvariant/methylation_aware_phasing.cc, lines 231–319 · score 0.52 · Wilcoxon rank sum
[3] § STAR★Methods › Quantification and statistical analysis ↔ deepvariant/methylation_aware_phasing.h, lines 46–108 · score 0.51 · Wilcoxon rank sum

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 2,378 lines · 80 KB · BSD-3-Clause · 1 match

# Copyright 2017 Google LLC.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""Postprocess output from call_variants to produce a VCF file."""
# TODO: Add type annotations to this module
import collections
import functools
import itertools
import json
import os
import tempfile
import time
from typing import Iterable, Iterator, Sequence
from absl import flags
from absl import logging
from google.protobuf import json_format
import numpy as np
import pysam
import tensorflow as tf
from deepvariant import calling_regions_utils
from deepvariant import dv_constants
from deepvariant import dv_utils
from deepvariant import dv_vcf_constants
from deepvariant import haplotypes
from deepvariant import logging_level
from deepvariant.protos import deepvariant_pb2
from deepvariant.python import merge_phased_reads as merge_phased_reads_lib
from deepvariant.python import postprocess_variants as postprocess_variants_lib
from deepvariant.small_model import inference as small_model_inference
from absl import app
import multiprocessing
from third_party.nucleus.io import sharded_file_utils
from third_party.nucleus.io import tabix
from third_party.nucleus.io import tfrecord
from third_party.nucleus.io import vcf
from third_party.nucleus.io.python import merge_variants
from third_party.nucleus.io.python import vcf_concat
from third_party.nucleus.protos import range_pb2
from third_party.nucleus.protos import reference_pb2
from third_party.nucleus.protos import variants_pb2
from third_party.nucleus.util import errors
from third_party.nucleus.util import genomics_math
from third_party.nucleus.util import proto_utils
from third_party.nucleus.util import ranges
from third_party.nucleus.util import struct_utils
from third_party.nucleus.util import variant_utils
from third_party.nucleus.util import variantcall_utils
_INFILE = flags.DEFINE_string(
'infile',
None,
(
'Required. Path(s) to CallVariantOutput protos in TFRecord format to '
'postprocess. These should be the complete set of outputs for '
'call_variants.py.'
),
)
_OUTFILE = flags.DEFINE_string(
'outfile',
None,
(
'Required. Destination path where we will write output variant calls in'
' VCF format.'
),
)
_REF = flags.DEFINE_string(
'ref',
None,
(
'Required. Genome reference in FAI-indexed FASTA format. Used to'
' determine the sort order for the emitted variants and the VCF header.'
),
)
_SMALL_MODEL_CVO_RECORDS = flags.DEFINE_string(
'small_model_cvo_records',
None,
(
'Optional. Path(s) to CallVariantOutput protos in TFRecord format that'
' were called by the small model to include in postprocess .'
),
)
_QUAL_FILTER = flags.DEFINE_float(
'qual_filter',
1.0,
'Any variant with QUAL < qual_filter will be filtered in the VCF file.',
)
_CNN_HOMREF_CALL_MIN_GQ = flags.DEFINE_float(
'cnn_homref_call_min_gq',
20.0,
(
'All CNN RefCalls whose GQ is less than this value will have ./.'
' genotype instead of 0/0.'
),
)
_MULT_ALLELIC_QUAL_FILTER = flags.DEFINE_float(
'multi_allelic_qual_filter',
1.0,
'The qual value below which to filter multi-allelic variants.',
)
_NONVARIANT_SITE_TFRECORD_PATH = flags.DEFINE_string(
'nonvariant_site_tfrecord_path',
None,
(
'Optional. Path(s) to the non-variant sites protos in TFRecord format'
' to convert to gVCF file. This should be the complete set of outputs'
' from the --gvcf flag of make_examples.py.'
),
)
_PHASED_READS_INPUT_PATH = flags.DEFINE_string(
'phased_reads_input_path',
None,
(
'Optional. Path to a TSV file containing phased read information, '
'typically an output of the make_examples step. This information '
'will be used to extend phase blocks in the output VCF.'
),
)
_CHECKPOINT_JSON = flags.DEFINE_string(
'checkpoint_json',
None,
'Optional. Path to the json file containing the flags for postprocessing.',
)
_PHASED_READS_SWITCHES_OUTPUT_PATH = flags.DEFINE_string(
'phased_reads_switches_output_path',
'/tmp/phased_reads_switches.tsv',
(
'Optional. Path to a TSV file containing switches information, '
'typically an output of the make_examples step. This information '
'will be used to extend phase blocks in the output VCF.'
),
)
_PHASED_READS_CORRECTED_OUTPUT_PATH = flags.DEFINE_string(
'phased_reads_corrected_output_path',
'/tmp/phased_reads_corrected.tsv',
(
'Optional. Path to a TSV file containing phased read information, '
'typically an output of the make_examples step. This information '
'will be used to extend phase blocks in the output VCF.'
),
)
_GVCF_OUTFILE = flags.DEFINE_string(
'gvcf_outfile',
None,
'Optional. Destination path where we will write the Genomic VCF output.',
)
_GROUP_VARIANTS = flags.DEFINE_boolean(
'group_variants',
True,
(
'If using vcf_candidate_importer and multi-allelic '
'sites are split across multiple lines in VCF, set to False so that '
'variants are not grouped when transforming CallVariantsOutput to '
'Variants.'
),
)
_VCF_STATS_REPORT = flags.DEFINE_boolean(
'vcf_stats_report',
False,
'Deprecated. Use vcf_stats_report.py instead.',
)
_SAMPLE_NAME = flags.DEFINE_string(
'sample_name',
None,
(
'Optional. If set, this will only be used if the sample name cannot be '
'determined from the CallVariantsOutput or non-variant sites protos.'
),
)
_USE_MULTIALLELIC_MODEL = flags.DEFINE_boolean(
'use_multiallelic_model',
False,
(
'If True, use a specialized model for genotype resolution of'
' multiallelic cases with two alts.'
),
)
_MULTIALLELIC_MODE = flags.DEFINE_enum(
'multiallelic_mode',
'product',
['min', 'product'],
'The fusion rule for merging probabilities in multiallelic calling.',
)
_DEBUG_OUTPUT_ALL_CANDIDATES = flags.DEFINE_enum(
'debug_output_all_candidates',
None,
['ALT', 'INFO'],
(
'Outputs all candidates considered by DeepVariant as additional ALT'
' alleles or as an INFO field. For ALT, filtered candidates are'
' assigned a GL=0 and added as ALTs alleles, but do not appear in any'
' sample genotypes. This flag is useful for debugging purposes.'
' ALT-mode is incompatible with the multiallelic caller.'
),
)
_ONLY_KEEP_PASS = flags.DEFINE_boolean(
'only_keep_pass', False, 'If True, only keep PASS calls.'
)
_HAPLOID_CONTIGS = flags.DEFINE_string(
'haploid_contigs',
None,
(
'Optional list of non autosomal chromosomes. For all listed chromosomes'
'HET probabilities are not considered. The list can be either comma '
'or space-separated.'
),
)
_CPUS = flags.DEFINE_integer(
'cpus',
multiprocessing.cpu_count(),
'Number of worker processes to use. Set --cpus < 2 to disable parallel'
' processing.',
short_name='j',
required=False,
)
_NUM_PARTITIONS = flags.DEFINE_integer(
'num_partitions',
0,
'Number of partitions to use for parallel or sequential processing. --Set'
' --num_partitions > --cpus to trade runtime for lower memory usage. Set'
' --num_partitions < 2 and --cpus < 2 to disable partitioning.',
required=False,
)
_PAR_REGIONS = flags.DEFINE_string(
'par_regions_bed',
None,
(
'Optional BED file containing Human Pseudoautosomal Region (PAR) '
'regions.'
'Variants within this region are unaffected by genotype reallocation '
'applied on regions supplied by --haploid_contigs flag.'
),
)
_REGIONS = flags.DEFINE_string(
'regions',
'',
(
'Optional. Space-separated list of regions we want to process. Elements'
' can be region literals (e.g., chr20:10-20) or paths to BED/BEDPE'
' files. This should match the flag passed to make_examples.py.'
),
)
_PROCESS_SOMATIC = flags.DEFINE_boolean(
'process_somatic',
False,
'Optional. If specified the input is treated as somatic.',
)
_PON_FILTERING = flags.DEFINE_string(
'pon_filtering',
None,
(
'Optional. Only used if --process_somatic is true. '
'A VCF file with Panel of Normals (PON) data.'
'If set, the output VCF will be filtered: any variants that appear in '
'PON will be marked with a PON filter, and PASS filter value will be '
'removed.'
),
)
_RESOLVE_CALL_VARIANTS_OUTPUTS_BY_MODEL = flags.DEFINE_bool(
'resolve_call_variants_outputs_by_model',
False,
'[Experimental] If true, postprocess_variants expects 2 CVO records for'
' each example, one from each model. One of the CVO record is chosen based'
' on the value of the `--small_model_gq_threshold` flag. This mirrors the'
' same behavior that normally happens in `make_examples` during inference.'
' The purpose is to explore the joint accuracy of the two'
' models as a function of GQ thresholds.',
)
_SMALL_MODEL_GQ_THRESHOLD = flags.DEFINE_integer(
'small_model_gq_threshold',
-1,
'[Experimental] The GQ threshold for accepting classifications from the'
' small model. Used when `--resolve_call_variants_outputs_by_model` is'
' true.',
)
# Some format fields are indexed by alt allele, such as AD (depth by allele).
# These need to be cleaned up if we remove any alt alleles. Any info field
# listed here will be have its values cleaned up if we've removed any alt
# alleles.
# Each tuple contains: field name, ref_is_zero.
_ALT_ALLELE_INDEXED_FORMAT_FIELDS = frozenset([
('AD', True),
('VAF', False),
('MF', True),
('MD', True),
('NAD', True),
('NAF', False),
])
# The number of places past the decimal point to round QUAL estimates to.
_QUAL_PRECISION = 7
# When this was set, it's about 20 seconds per log.
_LOG_EVERY_N = 100000
# When outputting all alt alleles, use placeholder value to indicate genotype
# will be soft-filtered.
_FILTERED_ALT_PROB = -9.0
# The number of genotype probabilities in a diploid sample.
_NUM_GENOTYPE_PROBABILITIES = 3
def _extract_single_sample_name(
record: deepvariant_pb2.CallVariantsOutput,
) -> str:
"""Returns the name of the single sample within the CallVariantsOutput file.
Args:
record: A deepvariant_pb2.CallVariantsOutput record.
Returns:
The name of the single individual in the first proto in the file.
Raises:
ValueError: There is not exactly one VariantCall in the proto or the
call_set_name of the VariantCall is not populated.
"""
variant = record.variant
call = variant_utils.only_call(variant)
name = call.call_set_name
if not name:
raise ValueError(
'Error extracting name: no call_set_name set: {}'.format(record)
)
return name
def _pysam_resolve_file_path(file_path: str) -> str:
"""Prepends a prefix to the file_path when accessing Google files.
Args:
file_path: str. Full path pointing a specific file to access with pysam.
Returns:
str. The full configured file path for pysam to open.
"""
# BEGN_INTERNAL
if (
file_path.startswith('/cns/')
or file_path.startswith('/placer/')
or file_path.startswith('/readahead/')
or file_path.startswith('/bigstore/')
):
return f'google:{file_path}'
# END_INTERNAL
return file_path
def most_likely_genotype(
predictions: Sequence[float], ploidy: int = 2, n_alleles: int = 2
) -> tuple[int, list[int]]:
"""Gets the most likely genotype from predictions.
From https://samtools.github.io/hts-specs/VCFv4.3.pdf:
Genotype Ordering. In general case of ploidy P and N alternate alleles (0 is
the REF and 1..N the alternate alleles), the ordering of genotypes for the
likelihoods can be expressed by the following pseudocode with as many nested
loops as ploidy:
* Note that we use inclusive for loop boundaries.
for a_P = 0 . . . N
for a_P-1 = 0 . . . aP
. . .
for a_1 = 0 . . . a2
println a1 a2 . . . aP
Alternatively, the same can be achieved recursively with the following
pseudocode:
Ordering (P , N , suffix =""):
for a in 0 . . . N
if (P == 1) println str (a) + suffix
if (P > 1) Ordering (P -1 , a, str (a) + suffix)
Examples:
* for P=2 and N=1, the ordering is 00,01,11
* for P=2 and N=2, the ordering is 00,01,11,02,12,22
* for P=3 and N=2, the ordering is 000,001,011,111,002,012,112,022,122,222
* for P=1, the index of the genotype a is a
* for P=2, the index of the genotype "a/b", where a <= b, is b(b + 1)/2 + a
* for P=2 and arbitrary N, the ordering can be easily derived from a
triangular matrix:
b / a 0 1 2 3
0 0
1 1 2
2 3 4 5
3 6 7 8 9
Args:
predictions: N element array-like. The real-space probabilities of each
genotype state for this variant. The number of elements in predictions is
related to ploidy and n_alleles is given by N = choose(ploidy + n_alleles
- 1, n_alleles -1) for more information see:
http://genome.sph.umich.edu/wiki/Relationship_between_Ploidy,_Alleles_and_Genotypes
ploidy: int >= 1. The ploidy (e.g., number of chromosomes) of this sample.
n_alleles: int >= 2. The number of alleles (ref + n_alts).
Returns:
Two values. The first is the index of the most likely prediction in
predictions. The second is a list of P elements with the VCF-style genotype
indices corresponding to this index. For example, with P = 2 and an index of
1, this returns the value (1, [0, 1]).
Raises:
NotImplementedError: if ploidy != 2 as this not yet implemented.
ValueError: If n_alleles < 2.
ValueError: If we cannot determine the genotype given prediction, n_alts,
and ploidy.
"""
# TODO: This can be memoized for efficiency.
if ploidy != 2:
raise NotImplementedError('Ploidy != 2 not yet implemented.')
if n_alleles < 2:
raise ValueError('n_alleles must be >= 2 but got', n_alleles)
# TODO: would be nice to add test that predictions has the right
# number of elements. But that would involve calculating the binomial
# coefficient of n_alleles and ploidy, which would be expensive. Probably
# need to memoize the whole function if we are going to add this.
index_of_max = np.argmax(predictions)
# This is the general case solution for fixed ploidy of 2 and arbitrary
# n_alleles. We should generalize this code to the arbitrary ploidy case when
# needed and memoize the mapping here.
index = 0
for h1 in range(0, n_alleles + 1):
for h2 in range(0, h1 + 1):
if index == index_of_max:
return index, [h2, h1]
index += 1
raise ValueError('No corresponding GenotypeType for predictions', predictions)
def uncall_gt_if_no_ad(variant: variants_pb2.Variant) -> None:
"""Converts genotype to "./." if sum(AD)=0."""
vcall = variant_utils.only_call(variant)
if sum(variantcall_utils.get_ad(vcall)) == 0:
# Set GT to ./.; GLs set to 0; GQ=0
vcall.genotype[:] = [-1, -1]
vcall.genotype_likelihood[:] = [0, 0]
variantcall_utils.set_gq(vcall, 0)
def uncall_homref_gt_if_lowqual(
variant: variants_pb2.Variant, min_homref_gq: float
) -> None:
"""Converts genotype to "./." if variant is CNN RefCall and has low GQ.
If the variant has "RefCall" filter (which means an example was created for
this site but CNN didn't call this as variant) and if the GQ is less than
the given min_homref_gq threshold, set the genotype of the variant proto
to "./.". See http://internal for more info.
Args:
variant: third_party.nucleus.protos.Variant proto.
min_homref_gq: float.
"""
vcall = variant_utils.only_call(variant)
if (
variant.filter == [dv_vcf_constants.DEEP_VARIANT_REF_FILTER]
and variantcall_utils.get_gq(vcall) < min_homref_gq
):
vcall.genotype[:] = [-1, -1]
variant.filter[:] = [dv_vcf_constants.DEEP_VARIANT_NO_CALL]
# TODO Implement ingration test to test phased output.
def maybe_phase_genotype(
variant: variants_pb2.Variant,
genotype: list[int],
) -> tuple[bool, list[int]]:
"""Phases the genotype if phase information is available.
The `ALT_PS` field contains phases assigned to each allele in range [0..2],
for HP tags 0,1,2. The length of this array is number of alleles + 1, since
REF is added implicitly.
For example, [1,2] means that REF allele is assigned phase 1 and ALT_1 allele
is assigned a phase 2. [2,1] means that REF allele is assigned phase 2 and
ALT_1 allele is assigned phase 1. [2,2,1,1] means REF is phase 2, ALT_1 is
phase 2, ALT_2 is phase 1, and ALT_3 is phase 1, etc.
Args:
variant: third_party.nucleus.protos.Variant proto.
genotype: list of ints. The genotype indices to be written to the VCF.
Returns:
is_phased: bool. Whether it was possible to phase the genotype.
genotype: genotype in the correct phasing order, if phased.
"""
if not (
variant_utils.get_info(variant, dv_constants.VARIANT_PHASE_SET)
and variant_utils.get_info(variant, dv_constants.PHASED_GENOTYPE)
):
return False, genotype
phase_info = [
p.int_value for p in variant.info[dv_constants.PHASED_GENOTYPE].values
]
if max(genotype) >= len(phase_info):
logging.warning(
(
'Genotype %s is out of range for phase info %s for variant %s. '
'Phasing was not applied.'
),
genotype,
phase_info,
variant,
)
return False, genotype
allele_1_haplotype = phase_info[genotype[0]]
allele_2_haplotype = phase_info[genotype[1]]
is_phased = (
0 not in (allele_1_haplotype, allele_2_haplotype)
and allele_1_haplotype != allele_2_haplotype
)
if is_phased:
genotype = [
genotype[allele_1_haplotype - 1],
genotype[allele_2_haplotype - 1],
]
return is_phased, genotype
def add_call_to_variant(
variant: variants_pb2.Variant,
predictions: Sequence[float],
qual_filter: float,
sample_name: str | None,
) -> variants_pb2.Variant:
"""Fills in Variant record using the prediction probabilities.
This functions sets the call[0].genotype, call[0].info['GQ'],
call[0].genotype_probabilities, variant.filter, and variant.quality fields of
variant based on the genotype likelihoods in predictions.
Args:
variant: third_party.nucleus.protos.Variant protobuf to be filled in with
info derived from predictions.
predictions: N element array-like. The real-space probabilities of each
genotype state for this variant.
qual_filter: float. If predictions implies that this isn't a reference call
and the QUAL of the prediction isn't larger than qual_filter variant will
be marked as FILTERed.
sample_name: str. The name of the sample to assign to the Variant proto
call_set_name field.
Returns:
A tuple of the Variant record and its phase information.
Raises:
ValueError: If variant doesn't have exactly one variant.call record.
"""
call = variant_utils.only_call(variant)
n_alleles = len(variant.alternate_bases) + 1
index, genotype = most_likely_genotype(predictions, n_alleles=n_alleles)
gq, variant.quality = compute_quals(predictions, index)
call.call_set_name = sample_name
call.is_phased, genotype = maybe_phase_genotype(
variant,
genotype,
)
if is_methylated(call):
mf = variantcall_utils.get_mf(call)
mt = variantcall_utils.determine_methylation_type(mf)
mi = variantcall_utils.get_mi(call)
variantcall_utils.set_mt(call, mt)
variantcall_utils.set_mi(call, mi)
variantcall_utils.set_gt(call, genotype)
variantcall_utils.set_gq(call, gq)
gls = [genomics_math.perror_to_bounded_log10_perror(gp) for gp in predictions]
variantcall_utils.set_gl(call, gls)
uncall_gt_if_no_ad(variant)
variant.filter[:] = dv_vcf_constants.compute_filter_fields(
variant, qual_filter
)
uncall_homref_gt_if_lowqual(variant, _CNN_HOMREF_CALL_MIN_GQ.value)
return variant
def compute_quals(
predictions: Sequence[float], prediction_index: int
) -> tuple[int, int]:
"""Computes GQ and QUAL values from a set of prediction probabilities.
Prediction probabilities are represented as a probability distribution over
the N genotype states (e.g., for 3 genotype states {HOM_REF, HET, HOM_VAR}).
Genotype Quality (or GQ) represents the PHRED scaled confidence in the
particular genotype assignment. Likewise the QUAL representes the PHRED scaled
confidence in variant as compared to reference, that is, P(NON_REF) / P(ALL)
which in the diploid genotype case is P(HET) + P(HOM_VAR) / P(ALL). These
quality scores are capped by _MAX_CONFIDENCE.
Args:
predictions: N element array-like. The real-space probabilities of each
genotype state for this variant.
prediction_index: int. The actual called genotype from the distribution.
Returns:
GQ and QUAL values for output in a Variant record.
"""
# GQ is prob(genotype) / prob(all genotypes)
# GQ is rounded to the nearest integer to comply with the VCF spec.
gq = int(
np.around(
genomics_math.ptrue_to_bounded_phred(predictions[prediction_index])
)
)
# QUAL is prob(variant genotype) / prob(all genotypes)
# Taking the min to avoid minor numerical issues than can push sum > 1.0.
# TODO: this is equivalent to the likely better implementation:
# genomics_math.perror_to_phred(max(predictions[0], min_ref_confidence))
# where min_ref_confidence is roughly 1.25e-10 (producing a qual of 99).
qual = genomics_math.ptrue_to_bounded_phred(min(sum(predictions[1:]), 1.0))
rounded_qual = round(qual, _QUAL_PRECISION)
return gq, rounded_qual
def expected_alt_allele_indices(num_alternate_bases: int) -> list[list[int]]:
"""Returns (sorted) expected list of alt_allele_indices, given #alt bases."""
num_alleles = num_alternate_bases + 1
alt_allele_indices_list = [
sorted(list(set(x) - {0}))
for x in itertools.combinations(range(num_alleles), 2)
]
# alt_allele_indices starts from 0, where 0 refers to the first alt allele.
# pylint: disable=g-complex-comprehension
return sorted([
[i - 1 for i in alt_allele_indices]
for alt_allele_indices in alt_allele_indices_list
])
# pylint: enable=g-complex-comprehension
def _check_alt_allele_indices(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
) -> bool:
"""Returns True if and only if the alt allele indices are valid."""
all_alt_allele_indices = sorted([
list(call_variants_output.alt_allele_indices.indices)
for call_variants_output in call_variants_outputs
])
if all_alt_allele_indices != expected_alt_allele_indices(
len(call_variants_outputs[0].variant.alternate_bases)
):
logging.warning(
(
'Alt allele indices found from call_variants_outputs for '
'variant %s is %s, which is invalid.'
),
call_variants_outputs[0].variant,
all_alt_allele_indices,
)
return False
return True
def variants_are_equal(
variant_1: variants_pb2.Variant,
variant_2: variants_pb2.Variant,
) -> bool:
"""Returns True if the variants are the same, ignoring the calls field.
The `calls` field is ignored because the `VariantCall` object might originate
from the small model (during make_examples) or the CNN model (during
call_variants); what matters is that the `Variant` objects themselves match.
Args:
variant_1: The first variant to compare.
variant_2: The second variant to compare.
Returns:
True if the variants are the same, otherwise False.
"""
variant_1_vars = json_format.MessageToDict(variant_1)
variant_2_vars = json_format.MessageToDict(variant_2)
del variant_1_vars['calls']
del variant_2_vars['calls']
if 'info' in variant_1_vars:
del variant_1_vars['info']
if 'info' in variant_2_vars:
del variant_2_vars['info']
return variant_1_vars == variant_2_vars
def is_valid_call_variants_outputs(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
) -> bool:
"""Returns True if the call_variants_outputs follows our assumptions.
Args:
call_variants_outputs: list of CallVariantsOutput to check.
Returns:
True if the sanity check passes.
"""
if not call_variants_outputs:
return True # An empty list is a degenerate case.
if (
not _check_alt_allele_indices(call_variants_outputs)
and not _RESOLVE_CALL_VARIANTS_OUTPUTS_BY_MODEL.value
):
return False
first_call, other_calls = call_variants_outputs[0], call_variants_outputs[1:]
# Sanity check that all call_variants_outputs have the same `variant`.
for call_to_check in other_calls:
if not variants_are_equal(first_call.variant, call_to_check.variant):
logging.warning(
(
'Expected all inputs to merge_predictions to have the '
'same `variant`, but getting %s and %s.'
),
first_call.variant,
call_to_check.variant,
)
return False
return True
def convert_call_variants_outputs_to_probs_dict(
canonical_variant: variants_pb2.Variant,
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
alt_alleles_to_remove: set[str],
debug_output_all_candidates: str | None = None,
) -> dict[tuple[str, str], list[float]]:
"""Converts a list of CallVariantsOutput to an internal allele probs dict.
Args:
canonical_variant: variants_pb2.Variant.
call_variants_outputs: list of CallVariantsOutput.
alt_alleles_to_remove: set of strings. Alleles to remove.
debug_output_all_candidates: If 'ALT', set low qual alleles to be
soft-filtered.
Returns:
Dictionary of {(allele1, allele2): list of probabilities},
where allele1 and allele2 are strings.
"""
flattened_dict = collections.defaultdict(list)
if not call_variants_outputs:
return flattened_dict
for call_variants_output in call_variants_outputs:
allele_set1 = frozenset([canonical_variant.reference_bases])
allele_set2 = frozenset(
canonical_variant.alternate_bases[index]
for index in call_variants_output.alt_allele_indices.indices
)
has_alleles_to_rm = bool(alt_alleles_to_remove.intersection(allele_set2))
if has_alleles_to_rm and debug_output_all_candidates != 'ALT':
continue
if has_alleles_to_rm:
# This block is run when debug_output_all_candidates=ALT
# It sets genotype likelihood to a placeholder value,
# which is later used to set GL=1.0 (prob=0).
p11, p12, p22 = (
_FILTERED_ALT_PROB,
_FILTERED_ALT_PROB,
_FILTERED_ALT_PROB,
)
else:
p11, p12, p22 = call_variants_output.genotype_probabilities
for set1, set2, p in [
(allele_set1, allele_set1, p11),
(allele_set1, allele_set2, p12),
(allele_set2, allele_set2, p22),
]:
for indices in itertools.product(set1, set2):
flattened_dict[indices].append(p)
return flattened_dict
def get_alt_alleles_to_remove(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
qual_filter: float,
) -> set[str]:
"""Returns all the alt alleles with quality below qual_filter.
Quality is defined as (1-p(ref/ref)). This removes all alt alleles whose
quality is below the filter value, with the exception that if the set of
removed alt alleles covers everything in the alternate_bases, the single alt
allele where the 1-p(ref/ref) is the highest is retained.
Args:
call_variants_outputs: list of CallVariantsOutput.
qual_filter: double. The qual value below which to filter variants.
Returns:
Set of strings: alt alleles to remove.
"""
alt_alleles_to_remove = set() # first alt is represented as 0.
if not qual_filter or not call_variants_outputs:
return alt_alleles_to_remove
max_qual, max_qual_allele = None, None
canonical_variant = call_variants_outputs[0].variant
for call_variants_output in call_variants_outputs:
# Go through the ones where alt_allele_indices has
# exactly one element. There are the pileup images that contains information
# like:
# p00, p01, p11
# or p00, p02, p22
# ...p00, p0N, pNN
if len(call_variants_output.alt_allele_indices.indices) == 1:
# From here, we want to see which ones of these alt alleles (1-N) that we
# can skip. We can use the concept of QUAL in VCF, and filter out ones
# where QUAL < FLAGS.qual_filter. This is because if QUAL is too low,
# it means it is unlikely this has a variant genotype.
_, qual = compute_quals(
call_variants_output.genotype_probabilities, prediction_index=0
)
alt_allele_index = call_variants_output.alt_allele_indices.indices[0]
# Keep track of one alt allele with the highest qual score.
if max_qual is None or max_qual < qual:
max_qual, max_qual_allele = (
qual,
canonical_variant.alternate_bases[alt_allele_index],
)
if qual < qual_filter:
alt_alleles_to_remove.add(
canonical_variant.alternate_bases[alt_allele_index]
)
# If all alt alleles are below `qual_filter`, keep at least one.
if len(alt_alleles_to_remove) == len(canonical_variant.alternate_bases):
alt_alleles_to_remove -= set([max_qual_allele])
return alt_alleles_to_remove
def is_methylated(call: variants_pb2.VariantCall) -> bool:
"""Determines if a VariantCall is methylated.
A variant is considered methylated if any of its methylation fractions
(MF) for the reference or alternate alleles is greater than 0.
Args:
call: A `variants_pb2.VariantCall` object.
Returns:
bool: True if any `methylation_fraction` value is greater than 0,
otherwise False.
"""
if 'MF' not in call.info:
return False
mf_values = variantcall_utils.get_mf(call)
if any(mf > 0 for mf in mf_values):
return True
return False
class AlleleRemapper:
"""Facilitates removing alt alleles from a Variant.
This class provides a one-to-shop for managing the information needed to
remove alternative alleles from Variant. It provides functions and properties
to get the original alts, the new alts, and asking if alleles (strings) or
indices (integers) should be retained or eliminated.
"""
def __init__(
self, original_alt_alleles: Sequence[str], alleles_to_remove: set[str]
):
self.original_alts = list(original_alt_alleles)
self.alleles_to_remove = set(alleles_to_remove)
def keep_index(
self, allele_index: int, ref_is_zero: bool | str = False
) -> bool:
if ref_is_zero:
return True if allele_index == 0 else self.keep_index(allele_index - 1)
else:
return self.original_alts[allele_index] not in self.alleles_to_remove
def retained_alt_alleles(self) -> Sequence[str]:
return [
alt for alt in self.original_alts if alt not in self.alleles_to_remove
]
def reindex_allele_indexed_fields(
self, variant: variants_pb2.Variant, fields: frozenset[tuple[str, bool]]
) -> None:
"""Updates variant.call fields indexed by ref + alt_alleles.
Args:
variant: Variant proto. We will update the info fields of the Variant.call
protos.
fields: Iterable of string. Each string should provide a key to an
alternative allele indexed field in VariantCall.info fields. Each field
specified here will be updated to remove values associated with alleles
no longer wanted according to this remapper object.
"""
for field_info in fields:
field = field_info[0]
ref_is_zero = field_info[1]
for call in variant.calls:
if field in call.info:
entry = call.info[field]
updated = [
v
for i, v in enumerate(entry.values)
if self.keep_index(i, ref_is_zero=ref_is_zero)
]
# We cannot do entry.values[:] = updated as the ListValue type "does
# not support assignment" so we have to do this grossness.
del entry.values[:]
entry.values.extend(updated)
def prune_alleles(
variant: variants_pb2.Variant, alt_alleles_to_remove: set[str]
) -> variants_pb2.Variant:
"""Remove the alt alleles in alt_alleles_to_remove from canonical_variant.
Args:
variant: variants_pb2.Variant.
alt_alleles_to_remove: iterable of str. Alt alleles to remove from variant.
Returns:
variants_pb2.Variant with the alt alleles removed from alternate_bases.
"""
# If we aren't removing any alt alleles, just return the unmodified variant.
if not alt_alleles_to_remove:
return variant
new_variant = variants_pb2.Variant()
new_variant.CopyFrom(variant)
# Cleanup any VariantCall.info fields indexed by alt allele.
remapper = AlleleRemapper(variant.alternate_bases, alt_alleles_to_remove)
remapper.reindex_allele_indexed_fields(
new_variant, _ALT_ALLELE_INDEXED_FORMAT_FIELDS
)
new_variant.alternate_bases[:] = remapper.retained_alt_alleles()
return new_variant
def get_multiallelic_distributions(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
pruned_alleles: set[str],
) -> np.ndarray:
"""Return 9 values for 3 distributions from given multiallelic CVOs.
This function is only called for sites with two alt alleles remaining after
pruning. However, call_variants_outputs contains CVOs from pruned and unpruned
alleles, so we ignore the CVOs containing alleles that were pruned.
Args:
call_variants_outputs: list of CVOs for a multiallelic site with exactly two
alts after pruning. For such a site, we would expect 3 CVOs (alt1, alt2,
alt1/2). However, there may be more than 3 CVOs if some alleles were
pruned at this site.
pruned_alleles: set of strings corresponding to pruned alleles. Used to
filter CVOs for pruned alleles.
Returns:
final_probs: array of shape (1, 9). The 9 values correspond to three model
output distributions. The first is from the image containing alt1, the
second is from the image for alt2, the third is from the image with both
alt1 and alt2.
"""
alt_allele_indices_to_probs = {}
first_alt_index = None
second_alt_index = None
# Find the CVOs with two alts, corresponding to the image with alt1 and alt2.
for cvo in call_variants_outputs:
indices = cvo.alt_allele_indices.indices[:]
curr_alleles = [cvo.variant.alternate_bases[i] for i in indices]
curr_alleles_pruned = any([a in pruned_alleles for a in curr_alleles])
# Ignore CVOs containing pruned alleles.
if len(indices) == 2 and not curr_alleles_pruned:
first_alt_index = min(indices)
second_alt_index = max(indices)
probs = cvo.genotype_probabilities[:]
alt_allele_indices_to_probs[(first_alt_index, second_alt_index)] = probs
# Find the single alt CVOs.
for cvo in call_variants_outputs:
if len(cvo.alt_allele_indices.indices[:]) == 1:
index = cvo.alt_allele_indices.indices[0]
if index == first_alt_index or index == second_alt_index:
probs = cvo.genotype_probabilities[:]
alt_allele_indices_to_probs[index] = probs
assert len(alt_allele_indices_to_probs) == 3
# Concatenate all probabilities into one array.
final_probs = np.array([
alt_allele_indices_to_probs[first_alt_index]
+ alt_allele_indices_to_probs[second_alt_index]
+ alt_allele_indices_to_probs[(first_alt_index, second_alt_index)]
])
return final_probs
@functools.lru_cache
def get_multiallelic_model(
use_multiallelic_model: bool,
) -> tf.keras.Model | None:
"""Loads and returns the model, which must be in saved model format.
Args:
use_multiallelic_model: if True, use a specialized model for genotype
resolution of multiallelic cases with two alts.
Returns:
A keras model instance if use_multiallelic_model, else None.
"""
if not use_multiallelic_model:
return None
curr_dir = os.path.dirname(__file__)
multiallelic_model_path = os.path.join(curr_dir, 'multiallelic_model')
return tf.keras.models.load_model(multiallelic_model_path, compile=False)
def normalize_predictions(predictions: Sequence[float]) -> Sequence[float]:
"""Normalize predictions and handle soft-filtered alt alleles."""
if sum(predictions) == 0:
predictions = [1.0] * len(predictions)
denominator = (
sum([i if i != _FILTERED_ALT_PROB else 0.0 for i in predictions]) or 1.0
)
normalized_predictions = [
i / denominator if i != _FILTERED_ALT_PROB else 0.0 for i in predictions
]
return normalized_predictions
def correct_nonautosome_probabilities(
probabilities: list[float],
variant: variants_pb2.Variant,
) -> Sequence[float]:
"""Recalculate probabilities for non-autosome heterozygous calls."""
n_alleles = len(variant.alternate_bases) + 1
# It is assumed that probabilities are stored in the specific order. See
# most_likely_genotype for details.
# Each heterozyhous probability is zeroed. For example, for biallelic case
# the probability of 0/1 genotype becomes zero.
index = 0
for h1 in range(0, n_alleles):
for h2 in range(0, h1 + 1):
if h2 != h1:
if len(probabilities) <= index:
raise ValueError("Probabilties array doesn't match alt alleles.")
probabilities[index] = 0
index += 1
new_sum = sum(probabilities) or 1.0
return list(map(lambda p: p / new_sum, probabilities))
def is_non_autosome(variant: variants_pb2.Variant) -> bool:
"""Returns True if variant is non_autosome."""
haploid_contigs_str = _HAPLOID_CONTIGS.value or ''
parts = haploid_contigs_str.split(',')
# pylint: disable=g-complex-comprehension
haploid_contigs = [item for part in parts for item in part.split()]
# pylint: enable=g-complex-comprehension
return haploid_contigs and variant.reference_name in haploid_contigs # pytype: disable=bad-return-type
def is_in_regions(
variant: variants_pb2.Variant, regions: ranges.RangeSet
) -> bool:
"""Returns True of variant overlaps one of the regions."""
if regions:
return regions.variant_overlaps(variant)
else:
return False
@functools.lru_cache
def get_par_regions() -> ranges.RangeSet | None:
"""Returns the cached par regions if specified, else None."""
if _PAR_REGIONS.value:
return ranges.RangeSet.from_bed(_PAR_REGIONS.value, enable_logging=False)
else:
return None
def resolve_call_variant_outputs_by_model(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
gq_threshold: float,
) -> Sequence[deepvariant_pb2.CallVariantsOutput]:
"""Picks one of the CVO pairs by model ID based on the GQ threshold.
In GQ debug mode, CVOs are written by both models. This means for every
candidate <> alt_allele_indices combination, there are 2 CVOs with class
probabilities from each model. For every alt_allele_indices set, one of the
CVOs is picked based on the GQ threshold, following the same logic as happens
in make_examples: if the GQ threshold is met, the CVO from the small model is
picked. Otherwise, the CVO from the deepvariant model is picked. This allows
multiple postprocess_variant commands to produce many VCFs from a single
deepvariant run, allowing for a detailed analysis of which GQ threshold
is optimal.
Args:
call_variants_outputs: list of CVOs for a given variant site.
gq_threshold: GQ threshold.
Returns:
A list of CVOs, one for each alt_allele_indices set from one of the models.
"""
filtered_by_gq = []
cvos_grouped_by_alt_allele_indices = itertools.groupby(
call_variants_outputs,
lambda x: x.alt_allele_indices,
)
for _, cvo_pair in cvos_grouped_by_alt_allele_indices:
# every group should have 2 CVOs, one for each model, which are sorted into
# alphabetical order by model ID: ["deepvariant", "small_model"]
deepvariant_call, small_model_call = sorted(
cvo_pair,
key=lambda x: variantcall_utils.get_model_id(x.variant.calls[0]),
)
if small_model_inference.passes_confidence_threshold(
small_model_call.genotype_probabilities, gq_threshold
):
filtered_by_gq.append(small_model_call)
else:
filtered_by_gq.append(deepvariant_call)
return filtered_by_gq
def merge_predictions(
call_variants_outputs: Sequence[deepvariant_pb2.CallVariantsOutput],
qual_filter: float | None = None,
multiallelic_model: tf.keras.Model | None = None,
debug_output_all_candidates: str | None = None,
) -> tuple[variants_pb2.Variant, Sequence[float]]:
"""Merges the predictions from the multi-allelic calls."""
# See the logic described in the class PileupImageCreator pileup_image.py
#
# Because of the logic above, this function expects all cases above to have
# genotype_predictions that we can combine from.
if _RESOLVE_CALL_VARIANTS_OUTPUTS_BY_MODEL.value:
call_variants_outputs = resolve_call_variant_outputs_by_model(
call_variants_outputs, _SMALL_MODEL_GQ_THRESHOLD.value
)
par_regions = get_par_regions()
if not call_variants_outputs:
raise ValueError('Expected 1 or more call_variants_outputs.')
if not is_valid_call_variants_outputs(call_variants_outputs):
raise ValueError('`call_variants_outputs` did not pass sanity check.')
first_call, other_calls = call_variants_outputs[0], call_variants_outputs[1:]
canonical_variant = first_call.variant
if not other_calls:
canonical_variant = variant_utils.simplify_variant_alleles(
canonical_variant
)
if is_non_autosome(canonical_variant) and not is_in_regions(
canonical_variant, par_regions
):
return canonical_variant, correct_nonautosome_probabilities(
list(first_call.genotype_probabilities), canonical_variant
)
return canonical_variant, first_call.genotype_probabilities
# Special handling of multiallelic variants
alt_alleles_to_remove = get_alt_alleles_to_remove(
call_variants_outputs, qual_filter
)
# flattened_probs_dict is only used with the multiallelic model
flattened_probs_dict = convert_call_variants_outputs_to_probs_dict(
canonical_variant,
call_variants_outputs,
alt_alleles_to_remove,
debug_output_all_candidates,
)
if debug_output_all_candidates == 'INFO':
struct_utils.add_string_field(
canonical_variant.info,
'CANDIDATES',
'|'.join(canonical_variant.alternate_bases),
)
if debug_output_all_candidates != 'ALT':
canonical_variant = prune_alleles(canonical_variant, alt_alleles_to_remove)
# Run alternate model for multiallelic cases.
num_alts = len(canonical_variant.alternate_bases)
if num_alts == 2 and multiallelic_model is not None:
# We have 3 CVOs for 2 alts. In this case, there are 6 possible genotypes.
cvo_probs = get_multiallelic_distributions(
call_variants_outputs, alt_alleles_to_remove
)
normalized_predictions = multiallelic_model(cvo_probs).numpy().tolist()[0]
elif _MULTIALLELIC_MODE.value == 'product':
# New logic: "overlap-count" with product fusion.
# 1. Collect information about each CVO's example.
example_info = []
original_variant = call_variants_outputs[0].variant
for cvo in call_variants_outputs:
example_alt_alleles = frozenset(
original_variant.alternate_bases[i]
for i in cvo.alt_allele_indices.indices
)
is_for_pruned_allele = bool(
alt_alleles_to_remove.intersection(example_alt_alleles)
)
if is_for_pruned_allele and debug_output_all_candidates != 'ALT':
continue
probs = (
(_FILTERED_ALT_PROB,) * _NUM_GENOTYPE_PROBABILITIES
if is_for_pruned_allele
else cvo.genotype_probabilities
)
example_info.append({'probs': probs, 'alts': example_alt_alleles})
# 2. Calculate raw probability for each possible genotype.
predictions = []
genotype_ordering = variant_utils.genotype_ordering_in_likelihoods(
canonical_variant
)
for _, _, allele1_str, allele2_str in genotype_ordering:
prob_list_for_genotype = []
for example in example_info:
# Check each allele of the diploid genotype independently against the
# example's alternate alleles. This correctly calculates an overlap of 2
# for homozygous alternate genotypes.
overlap = int(allele1_str in example['alts']) + int(
allele2_str in example['alts']
)
prob_list_for_genotype.append(example['probs'][overlap])
# 3. Fuse probabilities with product.
if _FILTERED_ALT_PROB in prob_list_for_genotype:
fused_prob = _FILTERED_ALT_PROB
else:
fused_prob = np.prod(prob_list_for_genotype)
predictions.append(fused_prob)
# 4. Normalize the final predictions.
normalized_predictions = normalize_predictions(predictions)
else:
def min_alt_filter(probs):
return min([x for x in probs if x != _FILTERED_ALT_PROB] or [0])
predictions = [
min_alt_filter(flattened_probs_dict[(m, n)])
for _, _, m, n in variant_utils.genotype_ordering_in_likelihoods(
canonical_variant
)
]
if sum(predictions) == 0:
predictions = [1.0] * len(predictions)
normalized_predictions = normalize_predictions(predictions)
# Note the simplify_variant_alleles call *must* happen after the predictions
# calculation above. flattened_probs_dict is indexed by alt allele, and
# simplify can change those alleles so we cannot simplify until afterwards.
canonical_variant = variant_utils.simplify_variant_alleles(canonical_variant)
if is_non_autosome(canonical_variant) and not is_in_regions(
canonical_variant, par_regions
):
return canonical_variant, correct_nonautosome_probabilities(
normalized_predictions, canonical_variant
)
else:
return canonical_variant, normalized_predictions
def should_filter(
variant: variants_pb2.Variant,
pon_vcf_reader: vcf.VcfReader,
padding_bases: int = 0,
) -> bool:
"""Returns True if the variant should be filtered based on PON."""
if pon_vcf_reader is None:
return False
query_region = ranges.make_range(
chrom=variant.reference_name,
start=variant.start - padding_bases,
end=variant.end + padding_bases,
)
pon_variants = list(pon_vcf_reader.query(query_region))
if not pon_variants:
return False
# TODO: Consider improving this logic to directly match the
# contig name, the position, and REF and ALT directly.
variant_key = variant_utils.variant_key(variant)
for pon_variant in pon_variants:
if variant_key == variant_utils.variant_key(pon_variant):
return True
return False
def add_pon_filter(
variant_generator: Iterator[variants_pb2.Variant],
pon_vcf_reader: vcf.VcfReader,
) -> Iterator[variants_pb2.Variant]:
for variant in variant_generator:
if dv_vcf_constants.DEEP_VARIANT_PASS in variant.filter and should_filter(
variant, pon_vcf_reader
):
variant.filter.remove(dv_vcf_constants.DEEP_VARIANT_PASS)
variant.filter.append(dv_vcf_constants.DEEP_VARIANT_PON)
yield variant
def write_variants_to_vcf(
variant_iterable: Iterator[variants_pb2.Variant],
output_vcf_path: str,
header: variants_pb2.VcfHeader,
):
"""Writes Variant protos to a VCF file.
Args:
variant_iterable: iterable. An iterable of sorted Variant protos.
output_vcf_path: str. Output file in VCF format.
header: VcfHeader proto. The VCF header to use for writing the variants.
"""
logging.info('Writing output to VCF file: %s', output_vcf_path)
with vcf.VcfWriter(
output_vcf_path, header=header, round_qualities=True
) as writer:
count = 0
for variant in variant_iterable:
if not _ONLY_KEEP_PASS.value or variant.filter == [
dv_vcf_constants.DEEP_VARIANT_PASS
]:
count += 1
if _PROCESS_SOMATIC.value:
writer.write_somatic(variant)
else:
writer.write(variant)
logging.log_every_n(
logging.INFO, '%s variants written.', _LOG_EVERY_N, count
)
logging.info('Total variants written: %s', count)
def _sort_grouped_variants(group: Sequence[deepvariant_pb2.CallVariantsOutput]):
return sorted(group, key=lambda x: sorted(x.alt_allele_indices.indices))
def _transform_call_variant_group_to_output_variant(
call_variant_group: Sequence[deepvariant_pb2.CallVariantsOutput],
qual_filter: float,
multi_allelic_qual_filter: float,
sample_name: str,
use_multiallelic_model: bool,
debug_output_all_candidates: str | None,
) -> variants_pb2.Variant:
"""Transforms a group of CalVariantOutput to VariantOutput.
The group of CVOs present in the call_variants_group are converted to the
Variant proto, with the following filters applied: 1) variants are omitted
if their quality is lower than the `qual_filter` threshold. 2) multi-allelic
variants omit individual alleles whose qualities are lower than the
`multi_allelic_qual_filter` threshold.
Args:
call_variant_group: list[CVO]. A group of CallVariantsOutput protos.
qual_filter: double. The qual value below which to filter variants.
multi_allelic_qual_filter: double. The qual value below which to filter
multi-allelic variants.
sample_name: str. Sample name to write to VCF file.
use_multiallelic_model: if True, use a specialized model for genotype
resolution of multiallelic cases with two alts.
debug_output_all_candidates: if 'ALT', output all alleles considered by
DeepVariant as ALT alleles.
Returns:
the Variant proto.
"""
multiallelic_model = get_multiallelic_model(
use_multiallelic_model=use_multiallelic_model
)
outputs = _sort_grouped_variants(call_variant_group)
canonical_variant, predictions = merge_predictions(
outputs,
multi_allelic_qual_filter,
multiallelic_model=multiallelic_model,
debug_output_all_candidates=debug_output_all_candidates,
)
return add_call_to_variant(
canonical_variant,
predictions,
qual_filter=qual_filter,
sample_name=sample_name,
)
def _transform_call_variants_output_to_variants(
input_sorted_tfrecord_path: str,
sample_name: str,
) -> Iterator[variants_pb2.Variant]:
"""Yields Variant protos in sorted order from CallVariantsOutput protos.
Args:
input_sorted_tfrecord_path: str. TFRecord format file containing sorted
CallVariantsOutput protos.
sample_name: str. Sample name use in the output VCF and gVCF.
Yields:
Variant protos in sorted order representing the CallVariantsOutput calls.
"""
for call_variant_group in group_call_variants_outputs(
input_sorted_tfrecord_path, _GROUP_VARIANTS.value
):
yield _transform_call_variant_group_to_output_variant(
call_variant_group,
_QUAL_FILTER.value,
_MULT_ALLELIC_QUAL_FILTER.value,
sample_name,
_USE_MULTIALLELIC_MODEL.value,
_DEBUG_OUTPUT_ALL_CANDIDATES.value,
)
def dump_variants_to_temp_file(
variant_protos: Iterator[variants_pb2.Variant],
) -> tempfile._TemporaryFileWrapper:
temp = tempfile.NamedTemporaryFile()
tfrecord.write_tfrecords(variant_protos, temp.name)
return temp
def group_call_variants_outputs(
input_sorted_tfrecord_path: str, group_variants: bool
) -> Iterator[Sequence[deepvariant_pb2.CallVariantsOutput]]:
"""Yields CallVariantOutputs grouped by their variant range.
Args:
input_sorted_tfrecord_path: str. TFRecord format file containing sorted
CallVariantsOutput protos.
group_variants: bool. If true, group variants that have same start and end
position.
"""
group_fn = None
if group_variants:
group_fn = lambda x: variant_utils.variant_range(x.variant)
for _, group in itertools.groupby(
tfrecord.read_tfrecords(
input_sorted_tfrecord_path, proto=deepvariant_pb2.CallVariantsOutput
),
group_fn,
):
yield list(group)
def _concat_vcf(
output_file: str, temp_vcf_files: Sequence[tempfile._TemporaryFileWrapper]
) -> None:
"""Concatenates a set of temp (g)VCF files."""
vcf_files_to_concat = [f.name for f in temp_vcf_files]
vcf_concat.concat(output_file, vcf_files_to_concat)
def process_contiguous_partition(
contiguous_range_set: Sequence[range_pb2.Range],
contigs: Sequence[reference_pb2.ContigInfo],
cvo_paths: Sequence[str],
temp_file_name: str,
sample_name: str,
) -> Iterator[variants_pb2.Variant]:
"""Postprocess all CVOs in the given partition and returns an iterator.
Args:
contiguous_range_set: set of contiguous ranges to load and transform.
contigs: all contigs from ref
cvo_paths: paths to all CVO files
temp_file_name: path to temp file to variants to.
sample_name: the sample name to use for the output VCF and gVCF.
Returns:
An iterator of processed variants.
"""
start_time = time.time()
num_cvo_records = postprocess_variants_lib.process_single_sites_tfrecords(
contigs,
cvo_paths,
temp_file_name,
contiguous_range_set,
)
if contiguous_range_set:
logging.info(
'Processing region %s:%s-%s:%s',
contiguous_range_set[0].reference_name,
contiguous_range_set[0].start,
contiguous_range_set[-1].reference_name,
contiguous_range_set[-1].end,
)
logging.info('CVO sorting took %s minutes', (time.time() - start_time) / 60)
if num_cvo_records == 0:
return iter([])
logging.info('Transforming call_variants_output to variants.')
independent_variants = _transform_call_variants_output_to_variants(
input_sorted_tfrecord_path=temp_file_name,
sample_name=sample_name,
)
variant_generator = haplotypes.maybe_resolve_conflicting_variants(
independent_variants, qual_filter=_QUAL_FILTER.value
)
logging.info('Processed %s variants.', num_cvo_records)
return variant_generator
def _yield_variants_from_temp_files(
temp_files: Sequence[tempfile._TemporaryFileWrapper],
) -> Iterable[variants_pb2.Variant]:
"""Yields variants from all the temp files in order.
Args:
temp_files: a list of NamedTemporaryFiles objects
Yields:
variants read in order from the given temp files.
"""
for temp_file in temp_files:
for variant in tfrecord.read_tfrecords(
temp_file.name, proto=variants_pb2.Variant
):
yield variant
def _decide_to_use_csi(contigs: Sequence[reference_pb2.ContigInfo]) -> bool:
"""Return True if CSI index is to be used over tabix index format.
If the length of any reference chromosomes exceeds 512M
(here we use 5e8 to keep a safety margin), we will choose csi
as the index format. Otherwise we use tbi as default.
Args:
contigs: list of contigs.
Returns:
A boolean variable indicating if the csi format is to be used or not.
"""
max_chrom_length = max([c.n_bases for c in contigs])
return max_chrom_length > 5e8
def build_index(vcf_file: str, csi: bool = False) -> None:
"""A helper function for indexing VCF files.
Args:
vcf_file: string. Path to the VCF file to be indexed.
csi: bool. If true, index using the CSI format.
"""
if csi:
tabix.build_csi_index(vcf_file, min_shift=14)
else:
tabix.build_index(vcf_file)
def get_cvo_paths(cvo_file_spec: str) -> list[str]:
"""Returns sharded filenames for the `cvo_file_spec` parameter."""
if sharded_file_utils.is_sharded_file_spec(cvo_file_spec):
# Input is already sharded, so dynamic sharding check is disabled.
paths = sharded_file_utils.maybe_generate_sharded_filenames(cvo_file_spec)
else:
# Input is expected to be dynamically sharded.
filename_resolver = cvo_file_spec.replace('.tfrecord.gz', '*')
all_files = sharded_file_utils.glob_list_sharded_file_patterns(
filename_resolver
)
filename_pattern = cvo_file_spec.replace(
'.tfrecord.gz', '@' + str(len(all_files)) + '.tfrecord.gz'
)
paths = sharded_file_utils.maybe_generate_sharded_filenames(
filename_pattern
)
# This check is to make sure all files we glob is exactly the same as the
# paths we create, otherwise we have multiple file patterns.
if sorted(all_files) != sorted(paths):
raise ValueError(
'Found multiple file patterns in input filename space: ',
cvo_file_spec,
)
return paths
def get_first_cvo_record(
paths: Sequence[str],
) -> deepvariant_pb2.CallVariantsOutput | None:
"""Returns the first record from the given paths."""
return dv_utils.get_one_example_from_examples_path(
','.join(paths), proto=deepvariant_pb2.CallVariantsOutput
)
def get_sample_name(cvo_paths: Sequence[str]) -> str:
"""Determines the sample name to be used for the output VCF and gVCF.
We check the following sources to determine the sample name and use the first
name available:
1) CallVariantsOutput
2) nonvariant site TFRecords
3) --sample_name flag
4) default sample name
Args:
cvo_paths: file paths to all CVO files.
Returns:
sample_name used when writing the output VCF and gVCF.
"""
record = get_first_cvo_record(cvo_paths)
gvcf_record = None
if _NONVARIANT_SITE_TFRECORD_PATH.value:
gvcf_record = dv_utils.get_one_example_from_examples_path(
_NONVARIANT_SITE_TFRECORD_PATH.value, proto=variants_pb2.Variant
)
if record is not None:
sample_name = _extract_single_sample_name(record)
logging.info(
'Using sample name from call_variants output. Sample name: %s',
sample_name,
)
if _SAMPLE_NAME.value:
logging.info('--sample_name is set but was not used.')
elif (
_NONVARIANT_SITE_TFRECORD_PATH.value and gvcf_record and gvcf_record.calls
):
sample_name = gvcf_record.calls[0].call_set_name
logging.info(
(
'call_variants output is empty, so using sample name from TFRecords'
' at --nonvariant_site_tfrecord_path. Sample name: %s'
),
sample_name,
)
if _SAMPLE_NAME.value:
logging.info('--sample_name is set but was not used.')
elif _SAMPLE_NAME.value:
sample_name = _SAMPLE_NAME.value
logging.info(
(
'call_variants output and nonvariant TFRecords are empty. Using'
' sample name set with --sample_name. Sample name: %s'
),
sample_name,
)
else:
sample_name = dv_constants.DEFAULT_SAMPLE_NAME
logging.info(
(
'Could not determine sample name and --sample_name is unset. Using'
' the default sample name. Sample name: %s'
),
sample_name,
)
return sample_name
def _merge_phasing_blocks(
input_path: str,
switches_output_path: str,
output_path: str,
) -> None:
"""Reads the phased reads TSV file and loads them into a dictionary.
Args:
input_path: path to the phased reads TSV file.
switches_output_path: path to the switches output TSV file.
output_path: path to the output TSV file.
"""
merger = merge_phased_reads_lib.Merger()
merger.load_from_files(input_path)
merger.merge_reads(switches_output_path)
merger.correct_and_print_read_stats(output_path)
def _load_phasing_info(
switches_output_path: str,
) -> dict[tuple[str, str], int]:
"""Loads the phasing info from the merge_reads output.
Args:
switches_output_path: path to the switches output TSV file.
Returns:
A map from (shard, region) to a boolean indicating whether the phasing
blocks in the shard and region need to be switched.
"""
phasing_info = {}
with open(switches_output_path, 'r') as f:
for line in f:
shard, region, switch_status = line.split('\t')
phasing_info[shard, region] = int(switch_status)
return phasing_info
def run_postprocess_variants_on_region(
output_vcf: str,
output_gvcf: str,
partition: Sequence[range_pb2.Range],
contigs: Sequence[reference_pb2.ContigInfo],
all_cvo_paths: Sequence[str],
header: variants_pb2.VcfHeader,
is_empty: bool,
sample_name: str,
emit_variants_as_tfrecords: bool,
output_tfrecord_file: str,
) -> None:
"""Runs postprocess_variants on the given partition.
If the partition is empty, we process all CVO records. If the partition is not
empty, we process only the CVO records in the partition.
Args:
output_vcf: path to the output VCF file.
output_gvcf: path to the output gVCF file.
partition: a list of nucleus.genomics.v1.Range protos.
contigs: all contigs from ref
all_cvo_paths: paths to all CVO files
header: the VCF header
is_empty: if the partition is empty.
sample_name: the sample name to use for the output VCF and gVCF.
emit_variants_as_tfrecords: if True, emit variants as TFRecords.
output_tfrecord_file: path to the output TFRecord file.
Returns:
None (the output is written to the output_vcf and output_gvcf files).
"""
temp = tempfile.NamedTemporaryFile()
start_time = time.time()
if not is_empty:
variant_generator = process_contiguous_partition(
partition,
contigs,
all_cvo_paths,
temp.name,
sample_name,
)
pon_reader = (
vcf.VcfReader(_PON_FILTERING.value) if _PON_FILTERING.value else None
)
variant_generator = add_pon_filter(variant_generator, pon_reader)
else:
logging.info('call_variants_output is empty. Writing out empty VCF.')
variant_generator = iter([])
logging.info(
'Processing variants (and writing to temporary files) took %s minutes',
(time.time() - start_time) / 60,
)
if emit_variants_as_tfrecords:
tfrecord.write_tfrecords(
variant_generator,
output_tfrecord_file,
compression_type='',
)
else:
emit_variants_to_vcf(
output_vcf,
output_gvcf,
header,
partition,
variant_generator,
)
temp.close()
def emit_variants_to_vcf(
output_vcf: str,
output_gvcf: str,
header: variants_pb2.VcfHeader,
partition: Sequence[range_pb2.Range],
variant_generator: Iterator[variants_pb2.Variant],
variant_tfrecord_path: str = '',
) -> None:
"""Writes variants either from an iterator or a TFRecord file to VCF/gVCF.
Args:
output_vcf: path to the output VCF file.
output_gvcf: path to the output gVCF file.
header: the VCF header
partition: a list of nucleus.genomics.v1.Range protos.
variant_generator: an iterator of variants to write to the VCF file.
variant_tfrecord_path: path to the variant TFRecord file.
"""
start_time = time.time()
if not _NONVARIANT_SITE_TFRECORD_PATH.value:
if _PROCESS_SOMATIC.value:
logging.info('Writing variants to somatic VCF.')
else:
logging.info('Writing variants to VCF.')
if variant_tfrecord_path:
variant_generator = tfrecord.read_tfrecords(
variant_tfrecord_path, variants_pb2.Variant
)
write_variants_to_vcf(
variant_iterable=variant_generator,
output_vcf_path=output_vcf,
header=header,
)
logging.info(
'VCF creation took %s minutes', (time.time() - start_time) / 60
)
else:
tmp_variant_file = None
if not variant_tfrecord_path:
tmp_variant_file = dump_variants_to_temp_file(variant_generator)
variant_tfrecord_path = tmp_variant_file.name
merge_variants.merge_and_write_variants_and_nonvariants(
_ONLY_KEEP_PASS.value,
variant_tfrecord_path,
tfrecord.expanded_paths_if_sharded(
_NONVARIANT_SITE_TFRECORD_PATH.value
),
_REF.value,
output_vcf,
output_gvcf,
header,
partition,
_PROCESS_SOMATIC.value,
)
if tmp_variant_file:
tmp_variant_file.close()
logging.info(
'VCF and gVCF creation took %s minutes.',
(time.time() - start_time) / 60,
)
def stitch_phase_sets(
*,
tfrecord_paths: Sequence[str],
switches_output_path: str,
output_tfrecord_paths: Sequence[str],
) -> None:
"""Stitches the phase sets of the variants in the tfrecord files."""
postprocess_variants_lib.stitch_phase_sets(
tfrecord_paths, switches_output_path, output_tfrecord_paths
)
def _process_partitions_in_parallel(
*,
contigs: Sequence[reference_pb2.ContigInfo],
all_cvo_paths: Sequence[str],
header: variants_pb2.VcfHeader,
is_empty: bool,
sample_name: str,
emit_variants_as_tfrecords: bool,
temp_vcf_files: Sequence[tempfile._TemporaryFileWrapper],
temp_gvcf_files: Sequence[tempfile._TemporaryFileWrapper],
temp_tfrecord_files: Sequence[tempfile._TemporaryFileWrapper],
temp_tfrecord_output_files: Sequence[tempfile._TemporaryFileWrapper],
partitions: Sequence[Sequence[range_pb2.Range]],
num_partitions: int,
):
"""Processes multiple partitions in parallel.
Args:
contigs: all contigs from ref
all_cvo_paths: paths to all CVO files
header: the VCF header
is_empty: if the partition is empty.
sample_name: the sample name to use for the output VCF and gVCF.
emit_variants_as_tfrecords: if True, emit variants as TFRecords.
temp_vcf_files: temporary VCF files.
temp_gvcf_files: temporary gVCF files.
temp_tfrecord_files: temporary TFRecord files.
temp_tfrecord_output_files: temporary TFRecord output files.
partitions: the partitions to process.
num_partitions: the number of partitions.
Returns:
None (the output is written to the output_vcf and output_gvcf files).
"""
logging.info(
'Running postprocess_variants with parallelism using %s CPUs over'
' %s partitions.',
_CPUS.value,
num_partitions,
)
with multiprocessing.Pool(_CPUS.value) as pool:
tasks = []
for task_id in range(num_partitions):
tasks.append(
(
temp_vcf_files[task_id].name,
temp_gvcf_files[task_id].name,
partitions[task_id],
contigs,
all_cvo_paths,
header,
is_empty,
sample_name,
emit_variants_as_tfrecords,
temp_tfrecord_files[task_id].name,
),
)
async_result = pool.starmap_async(run_postprocess_variants_on_region, tasks)
async_result.get()
if emit_variants_as_tfrecords:
stitch_phase_sets(
tfrecord_paths=[t.name for t in temp_tfrecord_files],
switches_output_path=_PHASED_READS_SWITCHES_OUTPUT_PATH.value,
output_tfrecord_paths=[t.name for t in temp_tfrecord_output_files],
)
tasks = []
for task_id in range(num_partitions):
tasks.append((
temp_vcf_files[task_id].name,
temp_gvcf_files[task_id].name,
header,
partitions[task_id],
iter([]),
temp_tfrecord_output_files[task_id].name,
))
async_result = pool.starmap_async(emit_variants_to_vcf, tasks)
async_result.get()
def _process_partitions_sequentially(
*,
contigs: Sequence[reference_pb2.ContigInfo],
all_cvo_paths: Sequence[str],
header: variants_pb2.VcfHeader,
is_empty: bool,
sample_name: str,
emit_variants_as_tfrecords: bool,
temp_vcf_files: Sequence[tempfile._TemporaryFileWrapper],
temp_gvcf_files: Sequence[tempfile._TemporaryFileWrapper],
temp_tfrecord_files: Sequence[tempfile._TemporaryFileWrapper],
temp_tfrecord_output_files: Sequence[tempfile._TemporaryFileWrapper],
partitions: Sequence[Sequence[range_pb2.Range]],
num_partitions: int,
):
"""Processes multiple partitions sequentially.
This mode minimizes the memory usage.
Args:
contigs: all contigs from ref
all_cvo_paths: paths to all CVO files
header: the VCF header
is_empty: if the partition is empty.
sample_name: the sample name to use for the output VCF and gVCF.
emit_variants_as_tfrecords: if true, emit variants as TFRecords.
temp_vcf_files: temporary VCF files.
temp_gvcf_files: temporary gVCF files.
temp_tfrecord_files: temporary TFRecord files.
temp_tfrecord_output_files: temporary TFRecord output files.
partitions: the partitions to process.
num_partitions: the number of partitions.
Returns:
None (the output is written to the output_vcf and output_gvcf files).
"""
logging.info(
'Running postprocess_variants sequentially over %s partitions.',
num_partitions,
)
for task_id in range(num_partitions):
run_postprocess_variants_on_region(
temp_vcf_files[task_id].name,
temp_gvcf_files[task_id].name,
partitions[task_id],
contigs,
all_cvo_paths,
header,
is_empty,
sample_name,
emit_variants_as_tfrecords,
temp_tfrecord_files[task_id].name,
)
if emit_variants_as_tfrecords:
stitch_phase_sets(
tfrecord_paths=[t.name for t in temp_tfrecord_files],
switches_output_path=_PHASED_READS_SWITCHES_OUTPUT_PATH.value,
output_tfrecord_paths=[t.name for t in temp_tfrecord_output_files],
)
for task_id in range(num_partitions):
emit_variants_to_vcf(
temp_vcf_files[task_id].name,
temp_gvcf_files[task_id].name,
header,
partitions[task_id],
iter([]),
temp_tfrecord_output_files[task_id].name,
)
def run_postprocessing_over_multiple_partitions(
*,
contigs: Sequence[reference_pb2.ContigInfo],
all_cvo_paths: Sequence[str],
header: variants_pb2.VcfHeader,
is_empty: bool,
sample_name: str,
) -> None:
"""Runs postprocessing over multiple partitions.
Args:
contigs: all contigs from ref
all_cvo_paths: paths to all CVO files
header: the VCF header
is_empty: if the partition is empty.
sample_name: the sample name to use for the output VCF and gVCF.
Returns:
None (the output is written to the output_vcf and output_gvcf files).
"""
calling_regions = calling_regions_utils.build_calling_regions(
contigs=contigs,
regions_to_include=calling_regions_utils.parse_regions_flag(
_REGIONS.value
),
regions_to_exclude=[],
ref_n_regions=[],
)
num_partitions = max(_NUM_PARTITIONS.value, _CPUS.value)
partitions = calling_regions_utils.partition_calling_regions(
calling_regions, num_partitions=num_partitions
)
temp_vcf_files = [
tempfile.NamedTemporaryFile(suffix='.gz') for _ in partitions
]
temp_gvcf_files = [
tempfile.NamedTemporaryFile(suffix='.gz') for _ in partitions
]
temp_tfrecord_files = [
tempfile.NamedTemporaryFile(suffix='.tfrecord') for _ in partitions
]
temp_tfrecord_output_files = [
tempfile.NamedTemporaryFile(suffix='.tfrecord') for _ in partitions
]
emit_variants_as_tfrecords = bool(_PHASED_READS_INPUT_PATH.value)
if _CPUS.value > 1:
_process_partitions_in_parallel(
contigs=contigs,
all_cvo_paths=all_cvo_paths,
header=header,
is_empty=is_empty,
sample_name=sample_name,
emit_variants_as_tfrecords=emit_variants_as_tfrecords,
temp_vcf_files=temp_vcf_files,
temp_gvcf_files=temp_gvcf_files,
temp_tfrecord_files=temp_tfrecord_files,
temp_tfrecord_output_files=temp_tfrecord_output_files,
partitions=partitions,
num_partitions=num_partitions,
)
else:
_process_partitions_sequentially(
contigs=contigs,
all_cvo_paths=all_cvo_paths,
header=header,
is_empty=is_empty,
sample_name=sample_name,
emit_variants_as_tfrecords=emit_variants_as_tfrecords,
temp_vcf_files=temp_vcf_files,
temp_gvcf_files=temp_gvcf_files,
temp_tfrecord_files=temp_tfrecord_files,
temp_tfrecord_output_files=temp_tfrecord_output_files,
partitions=partitions,
num_partitions=num_partitions,
)
_concat_vcf(_OUTFILE.value, temp_vcf_files)
if _NONVARIANT_SITE_TFRECORD_PATH.value:
_concat_vcf(_GVCF_OUTFILE.value, temp_gvcf_files)
for temp_vcf_file in temp_vcf_files:
temp_vcf_file.close()
for temp_gvcf_file in temp_gvcf_files:
temp_gvcf_file.close()
for temp_tfrecord_file in temp_tfrecord_files:
temp_tfrecord_file.close()
for temp_tfrecord_output_file in temp_tfrecord_output_files:
temp_tfrecord_output_file.close()
def run_postprocessing_without_partitioning(
*,
contigs: Sequence[reference_pb2.ContigInfo],
all_cvo_paths: Sequence[str],
header: variants_pb2.VcfHeader,
is_empty: bool,
sample_name: str,
) -> None:
"""Runs postprocessing without any partitioning.
Args:
contigs: all contigs from ref
all_cvo_paths: paths to all CVO files
header: the VCF header
is_empty: if the partition is empty.
sample_name: the sample name to use for the output VCF and gVCF.
Returns:
None (the output is written to the output_vcf and output_gvcf files).
"""
logging.info(
'Running postprocess_variants without parallelism or partitions.'
)
emit_variants_as_tfrecords = bool(_PHASED_READS_INPUT_PATH.value)
tmp_tfrecord_file = None
tmp_tfrecord_file_name = ''
if emit_variants_as_tfrecords:
tmp_tfrecord_file = tempfile.NamedTemporaryFile(suffix='.tfrecord')
tmp_tfrecord_file_name = tmp_tfrecord_file.name
run_postprocess_variants_on_region(
_OUTFILE.value,
_GVCF_OUTFILE.value,
[],
contigs,
all_cvo_paths,
header,
is_empty,
sample_name,
emit_variants_as_tfrecords,
tmp_tfrecord_file_name,
)
if emit_variants_as_tfrecords:
output_tfrecord_file = tempfile.NamedTemporaryFile(suffix='.tfrecord')
stitch_phase_sets(
tfrecord_paths=[tmp_tfrecord_file_name],
switches_output_path=_PHASED_READS_SWITCHES_OUTPUT_PATH.value,
output_tfrecord_paths=[output_tfrecord_file.name],
)
emit_variants_to_vcf(
_OUTFILE.value,
_GVCF_OUTFILE.value,
header,
[],
iter([]),
output_tfrecord_file.name,
)
output_tfrecord_file.close()
if tmp_tfrecord_file:
tmp_tfrecord_file.close()
def apply_flags_for_postprocessing(flags_obj):
"""Read flags for postprocessing from the model.example_info.json file.
Args:
flags_obj: The flag values object.
"""
if not flags_obj.checkpoint_json:
return
logging.info(
'Reading flags_for_postprocessing from %s', flags_obj.checkpoint_json
)
with tf.io.gfile.GFile(flags_obj.checkpoint_json, 'r') as fin:
flags_map = json.load(fin).get('flags_for_postprocessing', {})
if flags_map:
logging.info(
'Flags for postprocessing:\n%s',
'\n'.join([f'{k}: {v}' for k, v in flags_map.items()]),
)
for flag_name, flag_value in flags_map.items():
if flag_name not in flags_obj:
logging.warning(
'Flag "%s" from json is not defined as an application flag.',
flag_name,
)
continue
flag = flags_obj[flag_name]
if flag.present:
if flag.value != flag_value:
logging.warning(
'Flag %s is specified in model.example_info.json [%s] but '
'overridden by command line with value [%s]',
flag_name,
flag_value,
flag.value,
)
continue
flag.value = flag_value
def main(argv=()):
with errors.clean_commandline_error_exit():
apply_flags_for_postprocessing(flags.FLAGS)
if len(argv) > 1:
errors.log_and_raise(
'Command line parsing failure: postprocess_variants does not accept '
'positional arguments but some are present on the command line: '
'"{}".'.format(str(argv)),
errors.CommandLineError,
)
del argv # Unused.
if (not _NONVARIANT_SITE_TFRECORD_PATH.value) != (not _GVCF_OUTFILE.value):
errors.log_and_raise(
(
'gVCF creation requires both nonvariant_site_tfrecord_path and '
'gvcf_outfile flags to be set.'
),
errors.CommandLineError,
)
if (
_USE_MULTIALLELIC_MODEL.value
and _DEBUG_OUTPUT_ALL_CANDIDATES.value == 'ALT'
):
errors.log_and_raise(
(
'debug_output_all_candidates=ALT is incompatible with the '
'multiallelic model. Use INFO instead.'
),
errors.CommandLineError,
)
if _NUM_PARTITIONS.value > 0 and _NUM_PARTITIONS.value < _CPUS.value:
logging.warning(
'--num_partitions is less than --cpus. Setting --num_partitions to'
' --cpus=%s',
_CPUS.value,
)
proto_utils.uses_fast_cpp_protos_or_die()
logging_level.set_from_flag()
fasta_reader = pysam.FastaFile(
filename=_pysam_resolve_file_path(_REF.value)
)
contigs = []
for reference_index in range(fasta_reader.nreferences):
contigs.append(
reference_pb2.ContigInfo(
name=fasta_reader.references[reference_index],
n_bases=fasta_reader.lengths[reference_index],
pos_in_fasta=reference_index,
)
)
cvo_paths = get_cvo_paths(_INFILE.value)
small_model_cvo_paths = []
if _SMALL_MODEL_CVO_RECORDS.value:
small_model_cvo_paths = get_cvo_paths(_SMALL_MODEL_CVO_RECORDS.value)
all_cvo_paths = cvo_paths + small_model_cvo_paths
sample_name = get_sample_name(all_cvo_paths)
header = dv_vcf_constants.deepvariant_header(
contigs=contigs,
sample_names=[sample_name],
add_info_candidates=_DEBUG_OUTPUT_ALL_CANDIDATES.value == 'INFO',
include_model_id=_SMALL_MODEL_CVO_RECORDS.value is not None,
include_somatic_fields=_PROCESS_SOMATIC.value,
)
if _PROCESS_SOMATIC.value:
header.filters.append(
variants_pb2.VcfFilterInfo(
id=dv_vcf_constants.DEEP_VARIANT_GERMLINE,
description='Non somatic variants',
)
)
if _PON_FILTERING.value:
if not _PROCESS_SOMATIC.value:
raise ValueError(
'PON filtering is only supported for somatic variant calling.'
)
header.filters.append(
variants_pb2.VcfFilterInfo(
id=dv_vcf_constants.DEEP_VARIANT_PON,
description='Filtered by Panel of Normals (PON)',
)
)
if _PHASED_READS_INPUT_PATH.value:
logging.info(
'Attempting to merge phasing blocks from %s',
_PHASED_READS_INPUT_PATH.value,
)
logging.info(
'Writing switches to %s',
_PHASED_READS_SWITCHES_OUTPUT_PATH.value,
)
logging.info(
'Writing corrected reads to %s',
_PHASED_READS_CORRECTED_OUTPUT_PATH.value,
)
_merge_phasing_blocks(
_PHASED_READS_INPUT_PATH.value,
_PHASED_READS_SWITCHES_OUTPUT_PATH.value,
_PHASED_READS_CORRECTED_OUTPUT_PATH.value,
)
is_empty = get_first_cvo_record(all_cvo_paths) is None
# Run sequentially in the absence of multiple CPUs or partitions.
if _CPUS.value < 1 and _NUM_PARTITIONS.value < 1:
run_postprocessing_without_partitioning(
contigs=contigs,
all_cvo_paths=all_cvo_paths,
header=header,
is_empty=is_empty,
sample_name=sample_name,
)
else:
run_postprocessing_over_multiple_partitions(
contigs=contigs,
all_cvo_paths=all_cvo_paths,
header=header,
is_empty=is_empty,
sample_name=sample_name,
)
start_time = time.time()
use_csi = _decide_to_use_csi(contigs)
if str(_OUTFILE.value).endswith('.gz'):
build_index(_OUTFILE.value, use_csi)
if _NONVARIANT_SITE_TFRECORD_PATH.value and str(
_GVCF_OUTFILE.value
).endswith('.gz'):
build_index(_GVCF_OUTFILE.value, use_csi)
logging.info(
'Indexing VCF and gVCF took %s minutes.',
(time.time() - start_time) / 60,
)
if __name__ == '__main__':
flags.mark_flags_as_required(['infile', 'outfile', 'ref'])
logging.set_verbosity(logging.INFO)
logging.get_absl_logger().setLevel(logging.INFO)
app.run(main)

postprocess_variants.py at commit 45f2627, under BSD-3-Clause · at the source

Overview

Authors: Andrew B Munkacsi^1,2, Gabriella E Conway¹, Keri Multerer¹, Michael D Jackson¹, Lauren Csaki³, Tamayanthi Rajakumar¹, Jeffrey P Sheridan¹, Eliatan Niktab¹, Wuh-Liang Hwu^4,5, Yin-Hsiu Chien⁵, Mark Walterfang⁶, Cintya Del Rio Hernandez¹, Robin B Chan⁷, Bowen Zhou⁸, Renu Nandakumar⁷, Peixang Zhang³, Gilbert Di Paolo⁹, Robert A Maue¹⁰, Karen Reue³, Stephen L Sturley¹¹

ORCID iDs: Andrew B Munkacsi

Institutions: Victoria University of Wellington (New Zealand); University of California, Los Angeles (United States); China Medical University (Taiwan); China Medical University Hospital (Taiwan); National Taiwan University Hospital (Taiwan); The Royal Melbourne Hospital (Australia); Columbia University Irving Medical Center (United States); Columbia University (United States); Case Western Reserve University (United States); Denali Therapeutics (United States) (United States); Dartmouth College (United States); Barnard College (United States)

Journal: iScience, volume 29, issue 6, article 116121

Dates: received 26 November 2024; accepted 11 May 2026; published online 28 May 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1016/j.isci.2026.116121 · PMID 42256287 · PMCID PMC13233811 · OpenAlex W7162675976

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), mouse (organism)

Methods: Statistics, Smoothing, state filtering, decompositions, Preprocessing, Graphs

Keywords: Disease, Neurogenetics, Genomics, Lipidomics, Model organism

Topic: Lysosomal Storage Disorders Research (Physiology, Medicine), according to OpenAlex

Funding: National Institutes of Health (T32 HL07343); Dartmouth Geisel School of Medicine (200128); National Institute of Diabetes and Digestive and Kidney Diseases (DK54320); Charles H. Revson Foundation (T32 HL07343); NIDDK NIH HHS (R01 DK128898); National Niemann Pick Disease Foundation Inc; Victoria University of Wellington; Research for Life; Ara Parseghian Medical Research Foundation

Citations: not cited yet (Europe PMC); 94 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

google/deepvariant

License: BSD-3-Clause

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 45f2627504c59785ea2b88d0256a2ec347bce7b4, 18 March 2026

Languages: Python (205), C++ (170), C/C++ (98), Shell (18), Jupyter (1), C (1)

Size: 924 files, 493 scripts

Software Heritage: archived

Found in: the resources table

Holds: README, license file, environment (Dockerfile, Dockerfile.deepsomatic, Dockerfile.deeptrio, Dockerfile.pangenome_aware_deepvariant, Dockerfile.tpu-train, setup.py), tests, documentation, 1 notebook

Not found: CITATION.cff, continuous integration

Tools: TensorFlow (44 files), NumPy (28 files), pandas (6 files), pysam (4 files), Keras (3 files), Pillow (2 files), BCFtools (1 file), JAX (1 file), SAMtools (1 file), scikit-learn (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

495 files

__init__.py, Python, 29 lines
build-prereq.sh, Shell, 147 lines
build_and_test.sh, Shell, 87 lines
build_release_binaries.sh, Shell, 191 lines
deeptrio/dt_constants.py, Python, 60 lines
deeptrio/make_examples.py, Python, 447 lines
deeptrio/make_examples_test.py, Python, 1,521 lines
deeptrio/testdata.py, Python, 234 lines
deepvariant/__init__.py, Python, 29 lines
deepvariant/allele_frequency.py, Python, 421 lines
deepvariant/allele_frequency_test.py, Python, 570 lines
deepvariant/allelecounter.cc, C++, 1,011 lines
deepvariant/allelecounter.h, C/C++, 518 lines
deepvariant/allelecounter_test.cc, C++, 1,582 lines
deepvariant/alt_aligned_pileup_lib.cc, C++, 317 lines
deepvariant/alt_aligned_pileup_lib.h, C/C++, 96 lines
deepvariant/alt_aligned_pileup_lib_test.cc, C++, 501 lines
deepvariant/call_variants.py, Python, 1,116 lines
deepvariant/call_variants_test.py, Python, 400 lines
deepvariant/calling_regions_utils.py, Python, 162 lines
deepvariant/calling_regions_utils_test.py, Python, 327 lines
deepvariant/channels/allele_frequency_channel.cc, C++, 124 lines
deepvariant/channels/allele_frequency_channel.h, C/C++, 83 lines
deepvariant/channels/allele_sample_probability_channel.cc, C++, 103 lines
deepvariant/channels/allele_sample_probability_channel.h, C/C++, 71 lines
deepvariant/channels/avg_base_quality_channel.cc, C++, 98 lines
deepvariant/channels/avg_base_quality_channel.h, C/C++, 87 lines
deepvariant/channels/base_6ma_channel.cc, C++, 101 lines
deepvariant/channels/base_6ma_channel.h, C/C++, 85 lines
deepvariant/channels/base_differs_from_ref_channel.cc, C++, 69 lines
deepvariant/channels/base_differs_from_ref_channel.h, C/C++, 74 lines
deepvariant/channels/base_methylation_channel.cc, C++, 101 lines
deepvariant/channels/base_methylation_channel.h, C/C++, 84 lines
deepvariant/channels/base_quality_channel.cc, C++, 70 lines
deepvariant/channels/base_quality_channel.h, C/C++, 74 lines
deepvariant/channels/blank_channel.cc, C++, 58 lines
deepvariant/channels/blank_channel.h, C/C++, 71 lines
deepvariant/channels/channel.cc, C++, 47 lines
deepvariant/channels/channel.h, C/C++, 86 lines
deepvariant/channels/channel_utils.cc, C++, 51 lines
deepvariant/channels/channel_utils.h, C/C++, 60 lines
deepvariant/channels/gap_compressed_identity_channel.cc, C++, 116 lines
deepvariant/channels/gap_compressed_identity_channel.h, C/C++, 85 lines
deepvariant/channels/gc_content_channel.cc, C++, 102 lines
deepvariant/channels/gc_content_channel.h, C/C++, 84 lines
deepvariant/channels/haplotype_tag_channel.cc, C++, 112 lines
deepvariant/channels/haplotype_tag_channel.h, C/C++, 85 lines
deepvariant/channels/homopolymer_deletion_quality_channel.cc, C++, 74 lines
deepvariant/channels/homopolymer_deletion_quality_channel.h, C/C++, 76 lines
deepvariant/channels/homopolymer_indel_quality_channel.cc, C++, 187 lines
deepvariant/channels/homopolymer_indel_quality_channel.h, C/C++, 73 lines
deepvariant/channels/homopolymer_insertion_quality_channel.cc, C++, 74 lines
deepvariant/channels/homopolymer_insertion_quality_channel.h, C/C++, 76 lines
deepvariant/channels/homopolymer_weighted_channel.cc, C++, 123 lines
deepvariant/channels/homopolymer_weighted_channel.h, C/C++, 91 lines
deepvariant/channels/identity_channel.cc, C++, 107 lines
deepvariant/channels/identity_channel.h, C/C++, 86 lines
deepvariant/channels/insert_size_channel.cc, C++, 92 lines
deepvariant/channels/insert_size_channel.h, C/C++, 85 lines
deepvariant/channels/inter_homopolymer_insertion_quality_channel.cc, C++, 129 lines
deepvariant/channels/inter_homopolymer_insertion_quality_channel.h, C/C++, 83 lines
deepvariant/channels/is_homopolymer_channel.cc, C++, 114 lines
deepvariant/channels/is_homopolymer_channel.h, C/C++, 88 lines
deepvariant/channels/mapping_quality_channel.cc, C++, 70 lines
deepvariant/channels/mapping_quality_channel.h, C/C++, 73 lines
deepvariant/channels/read_base_channel.cc, C++, 77 lines
deepvariant/channels/read_base_channel.h, C/C++, 74 lines
deepvariant/channels/read_mapping_percent_channel.cc, C++, 102 lines
deepvariant/channels/read_mapping_percent_channel.h, C/C++, 87 lines
deepvariant/channels/read_supports_variant_channel.cc, C++, 120 lines
deepvariant/channels/read_supports_variant_channel.h, C/C++, 84 lines
deepvariant/channels/read_supports_variant_fuzzy_channel.cc, C++, 312 lines
deepvariant/channels/read_supports_variant_fuzzy_channel.h, C/C++, 81 lines
deepvariant/channels/read_supports_variant_fuzzy_channel_test.cc, C++, 163 lines
deepvariant/channels/strand_channel.cc, C++, 64 lines
deepvariant/channels/strand_channel.h, C/C++, 73 lines
deepvariant/channels/supplementary_alignment_channel.cc, C++, 68 lines
deepvariant/channels/supplementary_alignment_channel.h, C/C++, 69 lines
deepvariant/convert_to_saved_model.py, Python, 142 lines
deepvariant/dashboard_utils.py, Python, 176 lines
deepvariant/dashboard_utils_test.py, Python, 92 lines
deepvariant/data_providers.py, Python, 314 lines
deepvariant/data_providers_test.py, Python, 188 lines
deepvariant/direct_phasing.cc, C++, 1,006 lines
deepvariant/direct_phasing.h, C/C++, 381 lines
deepvariant/direct_phasing_test.cc, C++, 1,240 lines
deepvariant/distribution_functor.h, C/C++, 338 lines
deepvariant/distribution_functor_test.cc, C++, 210 lines
deepvariant/dv_config.py, Python, 562 lines
deepvariant/dv_constants.py, Python, 208 lines
deepvariant/dv_utils.py, Python, 445 lines
deepvariant/dv_utils_test.py, Python, 247 lines
deepvariant/dv_utils_using_clif.py, Python, 135 lines
deepvariant/dv_utils_using_clif_test.py, Python, 130 lines
deepvariant/dv_vcf_constants.py, Python, 227 lines
deepvariant/dv_vcf_constants_test.py, Python, 104 lines
deepvariant/environment_tests/env_smoke_test.py, Python, 51 lines
deepvariant/environment_tests/protobuf_implementation_test.py, Python, 45 lines
deepvariant/exclude_contigs.py, Python, 3,444 lines
deepvariant/exclude_contigs_test.py, Python, 48 lines
deepvariant/fast_pipeline.cc, C++, 268 lines
deepvariant/fast_pipeline.h, C/C++, 87 lines
deepvariant/fast_pipeline_utils.h, C/C++, 67 lines
deepvariant/haplotypes.py, Python, 539 lines
deepvariant/haplotypes_test.py, Python, 877 lines
deepvariant/keras_modeling.py, Python, 495 lines
deepvariant/keras_modeling_test.py, Python, 253 lines
deepvariant/labeler/__init__.py, Python, 29 lines
deepvariant/labeler/combined_labeler.py, Python, 125 lines
deepvariant/labeler/combined_labeler_test.py, Python, 137 lines
deepvariant/labeler/compare_labelers.py, Python, 136 lines
deepvariant/labeler/customized_classes_labeler.py, Python, 213 lines
deepvariant/labeler/customized_classes_labeler_test.py, Python, 310 lines
deepvariant/labeler/haplotype_labeler.py, Python, 1,309 lines
deepvariant/labeler/haplotype_labeler_test.py, Python, 2,204 lines
deepvariant/labeler/labeled_examples_to_vcf.py, Python, 234 lines
deepvariant/labeler/labeled_examples_to_vcf_test.py, Python, 97 lines
deepvariant/labeler/positional_labeler.py, Python, 249 lines
deepvariant/labeler/positional_labeler_test.py, Python, 243 lines
deepvariant/labeler/soft_labeler.py, Python, 1,466 lines
deepvariant/labeler/soft_labeler_test.py, Python, 167 lines
deepvariant/labeler/variant_labeler.py, Python, 290 lines
deepvariant/labeler/variant_labeler_test.py, Python, 299 lines
deepvariant/load_gbz_into_shared_memory.py, Python, 161 lines
deepvariant/logging_level.py, Python, 65 lines
deepvariant/make_examples.py, Python, 239 lines
deepvariant/make_examples_core.py, Python, 3,908 lines
deepvariant/make_examples_core_test.py, Python, 1,910 lines
deepvariant/make_examples_native.cc, C++, 895 lines
deepvariant/make_examples_native.h, C/C++, 321 lines
deepvariant/make_examples_native_test.cc, C++, 838 lines
deepvariant/make_examples_options.py, Python, 1,539 lines
deepvariant/make_examples_pangenome_aware_dv.py, Python, 398 lines
deepvariant/make_examples_somatic.py, Python, 339 lines
deepvariant/make_examples_somatic_test.py, Python, 156 lines
deepvariant/make_examples_test.py, Python, 1,503 lines
deepvariant/merge_phased_reads.cc, C++, 362 lines
deepvariant/merge_phased_reads.h, C/C++, 183 lines
deepvariant/merge_phased_reads_main.cc, C++, 64 lines
deepvariant/merge_phased_reads_test.cc, C++, 364 lines
deepvariant/methylation_aware_phasing.cc, C++, 483 lines, 1 match
deepvariant/methylation_aware_phasing.h, C/C++, 108 lines, 1 match
deepvariant/methylation_aware_phasing_test.cc, C++, 231 lines
deepvariant/metrics.py, Python, 106 lines
deepvariant/multisample_make_examples.py, Python, 270 lines
deepvariant/pileup_channel_lib.cc, C++, 526 lines
deepvariant/pileup_channel_lib.h, C/C++, 212 lines
deepvariant/pileup_channel_lib_test.cc, C++, 955 lines
deepvariant/pileup_image.py, Python, 139 lines
deepvariant/pileup_image_native.cc, C++, 531 lines
deepvariant/pileup_image_native.h, C/C++, 341 lines
deepvariant/pileup_image_native_test.cc, C++, 973 lines
deepvariant/pileup_image_test.py, Python, 785 lines
deepvariant/postprocess_variants.cc, C++, 327 lines
deepvariant/postprocess_variants.h, C/C++, 101 lines
deepvariant/postprocess_variants.py, Python, 2,378 lines, 1 match
deepvariant/postprocess_variants_test.cc, C++, 403 lines
deepvariant/postprocess_variants_test.py, Python, 2,284 lines
deepvariant/python/allelecounter_pybind.cc, C++, 93 lines
deepvariant/python/allelecounter_wrap_test.py, Python, 66 lines
deepvariant/python/clif_converters.cc, C++, 92 lines
deepvariant/python/clif_converters.h, C/C++, 54 lines
deepvariant/python/direct_phasing_pybind.cc, C++, 65 lines
deepvariant/python/make_examples_native_pybind.cc, C++, 109 lines
deepvariant/python/merge_phased_reads_pybind.cc, C++, 58 lines
deepvariant/python/methylation_aware_phasing_pybind.cc, C++, 57 lines
deepvariant/python/pileup_image_native_pybind.cc, C++, 113 lines
deepvariant/python/postprocess_variants_pybind.cc, C++, 58 lines
deepvariant/python/variant_calling_multisample_pybind.cc, C++, 85 lines
deepvariant/python/variant_calling_multisample_wrap_test.py, Python, 94 lines
deepvariant/python/variant_calling_pybind.cc, C++, 75 lines
deepvariant/python/variant_calling_wrap_test.py, Python, 93 lines
deepvariant/realigner/__init__.py, Python, 29 lines
deepvariant/realigner/debruijn_graph.cc, C++, 487 lines
deepvariant/realigner/debruijn_graph.h, C/C++, 192 lines
deepvariant/realigner/fast_pass_aligner.cc, C++, 993 lines
deepvariant/realigner/fast_pass_aligner.h, C/C++, 423 lines
deepvariant/realigner/fast_pass_aligner_test.cc, C++, 1,249 lines
deepvariant/realigner/python/debruijn_graph_pybind.cc, C++, 56 lines
deepvariant/realigner/python/debruijn_graph_wrap_test.py, Python, 362 lines
deepvariant/realigner/python/fast_pass_aligner_pybind.cc, C++, 80 lines
deepvariant/realigner/python/ssw_misc_test.py, Python, 87 lines
deepvariant/realigner/python/ssw_pybind.cc, C++, 115 lines
deepvariant/realigner/python/ssw_wrap_test.py, Python, 78 lines
deepvariant/realigner/python/window_selector_pybind.cc, C++, 57 lines
deepvariant/realigner/realigner.py, Python, 1,062 lines
deepvariant/realigner/realigner_test.py, Python, 970 lines
deepvariant/realigner/ssw.cc, C++, 78 lines
deepvariant/realigner/ssw.h, C/C++, 80 lines
deepvariant/realigner/ssw_test.cc, C++, 66 lines
deepvariant/realigner/utils.py, Python, 43 lines
deepvariant/realigner/window_selector.cc, C++, 212 lines
deepvariant/realigner/window_selector.h, C/C++, 73 lines
deepvariant/realigner/window_selector.py, Python, 238 lines
deepvariant/realigner/window_selector_test.py, Python, 600 lines
deepvariant/resources.py, Python, 164 lines
deepvariant/resources_test.py, Python, 153 lines
deepvariant/runtime_by_region_vis.py, Python, 680 lines
deepvariant/runtime_by_region_vis_test.py, Python, 168 lines
deepvariant/sample.py, Python, 75 lines
deepvariant/sampling_util.h, C/C++, 161 lines
deepvariant/sampling_util_test.cc, C++, 156 lines
deepvariant/show_examples.py, Python, 584 lines
deepvariant/show_examples_test.py, Python, 242 lines
deepvariant/small_model/__init__.py, Python, 29 lines
deepvariant/small_model/inference.py, Python, 194 lines
deepvariant/small_model/inference_test.py, Python, 260 lines
deepvariant/small_model/keras_config.py, Python, 194 lines
deepvariant/small_model/make_small_model_examples.py, Python, 823 lines
deepvariant/small_model/make_small_model_examples_test.py, Python, 1,299 lines
deepvariant/small_model/small_model_config.py, Python, 132 lines
deepvariant/small_model/train_utils_test.py, Python, 180 lines
deepvariant/stream_examples.cc, C++, 179 lines
deepvariant/stream_examples.h, C/C++, 114 lines
deepvariant/stream_examples_kernel.cc, C++, 341 lines
deepvariant/stream_examples_ops.cc, C++, 66 lines
deepvariant/testdata.py, Python, 375 lines
deepvariant/testing_utils.cc, C++, 170 lines
deepvariant/testing_utils.h, C/C++, 95 lines
deepvariant/train.py, Python, 773 lines
deepvariant/train_test.py, Python, 111 lines
deepvariant/utils.cc, C++, 102 lines
deepvariant/utils.h, C/C++, 65 lines
deepvariant/utils_test.cc, C++, 71 lines
deepvariant/variant_caller.py, Python, 486 lines
deepvariant/variant_caller_test.py, Python, 578 lines
deepvariant/variant_caller_trio_test.py, Python, 544 lines
deepvariant/variant_calling.cc, C++, 720 lines
deepvariant/variant_calling.h, C/C++, 251 lines
deepvariant/variant_calling_multisample.cc, C++, 1,581 lines
deepvariant/variant_calling_multisample.h, C/C++, 454 lines
deepvariant/variant_calling_multisample_somatic_test.cc, C++, 324 lines
deepvariant/variant_calling_multisample_test.cc, C++, 1,339 lines
deepvariant/variant_calling_multisample_trio_test.cc, C++, 1,138 lines
deepvariant/variant_calling_test.cc, C++, 1,475 lines
deepvariant/vcf_candidate_importer.py, Python, 78 lines
deepvariant/vcf_candidate_importer_test.py, Python, 161 lines
deepvariant/vcf_stats.py, Python, 394 lines
deepvariant/vcf_stats_report.py, Python, 103 lines
deepvariant/vcf_stats_test.py, Python, 431 lines
deepvariant/vcf_stats_vis.py, Python, 542 lines
deepvariant/vcf_stats_vis_test.py, Python, 236 lines
deepvariant/vendor/timer.py, Python, 259 lines
deepvariant/vendor/timer_test.py, Python, 43 lines
deepvariant/very_sensitive_caller.py, Python, 75 lines
deepvariant/very_sensitive_caller_test.py, Python, 130 lines
deepvariant/very_sensitive_caller_trio_test.py, Python, 127 lines
docs/visualizing_examples.ipynb, Jupyter, 144 lines
run-prereq.sh, Shell, 306 lines
scripts/create_golden.sh, Shell, 500 lines
scripts/create_golden_deeptrio.sh, Shell, 431 lines
scripts/docker_utils.sh, Shell, 191 lines
scripts/inference_deepsomatic.sh, Shell, 909 lines
scripts/inference_deeptrio.sh, Shell, 949 lines
scripts/inference_deepvariant.sh, Shell, 943 lines
scripts/inference_utils.sh, Shell, 68 lines
scripts/install_nvidia_docker.sh, Shell, 88 lines
scripts/install_singularity.sh, Shell, 68 lines
scripts/run_deepsomatic.py, Python, 964 lines
scripts/run_deepsomatic_test.py, Python, 539 lines
scripts/run_deeptrio.py, Python, 1,325 lines
scripts/run_deeptrio_test.py, Python, 1,061 lines
scripts/run_deepvariant.py, Python, 920 lines
scripts/run_oracle_inference.py, Python, 488 lines
scripts/run_oracle_inference_test.py, Python, 116 lines
scripts/run_pangenome_aware_deepvariant.py, Python, 861 lines
scripts/run_pangenome_aware_deepvariant_test.py, Python, 458 lines
settings.sh, Shell, 130 lines
setup.py, Python, 88 lines
third_party/nucleus/__init__.py, Python, 29 lines
third_party/nucleus/core/__init__.py, Python, 28 lines
third_party/nucleus/core/macros.h, C/C++, 48 lines
third_party/nucleus/core/python/__init__.py, Python, 28 lines
third_party/nucleus/core/python/statusor_examples_pybind.cc, C++, 65 lines
third_party/nucleus/core/python/statusor_examples_test.py, Python, 96 lines
third_party/nucleus/core/python/type_caster_nucleus_status.h, C/C++, 67 lines
third_party/nucleus/core/python/type_caster_nucleus_statusor.h, C/C++, 91 lines
third_party/nucleus/core/status.cc, C++, 201 lines
third_party/nucleus/core/status.h, C/C++, 271 lines
third_party/nucleus/core/status_matchers.h, C/C++, 77 lines
third_party/nucleus/core/statusor.h, C/C++, 274 lines
third_party/nucleus/core/statusor_examples.h, C/C++, 101 lines
third_party/nucleus/core/statusor_examples_test.py, Python, 86 lines
third_party/nucleus/core/statusor_test.cc, C++, 89 lines
third_party/nucleus/io/__init__.py, Python, 28 lines
third_party/nucleus/io/bed.py, Python, 169 lines
third_party/nucleus/io/bed_reader.cc, C++, 369 lines
third_party/nucleus/io/bed_reader.h, C/C++, 158 lines
third_party/nucleus/io/bed_reader_test.cc, C++, 173 lines
third_party/nucleus/io/bed_test.py, Python, 148 lines
third_party/nucleus/io/bed_writer.cc, C++, 123 lines
third_party/nucleus/io/bed_writer.h, C/C++, 106 lines
third_party/nucleus/io/bed_writer_test.cc, C++, 142 lines
third_party/nucleus/io/bedgraph.py, Python, 155 lines
third_party/nucleus/io/bedgraph_reader.cc, C++, 163 lines
third_party/nucleus/io/bedgraph_reader.h, C/C++, 103 lines
third_party/nucleus/io/bedgraph_reader_test.cc, C++, 89 lines
third_party/nucleus/io/bedgraph_test.py, Python, 77 lines
third_party/nucleus/io/bedgraph_writer.cc, C++, 86 lines
third_party/nucleus/io/bedgraph_writer.h, C/C++, 99 lines
third_party/nucleus/io/bedgraph_writer_test.cc, C++, 135 lines
third_party/nucleus/io/clif_postproc.py, Python, 163 lines
third_party/nucleus/io/converter.py, Python, 202 lines
third_party/nucleus/io/converter_test.py, Python, 102 lines
third_party/nucleus/io/example_writer.cc, C++, 175 lines
third_party/nucleus/io/example_writer.h, C/C++, 71 lines
third_party/nucleus/io/example_writer_test.cc, C++, 67 lines
third_party/nucleus/io/fasta.py, Python, 266 lines
third_party/nucleus/io/fasta_test.py, Python, 244 lines
third_party/nucleus/io/fastq.py, Python, 155 lines
third_party/nucleus/io/fastq_reader.cc, C++, 205 lines
third_party/nucleus/io/fastq_reader.h, C/C++, 122 lines
third_party/nucleus/io/fastq_reader_test.cc, C++, 112 lines
third_party/nucleus/io/fastq_test.py, Python, 89 lines
third_party/nucleus/io/fastq_writer.cc, C++, 103 lines
third_party/nucleus/io/fastq_writer.h, C/C++, 100 lines
third_party/nucleus/io/fastq_writer_test.cc, C++, 119 lines
third_party/nucleus/io/gbz_reader.cc, C++, 466 lines
third_party/nucleus/io/gbz_reader.h, C/C++, 162 lines
third_party/nucleus/io/genomics_io_noplugin_test.py, Python, 51 lines
third_party/nucleus/io/genomics_reader.py, Python, 255 lines
third_party/nucleus/io/genomics_reader_test.py, Python, 120 lines
third_party/nucleus/io/genomics_writer.py, Python, 214 lines
third_party/nucleus/io/gff.py, Python, 160 lines
third_party/nucleus/io/gff_reader.cc, C++, 324 lines
third_party/nucleus/io/gff_reader.h, C/C++, 119 lines
third_party/nucleus/io/gff_reader_test.cc, C++, 117 lines
third_party/nucleus/io/gff_test.py, Python, 102 lines
third_party/nucleus/io/gff_writer.cc, C++, 180 lines
third_party/nucleus/io/gff_writer.h, C/C++, 108 lines
third_party/nucleus/io/gff_writer_test.cc, C++, 127 lines
third_party/nucleus/io/gfile.cc, C++, 115 lines
third_party/nucleus/io/gfile.h, C/C++, 89 lines
third_party/nucleus/io/gfile.py, Python, 87 lines
third_party/nucleus/io/gfile_test.py, Python, 76 lines
third_party/nucleus/io/hts_path.cc, C++, 80 lines
third_party/nucleus/io/hts_path.h, C/C++, 59 lines
third_party/nucleus/io/hts_test.cc, C++, 73 lines
third_party/nucleus/io/hts_verbose.cc, C++, 41 lines
third_party/nucleus/io/hts_verbose.h, C/C++, 45 lines
third_party/nucleus/io/hts_verbose_test.cc, C++, 59 lines
third_party/nucleus/io/merge_variants.cc, C++, 233 lines
third_party/nucleus/io/merge_variants.h, C/C++, 80 lines
third_party/nucleus/io/merge_variants_test.cc, C++, 235 lines
third_party/nucleus/io/python/__init__.py, Python, 28 lines
third_party/nucleus/io/python/bed_reader_pybind.cc, C++, 82 lines
third_party/nucleus/io/python/bed_reader_wrap_test.py, Python, 127 lines
third_party/nucleus/io/python/bed_writer_pybind.cc, C++, 60 lines
third_party/nucleus/io/python/bed_writer_wrap_test.py, Python, 103 lines
third_party/nucleus/io/python/bedgraph_reader_pybind.cc, C++, 70 lines
third_party/nucleus/io/python/bedgraph_writer_pybind.cc, C++, 56 lines
third_party/nucleus/io/python/fastq_reader_pybind.cc, C++, 70 lines
third_party/nucleus/io/python/fastq_reader_wrap_test.py, Python, 84 lines
third_party/nucleus/io/python/fastq_writer_pybind.cc, C++, 58 lines
third_party/nucleus/io/python/fastq_writer_wrap_test.py, Python, 142 lines
third_party/nucleus/io/python/gbz_reader_pybind.cc, C++, 64 lines
third_party/nucleus/io/python/gff_reader_pybind.cc, C++, 72 lines
third_party/nucleus/io/python/gff_reader_wrap_test.py, Python, 94 lines
third_party/nucleus/io/python/gff_writer_pybind.cc, C++, 56 lines
third_party/nucleus/io/python/gff_writer_wrap_test.py, Python, 99 lines
third_party/nucleus/io/python/gfile_pybind.cc, C++, 75 lines
third_party/nucleus/io/python/hts_verbose_pybind.cc, C++, 61 lines
third_party/nucleus/io/python/hts_verbose_test.py, Python, 58 lines
third_party/nucleus/io/python/merge_variants_pybind.cc, C++, 65 lines
third_party/nucleus/io/python/reference_pybind.cc, C++, 117 lines
third_party/nucleus/io/python/reference_wrap_test.py, Python, 99 lines
third_party/nucleus/io/python/sam_reader_pybind.cc, C++, 88 lines
third_party/nucleus/io/python/sam_reader_wrap_test.py, Python, 199 lines
third_party/nucleus/io/python/sam_writer_pybind.cc, C++, 61 lines
third_party/nucleus/io/python/tabix_indexer_pybind.cc, C++, 55 lines
third_party/nucleus/io/python/tfrecord_reader_pybind.cc, C++, 61 lines
third_party/nucleus/io/python/tfrecord_writer_pybind.cc, C++, 50 lines
third_party/nucleus/io/python/tstring_clif_conversions.h, C/C++, 61 lines
third_party/nucleus/io/python/vcf_concat_pybind.cc, C++, 53 lines
third_party/nucleus/io/python/vcf_reader_pybind.cc, C++, 95 lines
third_party/nucleus/io/python/vcf_reader_wrap_test.py, Python, 257 lines
third_party/nucleus/io/python/vcf_writer_pybind.cc, C++, 59 lines
third_party/nucleus/io/python/vcf_writer_wrap_test.py, Python, 309 lines
third_party/nucleus/io/reader_base.cc, C++, 83 lines
third_party/nucleus/io/reader_base.h, C/C++, 259 lines
third_party/nucleus/io/reader_base_test.cc, C++, 244 lines
third_party/nucleus/io/reference.cc, C++, 546 lines
third_party/nucleus/io/reference.h, C/C++, 377 lines
third_party/nucleus/io/reference_test.cc, C++, 436 lines
third_party/nucleus/io/sam.py, Python, 383 lines
third_party/nucleus/io/sam_reader.cc, C++, 1,168 lines
third_party/nucleus/io/sam_reader.h, C/C++, 216 lines
third_party/nucleus/io/sam_reader_test.cc, C++, 703 lines
third_party/nucleus/io/sam_test.py, Python, 382 lines
third_party/nucleus/io/sam_utils.cc, C++, 141 lines
third_party/nucleus/io/sam_utils.h, C/C++, 97 lines
third_party/nucleus/io/sam_utils_test.cc, C++, 70 lines
third_party/nucleus/io/sam_writer.cc, C++, 720 lines
third_party/nucleus/io/sam_writer.h, C/C++, 122 lines
third_party/nucleus/io/sam_writer_test.cc, C++, 340 lines
third_party/nucleus/io/sharded_file_utils.py, Python, 316 lines
third_party/nucleus/io/sharded_file_utils_test.py, Python, 272 lines
third_party/nucleus/io/tabix.py, Python, 45 lines
third_party/nucleus/io/tabix_indexer.cc, C++, 62 lines
third_party/nucleus/io/tabix_indexer.h, C/C++, 46 lines
third_party/nucleus/io/tabix_indexer_test.cc, C++, 104 lines
third_party/nucleus/io/tabix_test.py, Python, 103 lines
third_party/nucleus/io/text_io_test.cc, C++, 179 lines
third_party/nucleus/io/text_reader.cc, C++, 106 lines
third_party/nucleus/io/text_reader.h, C/C++, 83 lines
third_party/nucleus/io/text_writer.cc, C++, 131 lines
third_party/nucleus/io/text_writer.h, C/C++, 82 lines
third_party/nucleus/io/tfrecord.py, Python, 190 lines
third_party/nucleus/io/tfrecord_reader.cc, C++, 87 lines
third_party/nucleus/io/tfrecord_reader.h, C/C++, 89 lines
third_party/nucleus/io/tfrecord_reader_test.cc, C++, 68 lines
third_party/nucleus/io/tfrecord_test.py, Python, 124 lines
third_party/nucleus/io/tfrecord_writer.cc, C++, 103 lines
third_party/nucleus/io/tfrecord_writer.h, C/C++, 81 lines
third_party/nucleus/io/variant_reader.cc, C++, 160 lines
third_party/nucleus/io/variant_reader.h, C/C++, 158 lines
third_party/nucleus/io/variant_reader_test.cc, C++, 200 lines
third_party/nucleus/io/vcf.py, Python, 391 lines
third_party/nucleus/io/vcf_concat.cc, C++, 47 lines
third_party/nucleus/io/vcf_concat.h, C/C++, 46 lines
third_party/nucleus/io/vcf_conversion.cc, C++, 1,243 lines
third_party/nucleus/io/vcf_conversion.h, C/C++, 322 lines
third_party/nucleus/io/vcf_naive_concat.c, C, 258 lines
third_party/nucleus/io/vcf_naive_concat.h, C/C++, 35 lines
third_party/nucleus/io/vcf_reader.cc, C++, 354 lines
third_party/nucleus/io/vcf_reader.h, C/C++, 216 lines
third_party/nucleus/io/vcf_reader_test.cc, C++, 646 lines
third_party/nucleus/io/vcf_roundtrip_test.cc, C++, 78 lines
third_party/nucleus/io/vcf_test.py, Python, 484 lines
third_party/nucleus/io/vcf_writer.cc, C++, 222 lines
third_party/nucleus/io/vcf_writer.h, C/C++, 134 lines
third_party/nucleus/io/vcf_writer_test.cc, C++, 755 lines
third_party/nucleus/pip_package/build_pip_package.sh, Shell, 126 lines
third_party/nucleus/pip_package/setup.py, Python, 154 lines
third_party/nucleus/platform/types.h, C/C++, 60 lines
third_party/nucleus/protos/__init__.py, Python, 28 lines
third_party/nucleus/testing/__init__.py, Python, 28 lines
third_party/nucleus/testing/gunit_extras_test.cc, C++, 50 lines
third_party/nucleus/testing/protobuf_implementation_test.py, Python, 55 lines
third_party/nucleus/testing/protocol-buffer-matchers.cc, C++, 416 lines
third_party/nucleus/testing/protocol-buffer-matchers.h, C/C++, 1,139 lines
third_party/nucleus/testing/tensorflow_smoke_test.py, Python, 48 lines
third_party/nucleus/testing/test_utils.cc, C++, 176 lines
third_party/nucleus/testing/test_utils.h, C/C++, 183 lines
third_party/nucleus/testing/test_utils.py, Python, 401 lines
third_party/nucleus/testing/test_utils_test.py, Python, 87 lines
third_party/nucleus/tools/__init__.py, Python, 29 lines
third_party/nucleus/tools/make_cram_testdata.sh, Shell, 50 lines
third_party/nucleus/tools/md_extract_cmd.sh, Shell, 47 lines
third_party/nucleus/util/__init__.py, Python, 28 lines
third_party/nucleus/util/cigar.py, Python, 210 lines
third_party/nucleus/util/cigar_test.py, Python, 192 lines
third_party/nucleus/util/errors.py, Python, 91 lines
third_party/nucleus/util/errors_test.py, Python, 80 lines
third_party/nucleus/util/genomics_math.py, Python, 254 lines
third_party/nucleus/util/genomics_math_test.py, Python, 239 lines
third_party/nucleus/util/math.cc, C++, 100 lines
third_party/nucleus/util/math.h, C/C++, 124 lines
third_party/nucleus/util/math_test.cc, C++, 164 lines
third_party/nucleus/util/port.h, C/C++, 51 lines
third_party/nucleus/util/proto_ptr.h, C/C++, 65 lines
third_party/nucleus/util/proto_utils.py, Python, 43 lines
third_party/nucleus/util/python/__init__.py, Python, 28 lines
third_party/nucleus/util/python/math_pybind.cc, C++, 58 lines
third_party/nucleus/util/python/math_wrap_test.py, Python, 76 lines
third_party/nucleus/util/python/type_caster_nucleus_proto_ptr.h, C/C++, 146 lines
third_party/nucleus/util/python/utils_pybind.cc, C++, 53 lines
third_party/nucleus/util/python/utils_wrap_test.py, Python, 95 lines
third_party/nucleus/util/ranges.py, Python, 797 lines
third_party/nucleus/util/ranges_test.py, Python, 796 lines
third_party/nucleus/util/samplers.h, C/C++, 83 lines
third_party/nucleus/util/samplers_test.cc, C++, 68 lines
third_party/nucleus/util/sequence_utils.py, Python, 121 lines
third_party/nucleus/util/sequence_utils_test.py, Python, 107 lines
third_party/nucleus/util/struct_utils.py, Python, 239 lines
third_party/nucleus/util/struct_utils_test.py, Python, 353 lines
third_party/nucleus/util/utils.cc, C++, 314 lines
third_party/nucleus/util/utils.h, C/C++, 339 lines
third_party/nucleus/util/utils.py, Python, 124 lines
third_party/nucleus/util/utils_test.cc, C++, 568 lines
third_party/nucleus/util/utils_test.py, Python, 168 lines
third_party/nucleus/util/variant_utils.py, Python, 1,078 lines
third_party/nucleus/util/variant_utils_test.py, Python, 1,002 lines
third_party/nucleus/util/variantcall_utils.py, Python, 486 lines
third_party/nucleus/util/variantcall_utils_test.py, Python, 395 lines
third_party/nucleus/util/vcf_constants.py, Python, 518 lines
third_party/nucleus/util/vcf_constants_test.py, Python, 268 lines
third_party/nucleus/util/vis.py, Python, 854 lines
third_party/nucleus/util/vis_test.py, Python, 504 lines
tools/build_absl.sh, Shell, 100 lines
tools/preprocess_truth.py, Python, 227 lines
tools/print_f1.py, Python, 88 lines
tools/shuffle_tfrecords_beam.py, Python, 297 lines
LICENSE, License, 26 lines
README.md, Text, 263 lines

At the source: github.com/google/deepvariant

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
493 scripts, each with its path and the digest of its content;
3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

figshare:28049867, at figshare; found in “Data and code availability”

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: figshare 28049867
it says that the data are available on request
it says that the code is available on request

Read it in the paper: doi.org/10.1016/j.isci.2026.116121.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 20 authors, 5 keywords, 9 funders, 92 references.

Cite

This paper

Munkacsi, A. B., Conway, G. E., Multerer, K., Jackson, M. D., Csaki, L., Rajakumar, T., Sheridan, J. P., Niktab, E., Hwu, W.-L., Chien, Y.-H., Walterfang, M., Del Rio Hernandez, C., Chan, R. B., Zhou, B., Nandakumar, R., Zhang, P., Di Paolo, G., Maue, R. A., Reue, K., & Sturley, S. L. (2026). Multi-omics reveal phosphatidic acid phosphatases modify Niemann-Pick type C disease severity. iScience, 29(6), 116121. https://doi.org/10.1016/j.isci.2026.116121

BibTeX

@article{munkacsi2026multi, author = {Munkacsi, Andrew B and Conway, Gabriella E and Multerer, Keri and Jackson, Michael D and Csaki, Lauren and Rajakumar, Tamayanthi and Sheridan, Jeffrey P and Niktab, Eliatan and Hwu, Wuh-Liang and Chien, Yin-Hsiu and Walterfang, Mark and Del Rio Hernandez, Cintya and Chan, Robin B and Zhou, Bowen and Nandakumar, Renu and Zhang, Peixang and Di Paolo, Gilbert and Maue, Robert A and Reue, Karen and Sturley, Stephen L}, title = {{Multi-omics reveal phosphatidic acid phosphatases modify Niemann-Pick type C disease severity}}, journal = {iScience}, year = {2026}, month = may, volume = {29}, number = {6}, pages = {116121}, publisher = {Elsevier}, issn = {2589-0042}, doi = {10.1016/j.isci.2026.116121}, url = {https://doi.org/10.1016/j.isci.2026.116121}, pmid = {42256287}, pmcid = {PMC13233811} }

RIS

TY - JOUR AU - Munkacsi, Andrew B AU - Conway, Gabriella E AU - Multerer, Keri AU - Jackson, Michael D AU - Csaki, Lauren AU - Rajakumar, Tamayanthi AU - Sheridan, Jeffrey P AU - Niktab, Eliatan AU - Hwu, Wuh-Liang AU - Chien, Yin-Hsiu AU - Walterfang, Mark AU - Del Rio Hernandez, Cintya AU - Chan, Robin B AU - Zhou, Bowen AU - Nandakumar, Renu AU - Zhang, Peixang AU - Di Paolo, Gilbert AU - Maue, Robert A AU - Reue, Karen AU - Sturley, Stephen L TI - Multi-omics reveal phosphatidic acid phosphatases modify Niemann-Pick type C disease severity T2 - iScience J2 - iScience PY - 2026 DA - 2026/05/28 VL - 29 IS - 6 SP - 116121 SN - 2589-0042 PB - Elsevier DO - 10.1016/j.isci.2026.116121 UR - https://doi.org/10.1016/j.isci.2026.116121 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.isci.2026.116121", "type": "article-journal", "title": "Multi-omics reveal phosphatidic acid phosphatases modify Niemann-Pick type C disease severity", "container-title": "iScience", "author": [ { "family": "Munkacsi", "given": "Andrew B" }, { "family": "Conway", "given": "Gabriella E" }, { "family": "Multerer", "given": "Keri" }, { "family": "Jackson", "given": "Michael D" }, { "family": "Csaki", "given": "Lauren" }, { "family": "Rajakumar", "given": "Tamayanthi" }, { "family": "Sheridan", "given": "Jeffrey P" }, { "family": "Niktab", "given": "Eliatan" }, { "family": "Hwu", "given": "Wuh-Liang" }, { "family": "Chien", "given": "Yin-Hsiu" }, { "family": "Walterfang", "given": "Mark" }, { "family": "Del Rio Hernandez", "given": "Cintya" }, { "family": "Chan", "given": "Robin B" }, { "family": "Zhou", "given": "Bowen" }, { "family": "Nandakumar", "given": "Renu" }, { "family": "Zhang", "given": "Peixang" }, { "family": "Di Paolo", "given": "Gilbert" }, { "family": "Maue", "given": "Robert A" }, { "family": "Reue", "given": "Karen" }, { "family": "Sturley", "given": "Stephen L" } ], "container-title-short": "iScience", "volume": "29", "issue": "6", "page": "116121", "DOI": "10.1016/j.isci.2026.116121", "PMID": "42256287", "PMCID": "PMC13233811", "ISSN": "2589-0042", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.isci.2026.116121", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 28 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: BCFtools, JAX, SAMtools, 6 other tools, mouse
[2] doi:10.1038/s41467-026-71790-5 [code]: Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification.
Journal: Nature communications
In common: BCFtools, pysam, SAMtools, 4 other tools, genetics / omics, mouse
[3] doi:10.1007/s00262-026-04390-3 [code]: Identification and prioritisation of tumour antigen candidates from 79 glioblastoma transcriptomes.
Journal: Cancer immunology, immunotherapy : CII
In common: BCFtools, pysam, SAMtools, 3 other tools, genetics / omics, 1 reference
[4] doi:10.1186/s13059-026-04177-w [code]: Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: pysam, SAMtools, TensorFlow, 4 other tools, genetics / omics, mouse
[5] doi:10.1038/s41592-026-03057-2 [code]: CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: pysam, Keras, TensorFlow, 4 other tools, genetics / omics, mouse
[6] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pysam, Keras, TensorFlow, 4 other tools, genetics / omics
[7] doi:10.1016/j.isci.2026.116945: Neurodegeneration in the olfactory system in Niemann Pick type C1 disease.
Journal: iScience
In common: mouse, 5 references
[8] doi:10.1038/s44320-026-00208-7 [code]: Interpretable deep generative ensemble learning for single-cell omics with Hydra.
Journal: Molecular systems biology
In common: pysam, Keras, TensorFlow, 3 other tools, 1 reference
[9] doi:10.1038/s44400-026-00094-8 [code]: Haplotype-resolved DNA methylation at the <i>APOE</i> locus identifies allele-specific epigenetic signatures relevant to Alzheimer's disease risk.
Journal: NPJ dementia
In common: BCFtools, pysam, scikit-learn, 2 other tools, genetics / omics, 1 reference
[10] doi:10.21203/rs.3.rs-9927928/v1 [code]: Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: BCFtools, pysam, SAMtools, 2 other tools, genetics / omics, mouse

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.