GermlineCNVCaller – GATK

Calls copy-number variants in germline samples given their counts and the output of DetermineGermlineContigPloidy

Category Copy Number Variant Discovery

Overview

Calls copy-number variants in germline samples given their counts and the corresponding output of DetermineGermlineContigPloidy. The former should be either HDF5 or TSV count files generated by CollectReadCounts; TSV files may be compressed (e.g., with bgzip), but must then have filenames ending with the extension .gz. See the documentation for the input argument for details on enabling streaming of indexed count files from Google Cloud Storage.

Introduction

Reliable detection of copy-number variation (CNV) from read-depth ("coverage" or "counts") data such as whole exome sequencing (WES), whole genome sequencing (WGS), and custom gene sequencing panels requires a comprehensive model to account for technical variation in library preparation and sequencing. The Bayesian model and the associated inference scheme implemented in GermlineCNVCaller includes provisions for inferring and explaining away much of the technical variation. Furthermore, CNV calling confidence is automatically adjusted for each sample and genomic interval.

The parameters of the probabilistic model for read-depth bias and variance estimation (hereafter, "the coverage model") can be automatically inferred by GermlineCNVCaller by providing a cohort of germline samples sequenced using the same sequencing platform and library preparation protocol (in case of WES, the same capture kit). We refer to this run mode of the tool as the COHORT mode. The number of samples required for the COHORT mode depends on several factors such as the sequencing depth, tissue type/quality and similarity in the cohort, and the stringency of following the library preparation and sequencing protocols. For WES and WGS samples, we recommend including at least 30 samples.

The parametrized coverage model can be used for CNV calling on future case samples provided that they are strictly compatible with the cohort used to generate the model parameters (in terms of tissue type(s), library preparation and sequencing protocols). We refer to this mode as the CASE run mode. There is no lower limit on the number of samples for running GermlineCNVCaller in CASE mode.

In both tool run modes, GermlineCNVCaller requires karyotyping and global read-depth information for all samples. Such information can be automatically generated by running DetermineGermlineContigPloidy on all samples, and passed on to GermlineCNVCaller by providing the ploidy output calls using the argument contig-ploidy-calls. The ploidy state of a contig is used as the baseline ("default") copy-number state of all intervals contained in that contig for the corresponding sample. All other copy-number states are treated as alternative states and get equal shares from the total alternative state probability (set using the p-alt argument).

Note that in rare cases the inference process can diverge which will lead to an error on the Python side. In such cases the inference will be automatically restarted one more time with a new random seed. If inference diverges the second time we suggest checking if input count or karyotype values or other inputs are abnormal (an example of abnormality is a count file containing mostly zeros).

Python environment setup

The computation done by this tool, aside from input data parsing and validation, is performed outside of the Java Virtual Machine and using the gCNV computational python module, namely gcnvkernel. It is crucial that the user has properly set up a python conda environment with gcnvkernel and its dependencies installed. If the user intends to run GermlineCNVCaller using one of the official GATK Docker images, the python environment is already set up. Otherwise, the environment must be created and activated as described in the main GATK README.md file.

Advanced users may wish to set the THEANO_FLAGS environment variable to override the GATK theano configuration. For example, by running THEANO_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk GermlineCNVCaller ..., users can specify the theano compilation directory (which is set to $HOME/.theano by default). See theano documentation at https://theano-pymc.readthedocs.io/en/latest/library/config.html.

Tool run modes

COHORT mode:

The tool will be run in COHORT mode using the argument run-mode COHORT. In this mode, coverage model parameters are inferred simultaneously with the CNV states. Depending on available memory, it may be necessary to run the tool over a subset of all intervals, which can be specified by -L; this can be used to pass a filtered interval list produced by FilterIntervals to mask intervals from modeling. The specified intervals must be present in all of the input count files. The output will contain two subdirectories, one ending with "-model" and the other with "-calls".

The model subdirectory contains a snapshot of the inferred parameters of the coverage model, which may be used later for CNV calling in one or more similarly-sequenced samples as mentioned earlier. Optionally, the path to a previously obtained coverage model parametrization can be provided via the model argument in COHORT mode, in which case, the provided parameters will be only used for model initialization and a new parametrization will be generated based on the input count files. Furthermore, the genomic intervals are set to those used in creating the previous parametrization and interval-related arguments will be ignored. Note that the newly obtained parametrization ultimately reflects the input count files from the last run, regardless of whether or not an initialization parameter set is given. If the users wishes to model coverage data from two or more cohorts simultaneously, all of the input counts files must be given to the tool at once.

The calls subdirectory contains sequentially-ordered subdirectories for each sample, each listing various sample-specific estimated quantities such as the probability of various copy-number states for each interval, a parametrization of the GC curve, sample-specific unexplained variance, read-depth, and loadings of coverage bias factors.

CASE mode:

The tool will be run in CASE mode using the argument run-mode CASE. The path to a previously obtained model directory must be provided via the model argument in this mode. The modeled intervals are then specified by a file contained in the model directory, all interval-related arguments are ignored in this mode, and all model intervals must be present in all of the input count files. The tool output in CASE mode is only the "-calls" subdirectory and is organized similarly to that in COHORT mode.

Note that at the moment, this tool does not automatically verify the compatibility of the provided parametrization with the provided count files. Model compatibility may be assessed a posteriori by inspecting the magnitude of sample-specific unexplained variance of each sample, and asserting that they lie within the same range as those obtained from the cohort used to generate the parametrization. This manual step is expected to be made automatic in the future.

Important Remarks

Choice of hyperparameters:: The quality of coverage model parametrization and the sensitivity/precision of germline CNV calling are sensitive to the choice of model hyperparameters, including the prior probability of alternative copy-number states (set using the p-alt argument), prevalence of active (i.e. CNV-rich) intervals (set via the p-active argument), the coherence length of CNV events and active/silent domains across intervals (set using the cnv-coherence-length and class-coherence-length arguments, respectively), and the typical scale of interval- and sample-specific unexplained variance (set using the interval-psi-scale and sample-psi-scale arguments, respectively). It is crucial to note that these hyperparameters are not universal and must be tuned for each sequencing protocol and properly set at runtime.
Running GermlineCNVCaller on a subset of intervals:: As mentioned earlier, it may be necessary to run the tool over a subset of all intervals depending on available memory. The number of intervals must be large enough to include a genomically diverse set of regions for reliable inference of the GC bias curve, as well as other bias factors. For WES and WGS, we recommend no less than 10000 consecutive intervals spanning at least 10 - 50 mb.
Memory requirements for the python subprocess ("gcnvkernel"):: The computation done by this tool, for the most part, is performed outside of JVM and via a spawned python subprocess. The Java heap memory is only used for loading sample counts and preparing raw data for the python subprocess. The user must ensure that the machine has enough free physical memory for spawning and executing the python subprocess. Generally speaking, the resource requirements of this tool scale linearly with each of the number of samples, the number of modeled intervals, the highest copy number state, the number of bias factors, and the number of knobs on the GC curve. For example, the python subprocess requires approximately 16GB of physical memory for modeling 10000 intervals for 100 samples, with 16 maximum bias factors, maximum copy-number state of 10, and explicit GC bias modeling.

Usage examples

COHORT mode:

 gatk GermlineCNVCaller \
   --run-mode COHORT \
   -L intervals.interval_list \
   --interval-merging-rule OVERLAPPING_ONLY \
   --contig-ploidy-calls path_to_contig_ploidy_calls \
   --input normal_1.counts.hdf5 \
   --input normal_2.counts.hdf5 \
   ... \
   --output output_dir \
   --output-prefix normal_cohort_run

CASE mode:

 gatk GermlineCNVCaller \
   --run-mode CASE \
   --contig-ploidy-calls path_to_contig_ploidy_calls \
   --model previous_model_path \
   --input normal_1.counts.hdf5 \
   ... \
   --output output_dir \
   --output-prefix normal_case_run

GermlineCNVCaller specific arguments

This table summarizes the command-line arguments that are specific to this tool. For more details on each argument, see the list further down below the table or click on an argument name to jump directly to that entry in the list.

Argument name(s)	Default value	Summary
Required Arguments
--contig-ploidy-calls		Input contig-ploidy calls directory (output of DetermineGermlineContigPloidy).
--input -I		Input paths for read-count files containing integer read counts in genomic intervals for all samples. All intervals specified via -L/-XL must be contained; if none are specified, then intervals must be identical and in the same order for all samples. If read-count files are given by Google Cloud Storage paths, have the extension .counts.tsv or .counts.tsv.gz, and have been indexed by IndexFeatureFile, only the specified intervals will be queried and streamed; this can reduce disk usage by avoiding the complete localization of all read-count files.
--output -O		Output directory. This will be created if it does not exist.
--output-prefix		Prefix for output filenames.
--run-mode		Tool run-mode.
Optional Tool Arguments
--active-class-padding-hybrid-mode	50000	If copy-number-posterior-expectation-mode is set to HYBRID, CNV-active intervals determined at any time will be padded by this value (in the units of bp) in order to obtain the set of intervals on which copy number posterior expectation is performed exactly.
--adamax-beta-1	0.9	Adamax optimizer first moment estimation forgetting factor.
--adamax-beta-2	0.99	Adamax optimizer second moment estimation forgetting factor.
--annotated-intervals		Input annotated-intervals file containing annotations for GC content in genomic intervals (output of AnnotateIntervals). All intervals specified via -L must be contained. This input should not be provided if an input denoising-model directory is given (the latter already contains the annotated-interval file).
--arguments_file		read one or more arguments files and add them to the command line
--caller-external-admixing-rate	1.0	Admixing ratio of new and old called posteriors (between 0 and 1; larger values implies using more of the new posterior and less of the old posterior) after convergence.
--caller-internal-admixing-rate	0.75	Admixing ratio of new and old called posteriors (between 0 and 1; larger values implies using more of the new posterior and less of the old posterior) for internal convergence loops.
--caller-update-convergence-threshold	0.001	Maximum tolerated calling update size for convergence.
--class-coherence-length	10000.0	Coherence length of CNV-active and CNV-silent domains (in the units of bp).
--cnv-coherence-length	10000.0	Coherence length of CNV events (in the units of bp).
--convergence-snr-averaging-window	500	Averaging window for calculating training signal-to-noise ratio (SNR) for convergence checking.
--convergence-snr-countdown-window	10	The number of ADVI iterations during which the SNR is required to stay below the set threshold for convergence.
--convergence-snr-trigger-threshold	0.1	The SNR threshold to be reached before triggering the convergence countdown.
--copy-number-posterior-expectation-mode	HYBRID	The strategy for calculating copy number posterior expectations in the coverage denoising model.
--depth-correction-tau	10000.0	Precision of read depth pinning to its global value.
--disable-annealing	false	(advanced) Disable annealing.
--disable-caller	false	(advanced) Disable caller.
--disable-sampler	false	(advanced) Disable sampler.
--enable-bias-factors	true	Enable discovery of bias factors.
--gc-curve-standard-deviation	1.0	Prior standard deviation of the GC curve from flat.
--gcs-max-retries -gcs-retries	20	If the GCS bucket channel errors out, how many times it will attempt to re-initiate the connection
--gcs-project-for-requester-pays		Project to bill when accessing "requester pays" buckets. If unset, these buckets cannot be accessed. User must have storage.buckets.get permission on the bucket being accessed.
--help -h	false	display the help message
--init-ard-rel-unexplained-variance	0.1	Initial value of ARD prior precisions relative to the scale of interval-specific unexplained variance.
--initial-temperature	1.5	Initial temperature (for DA-ADVI).
--interval-merging-rule -imr	ALL	Interval merging rule for abutting intervals
--interval-psi-scale	0.001	Typical scale of interval-specific unexplained variance.
--intervals -L		One or more genomic intervals over which to operate
--learning-rate	0.01	Adamax optimizer learning rate.
--log-emission-samples-per-round	50	Log emission samples drawn per round of sampling.
--log-emission-sampling-median-rel-error	0.005	Maximum tolerated median relative error in log emission sampling.
--log-emission-sampling-rounds	10	Log emission maximum sampling rounds.
--log-mean-bias-standard-deviation	0.1	Standard deviation of log mean bias.
--mapping-error-rate	0.01	Typical mapping error rate.
--max-advi-iter-first-epoch	5000	Maximum ADVI iterations in the first epoch.
--max-advi-iter-subsequent-epochs	200	Maximum ADVI iterations in subsequent epochs.
--max-bias-factors	5	Maximum number of bias factors.
--max-calling-iters	10	Maximum number of internal self-consistency iterations within each calling step.
--max-copy-number	5	Highest allowed copy-number state.
--max-training-epochs	50	Maximum number of training epochs.
--min-training-epochs	10	Minimum number of training epochs.
--model		Input denoising-model directory. In COHORT mode, this argument is optional; if provided,a new model will be built using this input model to initialize. In CASE mode, this argument is required and the denoising model parameters are set to this input model.
--num-gc-bins	20	Number of bins used to represent the GC-bias curves.
--num-thermal-advi-iters	2500	Number of thermal ADVI iterations (for DA-ADVI).
--p-active	0.01	Prior probability of treating an interval as CNV-active (in a CNV-active domains, all copy-number states are equally likely to be called).
--p-alt	1.0E-6	Total prior probability of alternative copy-number states (the reference copy-number is set to the contig integer ploidy)
--sample-psi-scale	1.0E-4	Typical scale of sample-specific correction to the unexplained variance.
--version	false	display the version number for this tool
Optional Common Arguments
--exclude-intervals -XL		One or more genomic intervals to exclude from processing
--gatk-config-file		A configuration file to use with the GATK.
--interval-exclusion-padding -ixp	0	Amount of padding (in bp) to add to each interval you are excluding.
--interval-padding -ip	0	Amount of padding (in bp) to add to each interval you are including.
--interval-set-rule -isr	UNION	Set merging approach to use for combining interval inputs
--QUIET	false	Whether to suppress job-summary info on System.err.
--tmp-dir		Temp directory to use.
--use-jdk-deflater -jdk-deflater	false	Whether to use the JdkDeflater (as opposed to IntelDeflater)
--use-jdk-inflater -jdk-inflater	false	Whether to use the JdkInflater (as opposed to IntelInflater)
--verbosity	INFO	Control verbosity of logging.
Advanced Arguments
--showHidden	false	display hidden arguments

Argument details

Arguments in this list are specific to this tool. Keep in mind that other arguments are available that are shared with other tools (e.g. command-line GATK arguments); see Inherited arguments above.

--active-class-padding-hybrid-mode

If copy-number-posterior-expectation-mode is set to HYBRID, CNV-active intervals determined at any time will be padded by this value (in the units of bp) in order to obtain the set of intervals on which copy number posterior expectation is performed exactly.

int 50000 [ [ -∞ ∞ ] ]

--adamax-beta-1

Adamax optimizer first moment estimation forgetting factor.

double 0.9 [ [ 0 1 ] ]

--adamax-beta-2

Adamax optimizer second moment estimation forgetting factor.

double 0.99 [ [ 0 1 ] ]

--annotated-intervals

Input annotated-intervals file containing annotations for GC content in genomic intervals (output of AnnotateIntervals). All intervals specified via -L must be contained. This input should not be provided if an input denoising-model directory is given (the latter already contains the annotated-interval file).

File null

--arguments_file

read one or more arguments files and add them to the command line

List[File] []

--caller-external-admixing-rate

Admixing ratio of new and old called posteriors (between 0 and 1; larger values implies using more of the new posterior and less of the old posterior) after convergence.

double 1.0 [ [ 0 ∞ ] ]

--caller-internal-admixing-rate

Admixing ratio of new and old called posteriors (between 0 and 1; larger values implies using more of the new posterior and less of the old posterior) for internal convergence loops.

double 0.75 [ [ 0 ∞ ] ]

--caller-update-convergence-threshold

Maximum tolerated calling update size for convergence.

double 0.001 [ [ 0 ∞ ] ]

--class-coherence-length

Coherence length of CNV-active and CNV-silent domains (in the units of bp).

double 10000.0 [ [ 0 ∞ ] ]

--cnv-coherence-length

Coherence length of CNV events (in the units of bp).

double 10000.0 [ [ 0 ∞ ] ]

--contig-ploidy-calls

Input contig-ploidy calls directory (output of DetermineGermlineContigPloidy).

R File null

--convergence-snr-averaging-window

Averaging window for calculating training signal-to-noise ratio (SNR) for convergence checking.

int 500 [ [ 0 ∞ ] ]

--convergence-snr-countdown-window

The number of ADVI iterations during which the SNR is required to stay below the set threshold for convergence.

int 10 [ [ 0 ∞ ] ]

--convergence-snr-trigger-threshold

The SNR threshold to be reached before triggering the convergence countdown.

double 0.1 [ [ 0 ∞ ] ]

--copy-number-posterior-expectation-mode

The strategy for calculating copy number posterior expectations in the coverage denoising model.

The --copy-number-posterior-expectation-mode argument is an enumerated type (CopyNumberPosteriorExpectationMode), which can have one of the following values:

MAP
EXACT
HYBRID

CopyNumberPosteriorExpectationMode HYBRID

--depth-correction-tau

Precision of read depth pinning to its global value.

double 10000.0 [ [ 0 ∞ ] ]

--disable-annealing

(advanced) Disable annealing.

boolean false

--disable-caller

(advanced) Disable caller.

boolean false

--disable-sampler

(advanced) Disable sampler.

boolean false

--enable-bias-factors

Enable discovery of bias factors.

boolean true

--exclude-intervals / -XL

One or more genomic intervals to exclude from processing
Use this argument to exclude certain parts of the genome from the analysis (like -L, but the opposite). This argument can be specified multiple times. You can use samtools-style intervals either explicitly on the command line (e.g. -XL 1 or -XL 1:100-200) or by loading in a file containing a list of intervals (e.g. -XL myFile.intervals).

List[String] []

--gatk-config-file

A configuration file to use with the GATK.

String null

--gc-curve-standard-deviation

Prior standard deviation of the GC curve from flat.

double 1.0 [ [ 0 ∞ ] ]

--gcs-max-retries / -gcs-retries

If the GCS bucket channel errors out, how many times it will attempt to re-initiate the connection

int 20 [ [ -∞ ∞ ] ]

--gcs-project-for-requester-pays

Project to bill when accessing "requester pays" buckets. If unset, these buckets cannot be accessed. User must have storage.buckets.get permission on the bucket being accessed.

String ""

--help / -h

display the help message

boolean false

--init-ard-rel-unexplained-variance

Initial value of ARD prior precisions relative to the scale of interval-specific unexplained variance.

double 0.1 [ [ 0 ∞ ] ]

--initial-temperature

Initial temperature (for DA-ADVI).

double 1.5 [ [ 0 ∞ ] ]

--input / -I

Input paths for read-count files containing integer read counts in genomic intervals for all samples. All intervals specified via -L/-XL must be contained; if none are specified, then intervals must be identical and in the same order for all samples. If read-count files are given by Google Cloud Storage paths, have the extension .counts.tsv or .counts.tsv.gz, and have been indexed by IndexFeatureFile, only the specified intervals will be queried and streamed; this can reduce disk usage by avoiding the complete localization of all read-count files.

R List[String] []

--interval-exclusion-padding / -ixp

Amount of padding (in bp) to add to each interval you are excluding.
Use this to add padding to the intervals specified using -XL. For example, '-XL 1:100' with a padding value of 20 would turn into '-XL 1:80-120'. This is typically used to add padding around targets when analyzing exomes.

int 0 [ [ -∞ ∞ ] ]

--interval-merging-rule / -imr

Interval merging rule for abutting intervals
By default, the program merges abutting intervals (i.e. intervals that are directly side-by-side but do not actually overlap) into a single continuous interval. However you can change this behavior if you want them to be treated as separate intervals instead.

The --interval-merging-rule argument is an enumerated type (IntervalMergingRule), which can have one of the following values:

ALL
OVERLAPPING_ONLY

IntervalMergingRule ALL

--interval-padding / -ip

Amount of padding (in bp) to add to each interval you are including.
Use this to add padding to the intervals specified using -L. For example, '-L 1:100' with a padding value of 20 would turn into '-L 1:80-120'. This is typically used to add padding around targets when analyzing exomes.

int 0 [ [ -∞ ∞ ] ]

--interval-psi-scale

Typical scale of interval-specific unexplained variance.

double 0.001 [ [ 0 ∞ ] ]

--interval-set-rule / -isr

Set merging approach to use for combining interval inputs
By default, the program will take the UNION of all intervals specified using -L and/or -XL. However, you can change this setting for -L, for example if you want to take the INTERSECTION of the sets instead. E.g. to perform the analysis only on chromosome 1 exomes, you could specify -L exomes.intervals -L 1 --interval-set-rule INTERSECTION. However, it is not possible to modify the merging approach for intervals passed using -XL (they will always be merged using UNION). Note that if you specify both -L and -XL, the -XL interval set will be subtracted from the -L interval set.

The --interval-set-rule argument is an enumerated type (IntervalSetRule), which can have one of the following values:

UNION: Take the union of all intervals
INTERSECTION: Take the intersection of intervals (the subset that overlaps all intervals specified)

IntervalSetRule UNION

--intervals / -L

One or more genomic intervals over which to operate

List[String] []

--learning-rate

Adamax optimizer learning rate.

double 0.01 [ [ 0 ∞ ] ]

--log-emission-samples-per-round

Log emission samples drawn per round of sampling.

int 50 [ [ 0 ∞ ] ]

--log-emission-sampling-median-rel-error

Maximum tolerated median relative error in log emission sampling.

double 0.005 [ [ 0 ∞ ] ]

--log-emission-sampling-rounds

Log emission maximum sampling rounds.

int 10 [ [ 0 ∞ ] ]

--log-mean-bias-standard-deviation

Standard deviation of log mean bias.

double 0.1 [ [ 0 ∞ ] ]

--mapping-error-rate

Typical mapping error rate.

double 0.01 [ [ 0 1 ] ]

--max-advi-iter-first-epoch

Maximum ADVI iterations in the first epoch.

int 5000 [ [ 0 ∞ ] ]

--max-advi-iter-subsequent-epochs

Maximum ADVI iterations in subsequent epochs.

int 200 [ [ 0 ∞ ] ]

--max-bias-factors

Maximum number of bias factors.

int 5 [ [ 0 ∞ ] ]

--max-calling-iters

Maximum number of internal self-consistency iterations within each calling step.

int 10 [ [ 0 ∞ ] ]

--max-copy-number

Highest allowed copy-number state.

int 5 [ [ 0 ∞ ] ]

--max-training-epochs

Maximum number of training epochs.

int 50 [ [ 0 ∞ ] ]

--min-training-epochs

Minimum number of training epochs.

int 10 [ [ 0 ∞ ] ]

--model

Input denoising-model directory. In COHORT mode, this argument is optional; if provided,a new model will be built using this input model to initialize. In CASE mode, this argument is required and the denoising model parameters are set to this input model.

File null

--num-gc-bins

Number of bins used to represent the GC-bias curves.

int 20 [ [ 2 ∞ ] ]

--num-thermal-advi-iters

Number of thermal ADVI iterations (for DA-ADVI).

int 2500 [ [ 0 ∞ ] ]

--output / -O

Output directory. This will be created if it does not exist.

R File null

COHORT
CASE

R RunMode null

--sample-psi-scale

Typical scale of sample-specific correction to the unexplained variance.

double 1.0E-4 [ [ 0 ∞ ] ]

--showHidden / -showHidden

display hidden arguments

boolean false

--tmp-dir

Temp directory to use.

GATKPath null

--use-jdk-deflater / -jdk-deflater

Whether to use the JdkDeflater (as opposed to IntelDeflater)

boolean false

--use-jdk-inflater / -jdk-inflater

Whether to use the JdkInflater (as opposed to IntelInflater)

boolean false

--verbosity / -verbosity

Control verbosity of logging.

The --verbosity argument is an enumerated type (LogLevel), which can have one of the following values:

ERROR
WARNING
INFO
DEBUG

LogLevel INFO

--version

display the version number for this tool

boolean false

Return to top

GATK version 4.2.2.0-SNAPSHOT built at Thu, 19 Aug 2021 09:49:28 -0700.

Genome Analysis Toolkit

Need Help?

Community Forum

Articles in this section

Category Copy Number Variant Discovery

Overview

Introduction

Python environment setup

Tool run modes

Important Remarks

Usage examples

GermlineCNVCaller specific arguments

Argument details

--exclude-intervals / -XL

--gcs-max-retries / -gcs-retries

--help / -h

--input / -I

--interval-exclusion-padding / -ixp

--interval-merging-rule / -imr

--interval-padding / -ip

--interval-set-rule / -isr

--intervals / -L

--output / -O

--showHidden / -showHidden

--use-jdk-deflater / -jdk-deflater

--use-jdk-inflater / -jdk-inflater

--verbosity / -verbosity

Related articles